Energy Stable Discontinuous Galerkin Methods for Maxwell's Equations in Nonlinear Optical Media

TL;DR

This paper develops energy-stable high-order discontinuous Galerkin methods for Maxwell's equations in nonlinear optical media, incorporating nonlinear responses and dispersion, with proven stability and demonstrated effectiveness through numerical simulations.

Contribution

It introduces energy-stable DG discretizations for nonlinear Maxwell's equations with novel strategies for nonlinear temporal discretization.

Findings

Energy stability of semi-discrete methods is established.

Error estimates are provided under certain nonlinearity restrictions.

Numerical simulations demonstrate the methods' effectiveness in modeling nonlinear wave phenomena.

Abstract

The propagation of electromagnetic waves in general media is modeled by the time-dependent Maxwell's partial differential equations (PDEs), coupled with constitutive laws that describe the response of the media. In this work, we focus on nonlinear optical media whose response is modeled by a system of first order nonlinear ordinary differential equations (ODEs), which include a single resonance linear Lorentz dispersion, and the nonlinearity comes from the instantaneous electronic Kerr response and the residual Raman molecular vibrational response. To design efficient, accurate, and stable computational methods, we apply high order discontinuous Galerkin discretizations in space to the hybrid PDE-ODE Maxwell system with several choices of numerical fluxes, and the resulting semi-discrete methods are shown to be energy stable. Under some restrictions on the strength of the nonlinearity, error estimates are also established. When we turn to fully discrete methods, the challenge to achieve provable stability lies in the temporal discretizations of the nonlinear terms. To overcome this, novel strategies are proposed to treat the nonlinearity in our model within the framework of the second-order leap-frog and implicit trapezoidal time integrators. The performance of the overall algorithms are demonstrated through numerical simulations of kink and antikink waves, and third-harmonic generation in soliton propagation.

Figures

Tables

Table 1. Table 1: A traveling kink and antikink wave: C ​ F ​ L 𝐶 𝐹 𝐿 CFL number.

$k$	Leap-frog scheme	Fully implicit scheme
1	$0.2/v$	5
2	1	10
3	2	20

Table 2. Table 2: A traveling kink and antikink wave: errors and orders of accuracy of E 𝐸 E . k = 1 𝑘 1 k=1 .

	$N$	Leap-frog scheme				Fully implicit scheme
		$L_2$ errors	order	$L_{\infty }$ error	order	$L_2$ errors	order	$L_{\infty }$ error	order
upwind flux	100	4.48E-04	–	1.80E-03	–	9.24E-03	–	2.38E-02	–
	200	7.68E-05	2.55	3.46E-04	2.38	3.75E-03	1.30	1.16E-02	1.04
	400	1.32E-05	2.54	5.92E-05	2.55	1.15E-03	1.71	4.20E-03	1.46
	800	2.75E-06	2.26	1.17E-05	2.33	3.01E-04	1.93	1.20E-03	1.80
	1600	6.57E-07	2.06	3.10E-06	1.92	7.59E-05	1.99	3.11E-04	1.95
central flux	100	1.07E-03	–	4.53E-03	–	9.13E-03	–	2.36E-02	–
	200	2.77E-04	1.94	1.31E-03	1.79	3.64E-03	1.33	1.18E-02	1.12
	400	7.12E-05	1.96	3.59E-04	1.87	1.10E-03	1.72	4.20E-03	1.49
	800	1.97E-05	1.85	1.05E-04	1.78	2.90E-04	1.93	1.23E-03	1.77
	1600	6.91E-06	1.51	3.51E-05	1.57	7.32E-05	1.98	3.24E-04	1.93
alternating flux I	100	1.15E-04	–	5.22E-03	–	9.38E-03	–	2.41E-02	–
	200	2.96E-05	1.96	1.13E-04	2.21	3.77E-03	1.31	1.16E-02	1.06
	400	8.91E-06	1.73	3.45E-05	1.71	1.15E-03	1.71	4.21E-03	1.46
	800	2.01E-06	2.15	9.48E-06	1.86	3.01E-04	1.93	1.21E-03	1.80
	1600	4.62E-07	2.12	2.03E-06	2.22	7.59E-05	1.99	3.11E-04	1.95
alternating flux II	100	9.41E-05	–	4.79E-04	–	9.39E-03	–	2.41E-02	–
	200	3.04E-05	1.63	1.49E-04	1.69	3.78E-03	1.31	1.16E-02	1.05
	400	8.36E-06	1.86	3.61E-05	2.04	1.15E-03	1.71	4.20E-03	1.47
	800	1.87E-06	2.16	8.79E-06	2.04	3.01E-04	1.93	1.21E-03	1.80
	1600	5.40E-07	1.79	2.65E-06	1.73	7.59E-05	1.99	3.11E-04	1.95

Table 3. Table 3: A traveling kink and antikink wave: errors and orders of accuracy of E 𝐸 E . k = 2 𝑘 2 k=2 .

	$N$	Leap-frog scheme				Fully implicit scheme
		$L_2$ errors	order	$L_{\infty }$ error	order	$L_2$ errors	order	$L_{\infty }$ error	order
upwind flux	100	2.49E-05	–	1.50E-04	–	3.66E-03	–	1.13E-02	–
	200	3.04E-06	3.04	1.88E-05	3.00	5.71E-04	2.68	2.21E-03	2.35
	400	3.69E-07	3.04	2.28E-06	3.04	7.29E-05	2.97	2.99E-04	2.89
	800	4.72E-08	2.96	3.04E-07	2.91	9.13E-06	3.00	3.77E-05	2.99
central flux	100	2.32E-05	–	1.05E-04	–	3.66E-03	–	1.13E-02	–
	200	2.97E-06	2.96	1.38E-05	2.92	5.72E-04	2.68	2.21E-03	2.35
	400	3.71E-07	3.00	1.75E-06	2.99	7.29E-05	2.97	2.99E-04	2.89
	800	4.62E-08	3.00	2.19E-07	2.99	9.13E-06	3.00	3.77E-05	2.99
alternating flux I	100	2.38E-05	–	1.04E-04	–	3.66E-03	–	1.13E-02	–
	200	2.98E-06	3.00	1.33E-05	2.96	5.71E-04	2.68	2.21E-03	2.35
	400	3.72E-07	3.00	1.70E-06	2.96	7.29E-05	2.97	2.99E-04	2.89
	800	4.62E-08	3.01	2.08E-07	3.04	9.13E-06	3.00	3.77E-05	2.99
alternating flux II	100	2.41E-05	–	1.25E-04	–	3.66E-03	–	1.13E-02	–
	200	3.03E-06	2.99	1.67E-05	2.91	5.72E-04	2.68	2.21E-03	2.35
	400	3.73E-07	3.02	1.98E-06	3.07	7.29E-05	2.97	2.99E-04	2.89
	800	4.66E-08	3.00	2.67E-07	2.89	9.13E-06	3.00	3.77E-05	2.99

Table 4. Table 4: A traveling kink and antikink wave: errors and orders of accuracy of E 𝐸 E . k = 3 𝑘 3 k=3 .

	$N$	Leap-frog scheme				Fully implicit scheme
		$L_2$ errors	order	$L_{\infty }$ error	order	$L_2$ errors	order	$L_{\infty }$ error	order
upwind flux	100	5.58E-06	–	2.41E-05	–	1.07E-03	–	3.93E-03	–
	200	3.50E-07	3.99	1.48E-06	4.03	7.00E-05	3.93	2.87E-04	3.77
	400	2.18E-08	4.00	1.05E-07	3.82	4.38E-06	4.00	1.81E-05	3.99
central flux	100	5.68E-06	–	2.59E-05	–	1.07E-03	–	3.93E-03	–
	200	3.62E-07	3.97	1.73E-06	3.91	7.00E-05	3.93	2.87E-04	3.77
	400	2.27E-08	4.00	1.09E-07	3.98	4.37E-06	4.00	1.80E-05	3.99
alternating flux I	100	5.60E-06	–	2.30E-05	–	1.07E-03	–	3.93E-03	–
	200	3.53E-07	3.99	1.44E-06	4.00	7.00E-05	3.93	2.87E-04	3.77
	400	2.19E-08	4.01	8.97E-08	4.00	4.38E-06	4.00	1.80E-05	3.99
alternating flux II	100	5.60E-06	–	2.35E-05	–	1.07E-03	–	3.93E-03	–
	200	3.53E-07	3.99	1.45E-06	4.02	7.00E-05	3.93	2.87E-04	3.77
	400	2.19E-08	4.01	9.01E-08	4.01	4.38E-06	4.00	1.81E-05	3.99

Table 5. Table 5: Soliton propagation: C ​ F ​ L 𝐶 𝐹 𝐿 CFL number.

Leap-frog scheme		Fully implicit scheme
Central/upwind flux	Alternating flux	Central/upwind flux	Alternating flux
0.05	0.1	0.3	0.5

Equations

$${\partial }_t\mathbf{B}+\nabla \times \mathbf{E}=0,\text{in}(0,T)\times \Omega ,$$

$${\partial }_t\mathbf{D}+{\mathbf{J}}_{\mathbf{s}}-\nabla \times \mathbf{H}=0,\text{in}(0,T)\times \Omega ,$$

$$\nabla \cdot \mathbf{B}=0,\nabla \cdot \mathbf{D}=\rho ,\text{in}(0,T)\times \Omega ,$$

$$\mathbf{D}={ϵ}_0({ϵ}_{\infty }\mathbf{E}+\mathbf{P}),\mathbf{B}={\mu }_0\mathbf{H},$$

$$\frac{{\partial }^2{\mathbf{P}}_{\mathrm{delay}}^{\mathrm{L}}}{\partial t^2}+\frac{1}{\tau }\frac{\partial {\mathbf{P}}_{\mathrm{delay}}^{\mathrm{L}}}{\partial t}+{\omega }_0^2{\mathbf{P}}_{\mathrm{delay}}^{\mathrm{L}}={\omega }_p^2\mathbf{E}.$$

$${\mathbf{P}}^{\mathrm{NL}}={\mathbf{P}}_{\mathrm{Kerr}}^{\mathrm{NL}}+{\mathbf{P}}_{\mathrm{delay}}^{\mathrm{NL}}=\underset{\mathrm{Kerr}}{\underbrace{a(1-\theta )\mathbf{E}|\mathbf{E}{|}^2}}+\underset{\mathrm{Raman}}{\underbrace{a\theta Q\mathbf{E}}}.$$

$$\frac{{\partial }^{2}Q}{\partial t^2}+\frac{1}{{\tau }_v}\frac{\partial Q}{\partial t}+{\omega }_v^{2}Q={\omega }_v^2|\mathbf{E}{|}^2,$$

$$\mathbf{D}={ϵ}_0({ϵ}_{\infty }\mathbf{E}+{\mathbf{P}}_{\mathrm{delay}}^{\mathrm{L}}+a(1-\theta )\mathbf{E}|\mathbf{E}{|}^2+a\theta Q\mathbf{E}).$$

$${\mu }_0{\partial }_{t}H$$

$${\partial }_{t}D$$

$${\partial }_{t}P$$

$${\partial }_{t}J$$

$${\partial }_{t}Q$$

$${\partial }_t\sigma$$

$$D={ϵ}_0({ϵ}_{\infty }E+P+a(1-\theta )E^3+a\theta QE),$$

$$\mathcal{E}={\int }_{\Omega }\left(\frac{{\mu }_0}{2}{H}^2+\frac{{ϵ}_0{ϵ}_{\infty }}{2}{E}^2+\frac{{ϵ}_0}{2{\omega }_p^2}{J}^2+\frac{{ϵ}_0{\omega }_0^2}{2{\omega }_p^2}{P}^2+\frac{{ϵ}_{0}a\theta }{4{\omega }_v^2}{\sigma }^2+\frac{{ϵ}_{0}a\theta }{2}Q{E}^2+\frac{3{ϵ}_{0}a(1-\theta )}{4}{E}^4+\frac{{ϵ}_{0}a\theta }{4}{Q}^2\right)dx,$$

$$\frac{d}{dt}\mathcal{E}=-\frac{{ϵ}_0}{{\omega }_p^2\tau }{\int }_{\Omega }{J}^{2}dx-\frac{{ϵ}_{0}a\theta }{2{\omega }_v^2{\tau }_v}{\int }_{\Omega }{\sigma }^{2}dx\le 0.$$

$$V_h^k=\{v:v{|}_{I_j}\in P^k(I_j),j=1,2,\cdots ,N\},$$

$${\mu }_0$$

$${\int }_{I_j}{\partial }_t{D}_h\varphi dx+{\int }_{I_j}{H}_h{\partial }_x\varphi dx-(\widetilde{H_h}{\varphi }^{-}{)}_{j+1/2}+(\widetilde{H_h}{\varphi }^{+}{)}_{j-1/2}=0,\forall \varphi \in V_h^k,$$

$${\partial }_t{P}_h=J_h,$$

$${\partial }_t{J}_h=-\left(\frac{1}{\tau }{J}_h+{\omega }_0^2{P}_h-{\omega }_p^2{E}_h\right),$$

$${\partial }_t{Q}_h={\sigma }_h,$$

$${\int }_{I_j}{\partial }_t{\sigma }_h\varphi dx=-{\int }_{I_j}\left(\frac{1}{{\tau }_v}{\sigma }_h+{\omega }_v^2{Q}_h-{\omega }_v^2{E}_h^2\right)\varphi dx,\forall \varphi \in V_h^k.$$

$${\int }_{I_j}{D}_h\varphi dx={\int }_{I_j}{ϵ}_0\left({ϵ}_{\infty }{E}_h+a(1-\theta )E_h^3+P_h+a\theta Q_h{E}_h\right)\varphi dx,\forall \varphi \in V_h^k.$$

$$\widehat{E_h}=\{E_h\},\widetilde{H_h}=\{H_h\},$$

$$\widehat{E_h}=E_h^{-},\widetilde{H_h}=H_h^{+};\widehat{E_h}=E_h^{+},\widetilde{H_h}=H_h^{-},$$

$$\widehat{E_h}=\{E_h\}+\frac{1}{2}\sqrt{\frac{{\mu }_0}{{ϵ}_0{ϵ}_{\infty }}}[H_h],\widetilde{H_h}=\{H_h\}+\frac{1}{2}\sqrt{\frac{{ϵ}_0{ϵ}_{\infty }}{{\mu }_0}}[E_h].$$

$$\frac{d}{dt}{\mathcal{E}}_h=-\frac{{ϵ}_0}{{\omega }_p^2\tau }{\int }_{\Omega }{J}_h^{2}dx-\frac{{ϵ}_{0}a\theta }{2{\omega }_v^2{\tau }_v}{\int }_{\Omega }{\sigma }_h^{2}dx\le 0,$$

$$\frac{d}{dt}{\mathcal{E}}_h=-\frac{{ϵ}_0}{{\omega }_p^2\tau }{\int }_{\Omega }{J}_h^{2}dx-\frac{{ϵ}_{0}a\theta }{2{\omega }_v^2{\tau }_v}{\int }_{\Omega }{\sigma }_h^{2}dx-\frac{1}{2}\sqrt{\frac{{\mu }_0}{{ϵ}_0{ϵ}_{\infty }}}\sum _{j=1}^N[H_h{]}_{j+1/2}^2-\frac{1}{2}\sqrt{\frac{{ϵ}_0{ϵ}_{\infty }}{{\mu }_0}}\sum _{j=1}^N[E_h{]}_{j+1/2}^2\le 0,$$

$${\mathcal{E}}_h={\int }_{\Omega }\left(\frac{{\mu }_0}{2}{H}_h^2+\frac{{ϵ}_0{ϵ}_{\infty }}{2}{E}_h^2+\frac{{ϵ}_0}{2{\omega }_p^2}{J}_h^2+\frac{{ϵ}_0{\omega }_0^2}{2{\omega }_p^2}{P}_h^2+\frac{{ϵ}_{0}a\theta }{4{\omega }_v^2}{\sigma }_h^2+\frac{{ϵ}_{0}a\theta }{2}{Q}_h{E}_h^2+\frac{3{ϵ}_{0}a(1-\theta )}{4}{E}_h^4+\frac{{ϵ}_{0}a\theta }{4}{Q}_h^2\right)dx$$

$${\int }_{\Omega }({\mu }_0{H}_h{\partial }_t{H}_h+E_h{\partial }_t{D}_h)dx+\sum _{j=1}^N{\int }_{I_j}{\partial }_x(E_h{H}_h)dx+\sum _{j=1}^N(\widehat{E_h}[H_h]+\widetilde{H_h}[E_h]{)}_{j-1/2}=0.$$

$$\widehat{E_h}[H_h]+\widetilde{H_h}[E_h]=[E_h{H}_h],$$

$$\widehat{E_h}[H_h]+\widetilde{H_h}[E_h]=[E_h{H}_h]+\frac{1}{2}\sqrt{\frac{{\mu }_0}{{ϵ}_0{ϵ}_{\infty }}}[H_h{]}^2+\frac{1}{2}\sqrt{\frac{{ϵ}_0{ϵ}_{\infty }}{{\mu }_0}}[E_h{]}^2,$$

$${\int }_{\Omega }({\mu }_0{H}_h{\partial }_t{H}_h+E_h{\partial }_t{D}_h)dx=M(E_h,H_h),$$

$$M(E_{h},H_{h}):=\left\{\begin{array}[]{ll}0,&\mbox{for flux choices \eqref{eq:flux:c} and \eqref{eq:flux:a}};\\ -\frac{1}{2}\sqrt{\frac{\mu_{0}}{\epsilon_{0}\epsilon_{\infty}}}\sum_{j=1}^{N}[H_{h}]_{j-1/2}^{2}-\frac{1}{2}\sqrt{\frac{\epsilon_{0}\epsilon_{\infty}}{\mu_{0}}}\sum_{j=1}^{N}[E_{h}]_{j-1/2}^{2},&\mbox{for flux choice \eqref{eq:flux:u} },\end{array}\right.$$

Full text

1 Introduction

Nonlinear optics is the study of the behavior of light in nonlinear media. This field has developed into a significant branch of physics since the introduction of intense lasers with high peak powers. In nonlinear media, the material response depends nonlinearly on the optical field, and many interesting physical phenomena, such as frequency mixing and second/third-harmonic generation have been observed and harnessed for practical applications. We refer to classical textbooks [3, 5, 37] for a more detailed review of the field of nonlinear optics.

Our interest here is in the development of novel numerical schemes for the Maxwell’s equations in nonlinear optical media. Relative to the widely used asymptotic and paraxial wave models derived from Maxwell’s equations, such as nonlinear Schrödinger equation (NLS) and beam propagation method (BPM) [3, 5], simulations of the nonlinear Maxwell’s system in the time domain are more computationally intensive. However, these simulations have the advantage of being substantially more robust because they directly solve for fundamental quantities, the electromagnetic fields in space and time. These simulations also avoid the simplifying assumptions that lead to conventional asymptotic and paraxial propagation analyses, and are able to treat interacting waves at different frequencies directly [28]. Recent optics and photonics research has focused on phenomena at smaller and smaller length scales or multiple spatial scales. For such phenomena simulating the full Maxwell PDE models is important in adequately capturing useful optical effects [20, 28, 41].

When Maxwell’s equations are considered to model the electromagnetic (EM) waves propagating through a nonlinear optical medium, the medium response is described by constitutive laws that relate the electric field ${\mathbf{E}}$ and the electric flux density ${\mathbf{D}}$ through the polarization ${\mathbf{P}}$ of the medium. In this work, we focus on a macroscopic phenomenological description of the polarization, which comprises both linear and nonlinear responses. Specifically, the linear response is modeled by a single resonance Lorentz dispersion, while the nonlinear response is cubic and incorporates the instantaneous Kerr effect and the delayed nonlinear Lorentz dispersion called Raman scattering. Within this description, we will follow the auxiliary differential equation (ADE) approach, where the linear and nonlinear Lorentz dispersion is represented through a set of ODEs, describing the time evolution of ${\mathbf{P}}$ (hence of ${\mathbf{D}}$) forced by ${\mathbf{E}}$, appended to Maxwell’s equations. An alternative representation is via a recursive convolution method, where ${\mathbf{D}}$ is computed from ${\mathbf{E}}$ through a time convolution integral [40].

In the literature, finite difference time domain (FDTD) based, finite element (FEM) based, pseudospectral based methods, finite volume (FV) based, among others, are available for the integration of the full Maxwell’s equations in nonlinear media, along with appended ODEs for the material response. The Yee scheme [43] is an FDTD method for Maxwell’s equations that has long been one of the gold standards for numerical simulation of Maxwell’s equations in the time domain, especially for linear problems [40]. Maxwell’s equations in a linear Lorentz medium with a nonlinear Kerr response are investigated in [25, 39], while in [19], additional effects due to Raman scattering are studied through a 1D FDTD analysis. More references for linear and nonlinear Lorentz dispersion, can be found in [4, 22, 38] for the 1D case, and in [17, 29, 44] for 2D and 3D cases. Yee based FDTD approaches result in second order schemes which accumulate significant errors over long time modeling of wave propagation [12, 13]. While higher order FDTD methods can alleviate this issue, they can be cumbersome in modeling complex geometries. On the other hand, though the FEMs are well suited for modeling complex geometries, they have not been well developed for nonlinear Maxwell models. FEM analysis for some nonlinear models can be found in [16]. In [42] a pseudospectral spatial domain (PSSD) approach is presented for linear Lorentz dispersion and nonlinear Kerr response, and in [30] optical carrier wave shock is studied using the PSSD technique. FV based methods for nonlinear Kerr media are addressed in [1, 14] in which the Maxwell-Kerr model is approached as a hyperbolic system and approximated by a Godunov scheme, and a third order Roe solver, respectively, in one and two spatial dimensions.

In this work, we use high order discontinuous Galerkin (DG) methods for the spatial discretization of our nonlinear Maxwell models. This is motivated by various properties of DG methods, including high order accuracy, excellent dispersive and dissipative properties in standard wave simulations, flexibility in adaptive implementation and high parallelization, and suitability for complicated geometry (e.g., [10, 24]). DG methods differ from classical finite element methods in their use of piecewise smooth approximate functions, while inter-element communication is achieved through the use of numerical fluxes, which are consistent with the physical fluxes and play a vital role in accuracy, stability, energy conservation, and computational efficiency. For the nonlinear Maxwell models that we consider, with the numerical fluxes chosen to be either central or alternating, the solutions to the semi-discrete DG methods satisfy an energy equation just as the exact solutions do, hence the methods are energy stable, even in the presence of both the Kerr and Raman nonlinear effects. Another dissipative flux, inspired by the upwind flux for Maxwell’s equations in a linear nondispersive dielectric [23], called “upwind flux” in this paper, is also considered with the respective energy stability established. For the semi-discrete methods with all three types of numerical fluxes, error estimates are carried out under some additional assumptions on the strength of the nonlinearity in the underlying model.

In addition to the error estimates, the nonlinearity in the model poses challenges to the design of fully discrete schemes with provable energy stability. As one major contribution, we propose in this work a novel strategy to discretize the nonlinear terms within the commonly used second-order leap-frog and implicit trapezoidal temporal discretizations. The resulting fully discrete methods are proved to be stable. More specifically, the method with the modified leap-frog time discretization is conditionally stable under a CFL condition, which is the same as the one for Maxwell’s equations without Kerr, linear and nonlinear Lorentz dispersion; while the fully implicit method with the modified trapezoidal temporal discretization is unconditionally stable. In both cases, we find it important, at least from the theoretical point of view, to discretize the ODE part of the system implicitly. To our best knowledge, the temporal discretizations that are adapted to nonlinear models and with provable stability are not yet available. In the present work, the methods and numerical verification are presented for the model in one dimension, and their extension to higher dimension will be explored in a separate paper.

DG methods have grown to be broadly adopted for EM simulations in the past two decades. They have been developed and analyzed for time dependent linear models, including Maxwell’s equations in free space (e.g., [8, 11, 23]), dispersive media (e.g., [18, 27, 32, 36]), as well as metamaterials (e.g., [7, 33, 34, 35]). However, there exists only limited study for DG methods for nonlinear Maxwell models. For example, in [2, 15], Kerr nonlinearity is investigated, where the entire Maxwell PDE-ODE system is cast as a nonlinear hyperbolic conservation law, for which DG methods have long been known for their success. A relaxed version of the Kerr model, called the Kerr-Debye model, was examined in [26], where a second-order asymptotic-preserving and positivity-preserving DG scheme is designed and analyzed.

The rest of this paper is organized as follows. In Section 2, Maxwell’s equations in an optical medium with a nonlinear dispersive response are introduced. In Section 3, DG spatial discretizations are formulated, where energy stability is established for the resulting semi-discrete schemes. Error estimates are further carried out in Section 4. In Section 5, temporal discretizations are presented within the framework of the second-order leap-frog and trapezoidal method, with a novel treatment of the nonlinear terms in the models aimed at obtaining energy stability for the fully discrete schemes. The performance of the overall algorithms are demonstrated in Section 6 through numerical simulations of the propagating kink and antikink waves, and third harmonic generation in soliton propagation. Finally, concluding remarks are made in Section 7.

2 Physical Model: Maxwell’s Equations and Polarization

We begin with the Maxwell’s equations, that govern the time evolution of the electric field $\mathbf{E}$ and magnetic field $\mathbf{H}$ in a non-magnetic nonlinear optical medium,

[TABLE]

along with initial and boundary data in the domain $\Omega\subset\mathbb{R}^{d},d=1,2,3$. The variable $\mathbf{J}_{s}$ is the source current density, and $\rho$ is the charge density. The electric flux density $\mathbf{D}$ and the magnetic induction $\mathbf{B}$ are related to the electric and magnetic field, respectively, via the constitutive laws

[TABLE]

where $\mathbf{P}$ is the polarization. The dielectric parameters are $\epsilon_{0}$, the electric permittivity of free space, $\epsilon_{\infty}$, the relative electric permittivity in the limit of the infinite frequency, and $\mu_{0}$, the magnetic permeability of free space. We will assume here that all model parameters are constant, and the material is isotropic. The term $\epsilon_{0}\epsilon_{\infty}{\mathbf{E}}$ captures the linear instantaneous response of the material to the EM fields.

To model the linear and nonlinear dispersion in the material we use the auxiliary differential equation (ADE) approach as presented in [19, 40]. The linear (L) delayed or retarded response of the material to the EM field is captured in the polarization, $\mathbf{P}$, via a linear single resonance Lorentz response, which, in the form of a second order ODE, is given as,

[TABLE]

Here $\omega_{0}$ and $\omega_{p}$ are the resonance and plasma frequencies of the medium, respectively, and $\tau^{-1}$ is a damping constant. In addition, $\omega_{p}^{2}=(\epsilon_{s}-\epsilon_{\infty})\omega_{0}^{2}$, with $\epsilon_{s}$ as the relative permittivity at zero frequency.

For pulse widths that are sufficiently short (for e.g., shorter than 1 pico-second (ps) for Silica) [25], the nonlinear response has an instantaneous as well as a delayed component. For the nonlinear (NL) response of the medium, we will consider a cubic Kerr-type instantaneous response, and a retarded Raman molecular vibrational response called Raman scattering. The Kerr effect is a phenomenon in which the refractive index of a material changes proportionally to the square of the applied electric field. Raman scattering arises from the electric field induced changes in the internal nuclear vibrations on time scales $\approx$ 1 to 100 femto-seconds (fs) [21], and is modeled by a nonlinear single resonance Lorentz delayed response. The two nonlinear responses are given as

[TABLE]

Here $a$ is a third order coupling constant, $\theta$ parameterizes the relative strength of the instantaneous electronic Kerr and retarded Raman molecular vibrational responses, and $Q$ describes the natural molecular vibrations within the dielectric material that has frequency many orders of magnitude less than the optical wave frequency, responding to the field intensity. The time evolution of $Q$ is given by the following ODE,

[TABLE]

where $\omega_{v}$ is the resonance frequency of the vibration, and $\tau_{v}^{-1}$ a damping constant. This is essentially a model for a simple linear oscillator, but coupled to the nonlinear field intensity $|{\mathbf{E}}|^{2}$.

Taking into account all the effects discussed above, the constitutive law for the electric flux density is given by

[TABLE]

With this, the mathematical model for EM wave propagation in this nonlinear optical medium will be given as a PDE-ODE system (1)-(5).

In the present work, the first order form of the second order ODEs, (3) and (4), will be adopted, and we will focus our investigation on the model in one spatial dimension, as below,

[TABLE]

with the constitutive law

[TABLE]

where $P=P^{L}_{\mathrm{delay}}$. In this model, we assume uniformity of all the vector fields in the $y$ and $z$ directions. Thus, all derivatives with respect to $y$ and $z$ in the curl and divergence operators are set to zero. All field quantities are represented by a single scalar component. The scalar magnetic field $H$ (hence $B$) represents the 2nd (or the 3rd) component of the vector magnetic field ${\mathbf{H}}$, and the scalar electric flux density $D$ (hence $E$) represents the 3rd (or the 2nd) component of ${\mathbf{D}}$ (hence ${\mathbf{E}}$). Gauss’s laws (1c) only involve the $x$ derivatives of the 1st components of ${\mathbf{B}}$ and ${\mathbf{D}}$, and therefore they are decoupled from the one-dimensional model (6)-(7) and become irrelevant. Under the assumption of periodic boundary conditions, the energy $\mathcal{E}=\mathcal{E}(t)$ of the system (6), defined as

[TABLE]

satisfies the following relation,

[TABLE]

Note that $\mathcal{E}(t)$ is guaranteed non-negative only when $\theta\in[0,\frac{3}{4}]$.

3 Semi-discrete Scheme: Discontinuous Galerkin Method

In this section, we introduce a semi-discrete DG method in space for the one dimensional model problem (6) - (7). For simplicity, periodic boundary conditions are considered in $x$ direction. (See Sect. 6.2 and the appendix for some more general boundary conditions.) Let $\Omega=[x_{L},x_{R}]$ be the computational domain, for which a mesh, $x_{L}=x_{1/2}<x_{3/2}<\cdots<x_{N+1/2}=x_{R}$, is introduced. Let $I_{j}=[x_{j-1/2},x_{j+1/2}]$ be a mesh element, with $x_{j}=\frac{1}{2}(x_{j-\frac{1}{2}}+x_{j+\frac{1}{2}})$ as its center, $h_{j}=x_{j+\frac{1}{2}}-x_{j-\frac{1}{2}}$ as its length, and $h=\max_{1\leq j\leq N}h_{j}$ as the largest meshsize. We now define a finite dimensional discrete space,

[TABLE]

which consists of piecewise polynomials of degree up to $k$ with respect to the mesh. For any $v\in V_{h}^{k}$, let $v^{+}_{j+\frac{1}{2}}$ (resp. $v^{-}_{j+\frac{1}{2}}$) denote the limit value of $v$ at $x_{j+\frac{1}{2}}$ from the element $I_{j+1}$ (resp. $I_{j}$), $[v]_{j+\frac{1}{2}}=v^{+}_{j+\frac{1}{2}}-v^{-}_{j+\frac{1}{2}}$ denote its jump, and $\{v\}_{j+\frac{1}{2}}=\frac{1}{2}(v^{+}_{j+\frac{1}{2}}+v^{-}_{j+\frac{1}{2}})$ be its average, again at $x_{j+\frac{1}{2}}$. The mesh is assumed to be quasi-uniform, namely, there exists a positive constant $\delta$, such that $\frac{h}{\min_{j}h_{j}}<\delta$, as the mesh is refined.

The semi-discrete DG method for the system (6) - (7) is formulated as follows: find $H_{h}(t,\cdot)$, $D_{h}(t,\cdot)$, $E_{h}(t,\cdot)$, $P_{h}(t,\cdot)$, $J_{h}(t,\cdot)$, $Q_{h}(t,\cdot)$, $\sigma_{h}(t,\cdot)\in V_{h}^{k}$, such that $\forall j,$

[TABLE]

The constitutive law is imposed via the $L^{2}$ projection, namely,

[TABLE]

Both terms $\widehat{E_{h}}$ and $\widetilde{H_{h}}$ are numerical fluxes. In this work, we take either central fluxes,

[TABLE]

one of the following alternating flux pair

[TABLE]

or the dissipative flux inspired by the upwind flux for the Maxwell system without Kerr, linear Lorentz and Raman effects,

[TABLE]

In the rest of the paper, we will call (15) as the upwind flux. It is known that the choice of numerical fluxes is important for the properties of the schemes, such as in numerical stability, accuracy, and even computational efficiency (see Sections 5 and 6 for more discussions). We emphasize that (11c)-(11e) hold in strong sense. In the theorem below, we establish stability of the semi-discrete DG scheme which is consistent with the energy stability (8)-(9) of the PDE-ODE system (6)-(7).

Theorem 1** (Semi-discrete stability).**

Under the assumption of periodic boundary conditions, the semi-discrete DG scheme (11)-(12) with central and alternating fluxes, (13) and (14), satisfies

[TABLE]

and the DG scheme with the upwind flux (15) satisfies

[TABLE]

*where *

[TABLE]

is the discrete energy. Moreover, $\mathcal{E}_{h}\geq 0$ when $\theta\in[0,\frac{3}{4}].$

Proof.

Let $\phi=H_{h}$ in (11a), $\phi=E_{h}$ in (11b) and sum up the two equalities over all elements, we obtain

[TABLE]

Note that with both central and alternating fluxes, (13) and (14), we have

[TABLE]

and with the upwind flux (15), we have

[TABLE]

while $\sum_{j=1}^{N}\int_{I_{j}}{\partial_{x}(}E_{h}H_{h})dx=-\sum_{j=1}^{N}[E_{h}H_{h}]_{j-1/2}$, therefore (17) becomes

[TABLE]

with

[TABLE]

which is non-positive. Differentiating (12) with respect to time, and substituting it into the equation (20), we obtain

[TABLE]

which, with (11c), is equivalent to

[TABLE]

By (11d) and (11c),

[TABLE]

This gives the relation

[TABLE]

Similarly, we take $\phi=\sigma_{h}$ in (11f), sum up over all elements, use (11e), and obtain

[TABLE]

which yields

[TABLE]

On the other hand,

[TABLE]

Combining the results in (23), (25), (27), (28), we now have

[TABLE]

This becomes

[TABLE]

with the discrete energy $\mathcal{E}_{h}$ defined in (16), which is guaranteed to be nonnegative as long as $\theta\in[0,\frac{3}{4}]$. Since all model parameters are positive, clearly we have $\frac{d}{dt}\mathcal{E}_{h}\leq 0.$ ∎

We have demonstrated in the theorem above that the DG scheme with appropriate flux choices can successfully maintain the energy stability of the original system on the semi-discrete level.

4 Semi-discrete Scheme: Error Estimates

In this section, we will establish the error estimates of the semi-discrete scheme, formulated in Section 3, up to a given time $T<\infty$. The following projections, $\pi_{h}$ (defined from $L^{2}(\Omega)$ onto $V_{h}^{k}$) and $\pi_{h}^{\pm}$ (defined from $H^{1}(\Omega)$ onto $V_{h}^{k}$), will be used in the analysis.

$L^{2}$ projection $\pi_{h}$: $\pi_{h}w\in V_{h}^{k}$, such that $\forall j$

[TABLE] 2. 2.

Gauss-Radau projection $\pi_{h}^{-}$: $\pi_{h}^{-}w\in V_{h}^{k}$, such that $\forall j$

[TABLE]

and $(\pi_{h}^{-}w)^{-}_{j+\frac{1}{2}}=w_{j+\frac{1}{2}}^{-}$. 3. 3.

Gauss-Radau projection $\pi_{h}^{+}$: $\pi_{h}^{+}w\in V_{h}^{k}$, such that $\forall j$

[TABLE]

and $(\pi_{h}^{+}w)^{+}_{j-\frac{1}{2}}=w_{j-\frac{1}{2}}^{+}.$

These projections are commonly used in analyzing DG methods, and the following approximation property and estimate can be easily established [9]:

[TABLE]

with $\Pi_{h}=\pi_{h},\pi_{h}^{-}$ or $\pi_{h}^{+}$, and

[TABLE]

Here $w-\Pi_{h}w$ represents the projection error. In (32)-(35), $\|\cdot\|$, $\|\cdot\|_{\infty}$, and $\|\cdot\|_{H^{k+1}}$ stand for the $L^{2}$-norm, $L^{\infty}$-norm, and $H^{k+1}$-norm on $\Omega$, respectively. And $\|w\|_{W^{1,\infty}}=(\|w\|_{\infty}^{2}+\|\frac{dw}{dx}\|_{\infty}^{2})^{1/2}$ is the $W^{1,\infty}$-norm for $W^{1,\infty}(\Omega)$. The constant $C_{\star}$ depends on $k$ but not on $h$ or $w$. Throughout the paper, $C_{\star}$ denotes a generic constant which may depend on $k$ and mesh parameter $\delta$. If we want to emphasize the sole dependence of $k$, this generic constant will be denoted by $C_{k}$, which usually is computable. $C$ is another generic constant, which is independent of $h$, but may depend on $k$, mesh parameter $\delta$, and some Sobolev norms of the exact solution of (7) up to time $T$. There is one more generic constant $C_{\text{model}}$, which depends on some or all model parameters. Each type of generic constants may take different values at different occurrences. In the analysis, the following inverse inequality will also be needed,

[TABLE]

A direct consequence of (35) is $\|w-\Pi_{h}w\|_{\infty}\leq C_{k}\|w\|_{\infty},\forall w\in W^{1,\infty}(\Omega)$ when $\Pi_{h}=\pi_{h}^{\pm}$ (or $\forall w\in L^{\infty}(\Omega)$ when $\Pi_{h}=\pi_{h}$).

We start with decomposing the error in $E$ into two parts, $E-E_{h}=\eta_{E}-\zeta_{E}$, with $\eta_{E}=E-\pi_{h}^{E}E$, $\zeta_{E}=E_{h}-\pi_{h}^{E}E$, where $\pi_{h}^{E}$ is a projection operator onto $V_{h}^{k}$. We later also use $E-\eta_{E}=\pi_{h}^{E}E$. Similarly, one can define the decomposition of errors in other quantities, namely $u-u_{h}=u-\pi_{h}^{u}-(u_{h}-\pi_{h}^{u})=\eta_{u}-\zeta_{u}$, with $u$ being $H,P,Q,J,\sigma$. In the analysis, $\pi_{h}^{u}$ is taken to be $\pi_{h}$, the $L^{2}$-projection, for $u=E,H,P,Q,J,\sigma$, except for the following two cases: when the numerical fluxes are alternating, we take

[TABLE]

while with the upwind flux, we use

[TABLE]

See [6] (such as Lemma 2.4) for the properties of such vector-form projection operators. For the a priori error estimate in next theorem, we assume the following regularity for the exact solutions,

[TABLE]

and

[TABLE]

where the former is standard for error analysis of linear models, and the latter are needed to treat nonlinearity.

Theorem 2** (Error estimates of semi-discrete scheme).**

Assuming the periodic boundary condition and the exact solutions being as regular as (40)-(41), under the conditions that $\theta\in[0,\frac{1}{4})$ and the strength of nonlinearity is sufficiently small, the following error estimates hold for the semi-discrete DG scheme (11)-(12) with flux choices (13), (14), or (15)

[TABLE]

where

[TABLE]

Proof.

With the numerical fluxes being consistent, the proposed semi-discrete scheme is consistent. That is, (11) holds if the numerical solutions are replaced by the exact ones, while the test functions are still taken from $V_{h}^{k}$. From this, one can get the error equations,

[TABLE]

coupled with

[TABLE]

Now we take $\phi=\zeta_{H}$ in (44a), $\phi=\zeta_{E}$ in (44b), $\phi=\zeta_{\sigma}$ in (44e). We then differentiate (45) in time $t$, and take $\phi=\zeta_{E}$. Following similar steps as in the proof of Theorem 1, we get

[TABLE]

Here the non-positive term $M(\cdot,\cdot)$ is defined in (21). The four terms on the right are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Next we will take two steps to estimate the left and the right hand side of (4), respectively.

Step 1: Compared with the discrete energy in the stability analysis, the terms in the second row of the left hand side of (4) are new, and they arise from the discretizations of nonlinear terms. With arbitrarily chosen constant parameters $\rho_{\text{err}}\in(0,1)$ and $\kappa_{\text{err}}\in(0,1)$, we have

[TABLE]

under the conditions:

[TABLE]

with a computable constant $C_{k}$ from (35). Indeed, under Condition 1,

[TABLE]

holds, while Condition 2 is to ensure

[TABLE]

Step 2: Next we will estimate $\Lambda_{j}$, $j=1,\cdots,4$. Cauchy-Schwartz inequality, Young’s inequality, as well as approximation result and estimate in (32)-(35) will be used repeatedly. For $\Lambda_{1}$,

[TABLE]

We here have used ${\partial_{t}\eta}_{H}={\partial_{t}H}-\pi_{h}^{H}{\partial_{t}H}$. As for $\Lambda_{2}$, with the choice of the projection operators $\pi_{h}^{E}$ and $\pi_{h}^{H}$, one has

[TABLE]

for alternating and upwind flux; while for central flux, we have

[TABLE]

For $\Lambda_{3}$,

[TABLE]

Using the approximation property and estimate in (32)-(35), as well as the boundedness of $E$, we have

[TABLE]

hence get

[TABLE]

Term $\Lambda_{4}$ is relatively subtle, and we will proceed as follows.

[TABLE]

The constant parameters $\alpha_{j}$, $j=1,\cdots,7$, are chosen so that

[TABLE]

[TABLE]

We then further restrict the strength of the nonlinearity such that

[TABLE]

and this, with the estimate (35), can be ensured under the condition

[TABLE]

Using (60)-(53) and applying (32)-(35), we are able to bound $\Lambda_{4}$

[TABLE]

Now we can combine (4), (4), (54), (55), (56), (58), (64), and reach

[TABLE]

where $r$ is specified in (43). Finally, we apply Gronwall inequality and the estimation of projection errors in (32), and conclude that

[TABLE]

under the Conditions 1-3. Note that $\rho_{\text{err}}$ is arbitrary in $(0,1)$, Condition 1 essentially implies $\theta\in[0,\frac{1}{4})$, while Conditions 2-3 require the smallness of the strength of the nonlinearity. In our final result (42), we no longer carry the two parameters $\kappa_{\text{err}}$ and $\rho_{\text{err}}$. ∎

5 Fully Discrete Scheme and Energy Analysis

In this section, we focus on fully discrete schemes for the nonlinear PDE-ODE system (6). A particular focus will be on designing temporal discretizations, with which the fully discrete methods have provable energy stability. This turns out to be a nontrivial task for the nonlinear model examined in this work. Common choices, such as the second order leap-frog or implicit trapezoidal method, may not yield provable stability results as for the linear models. The main difficulties arise from the nonlinear Kerr and Raman terms. What we will develop in this section can be understood as novel modifications of leap-frog or implicit trapezoidal method, in the presence of these nonlinear effects. The proposed temporal discretizations are still of formal second order accuracy. We will establish the energy stability for the resulting fully discrete methods. The time discretizations developed here can be used not only in conjunction with DG spatial discretizations, but also with other type discretizations, and this will be addressed in our future work.

We design two second-order time schemes, both implicit in the ODE parts. The first scheme uses the leap-frog staggered in time for the PDE parts. Given $u^{n}_{h}(\cdot)\in V_{h}^{k}$ at $t^{n}$, with $u=H,D,E,P,Q,J,\sigma$, we look for $u^{n+1}_{h}(\cdot)\in V_{h}^{k}$ at $t^{n+1}=t^{n}+\Delta t$, with $u=H,D,E,P,Q,J,\sigma$, satisfying $\forall j$

[TABLE]

The flux terms in the scheme have no ambiguity for the central and alternating fluxes (13)-(14) with $\widehat{\widehat{E^{n}_{h}}}=\widehat{E_{h}^{n}}$, and their expressions are omitted for brevity. For the upwind flux (15), the flux terms should be defined as

[TABLE]

as in the standard leap-frog formulations. Notice that the scheme is implicit for the upwind flux, but for the alternating and central fluxes, the implicit part is only on the ODEs which can be locally solved in each element. In practice, we use a Newton’s method to obtain $E_{h}^{n+1},Q_{h}^{n+1},\sigma_{h}^{n+1},P_{h}^{n+1},J_{h}^{n+1}$ from (67b)-(67h). The main novelty of the formulation is the introduction of $Y_{h}^{n}$ in (67d) as an auxiliary variable to approximate $Y=E^{3}.$ This is motivated by the fact that $dY=3E^{2}dE$ and is defined to achieve energy stability of the fully discrete scheme as shown in Theorem 3. One does not need to store $Y_{h}^{n}$, instead only its temporal difference is needed to be substituted into (67b). Another change in the scheme for stability consideration is the discretization of $E^{2}$ term in (67f) as $E_{h}^{n}E_{h}^{n+1}$. This is motivated by theoretical analysis as shown in the proof of Theorem 3.

Similarly, our second formulation, which is a fully implicit scheme writes

[TABLE]

The scheme is of second order accuracy in time. The flux terms are defined according to their semi-discrete counterparts. For example, with the upwind flux (15), the flux terms are

[TABLE]

Theorem 3** (Fully discrete stability).**

Assuming the periodic boundary condition, then the fully discrete scheme (67) with central and alternating fluxes, (13) and (14), satisfies

[TABLE]

where

[TABLE]

is the discrete energy. In addition, $\mathcal{E}_{h}\geq 0$ if $\theta\in[0,\frac{3}{4}]$ and the CFL condition $\frac{\Delta t}{h}\leq C_{\star}\sqrt{\mu_{0}\epsilon_{0}\epsilon_{\infty}}$ is satisfied.

The fully discrete scheme (67) with the upwind flux (68) satisfies

[TABLE]

where

[TABLE]

is the discrete energy. In addition, $\mathcal{E}_{h}\geq 0$ if $\theta\in[0,\frac{3}{4}]$ and the CFL condition $\frac{\Delta t}{h}\leq\frac{C_{\star}\mu_{0}\min(1,\sqrt{\frac{\epsilon_{0}\epsilon_{\infty}}{\mu_{0}}})}{(\sqrt{2}+\min(1,\sqrt{\frac{\epsilon_{0}\epsilon_{\infty}}{\mu_{0}}}))}$ is satisfied.

Similarly, the fully discrete scheme (69) with central and alternating fluxes, (13) and (14), satisfies

[TABLE]

and that with the upwind flux (70) satisfies

[TABLE]

where

[TABLE]

It is non-negative when $\theta\in[0,\frac{3}{4}]$. In other words, the scheme (69) is unconditionally stable for all three flux choices.

Proof.

We will only prove the results for scheme (67), while the proof for the scheme (69) shares great similarity and is omitted.

Apply two time steps of (67a) and (67i), we have

[TABLE]

where

[TABLE]

Let $\phi=E_{h}^{n+1}+E_{h}^{n}$ in (67b), $\phi=H_{h}^{n+1/2}$ in (78) and sum up the two equalities over all elements, we obtain

[TABLE]

by the identity (18) in the proof of Theorem 1.

Using (67c), we obtain

[TABLE]

and here (67d) is used. By (67g) and (67h),

[TABLE]

On the other hand,

[TABLE]

By (67f) and (67e), we have

[TABLE]

Substituting (5)-(5) into (80), we have shown the energy stability (71) for the scheme (67) with the alternating and central fluxes, where the discrete energy $\mathcal{E}_{h}^{n}$ is defined in (3). When the flux is upwind, we instead get the energy stability in (73), with the discrete energy $\mathcal{E}_{h}^{n}$ defined in (74).

The final step is to find the conditions to guarantee the discrete energy to be non-negative. We define two operators

[TABLE]

From (67a), we have

[TABLE]

and similarly by (67i)

[TABLE]

Therefore,

[TABLE]

In the last equality, (91) has been used.

With the alternating or central fluxes, the last term on the right of (92) vanishes, hence

[TABLE]

On the other hand, using inverse inequality (36), one gets

[TABLE]

In addition, we take $\phi=H_{h}^{n}-H_{h}^{n-1/2}$ in (91), and reach

[TABLE]

and therefore

[TABLE]

This estimate guarantees $\mathcal{E}_{h}^{n}$ being non-negative if $\theta\in[0,\frac{3}{4}]$ and $\frac{\Delta t^{2}C_{\star}^{2}}{8h^{2}\mu_{0}}\leq\frac{\epsilon_{0}\epsilon_{\infty}}{2}$, i.e. the CFL condition $\frac{\Delta t}{h}\leq\frac{2}{C_{\star}}\sqrt{\mu_{0}\epsilon_{0}\epsilon_{\infty}}$ is satisfied. This condition can be also written as $\frac{\Delta t}{h}\leq C_{\star}\sqrt{\mu_{0}\epsilon_{0}\epsilon_{\infty}}$. Here the generic constant $C_{\star}$ depends on $k$ and mesh regularity parameter $\delta$ and is independent of $h$.

With the upwind flux in (68), using the definitions of $\widehat{\widehat{\mathcal{H}}}$ and $\widehat{\mathcal{H}}$, the last term on the right of (92) becomes

[TABLE]

Now taking into account the jump terms in the discrete energy (74), and with (92)(96), we have

[TABLE]

An analogue of (95) for the upwind flux can be obtained

[TABLE]

We can choose $\Delta t$ small enough so that $A=\frac{\Delta tC_{\star}}{2\mu_{0}h}$ satisfies $\frac{A}{1-A}\leq\frac{1}{\sqrt{2}}\min(1,\sqrt{\frac{\epsilon_{0}\epsilon_{\infty}}{\mu_{0}}})$, then

[TABLE]

Hence

[TABLE]

and

[TABLE]

Combined with (97) and (74), we have shown that $\mathcal{E}_{h}^{n}$ in (74) is nonnegative, provided $\theta\in[0,\frac{3}{4}]$ and the CFL condition $\frac{\Delta t}{h}\leq\frac{2\mu_{0}\min(1,\sqrt{\frac{\epsilon_{0}\epsilon_{\infty}}{\mu_{0}}})}{C_{\star}(\sqrt{2}+\min(1,\sqrt{\frac{\epsilon_{0}\epsilon_{\infty}}{\mu_{0}}}))}$ is satisfied. This condition can also be written as $\frac{\Delta t}{h}\leq\frac{C_{\star}\mu_{0}\min(1,\sqrt{\frac{\epsilon_{0}\epsilon_{\infty}}{\mu_{0}}})}{(\sqrt{2}+\min(1,\sqrt{\frac{\epsilon_{0}\epsilon_{\infty}}{\mu_{0}}}))}$. ∎

From the proof, one can see that the fully discrete scheme with the leap-frog temporal discretization is conditionally stable, under a CFL condition that is the same as the one for Maxwell’s equations without Kerr, linear Lorentz and Raman effects; while the fully implicit scheme is unconditionally stable.

6 Numerical Results

In this section, we demonstrate the behavior of the fully discrete schemes through two numerical examples. The simulations are performed on the rescaled equations with the scaling chosen as follows: let the reference time scale be $t_{0}$, and reference space scale be $x_{0}$ with $x_{0}=ct_{0}$ and $c=1/\sqrt{\mu_{0}\epsilon_{0}}$. Henceforth, the rescaled fields and constants are defined based on a reference electric field $E_{0}$ as follows,

[TABLE]

where for simplicity, we have used the same notation to denote the scaled and original variables. In summary, we arrive at the dimensionless Maxwell’s equations

[TABLE]

Correspondingly, the energy $\mathcal{E}(t)$

[TABLE]

satisfies the relation

[TABLE]

For the rescaled system (98), all fluxes retain their original definition except for the upwind flux, which is modified to

[TABLE]

We further refer to one alternating flux

[TABLE]

as alternating flux I, and

[TABLE]

as alternating flux II. For numerical simulations in this section, we use a uniform mesh with size $h_{j}=h=(x_{R}-x_{L})/N$ for all $j$. When solving the nonlinear system, we employ a Jacobian-free Newton-Krylov solver [31] with absolute error threshold $\epsilon=10^{-10}$.

6.1 Kink shape solutions

The first numerical test we consider is originally discussed in [39], where a traveling wave solution was constructed for the instantaneous intensity-dependent Kerr response neglecting the influence of damping, i.e., $\theta=0$, $\tau=\infty$ in (98). This yields a simplified system

[TABLE]

We use this example as an accuracy test. As shown in [39], we can find a traveling wave solution $E(x,t)=E(\xi)$, where $\xi=x-vt$, and similarly for other variables $H$, $D$, $P$ and $J$. Here, $E(\xi)$ is comprised of a kink and antikink wave, and is solved based on the following ODE

[TABLE]

The parameters are

[TABLE]

Here, $\omega_{0}$ and $\Phi(0)$ are carefully chosen such that $E$ and $\Phi$ are both $6$-periodic. The approximate solution of (101), as shown in Figure 1(a), is obtained with 160000 grid points by a third order Runge-Kutta method. This serves as the initial condition for the electric field in the Maxwell’s system (100). Furthermore, with the help of (100) and the property that all variables are traveling waves, we can obtain the initial conditions for other variables:

[TABLE]

Numerical results are provided at $t=6/v$, when the wave moves back to the same position as the initial condition. Time steps are chosen as

[TABLE]

to guarantee $(k+1)$-th order accuracy in time, and the CFL numbers are listed in Table 1. Since leap-frog scheme uses staggered time, we set CFL$=0.2/v$ for $k=1$, such that the last time step has full length $\Delta t$. This can help us avoid the influence on accuracy caused by time step changes. When $k=2$ or $3$, we do not need to do this because the time steps are already pretty small. Note that the time steps of the fully implicit scheme are taken to be much larger.

We list the errors and orders of accuracy of $E$ in Tables 2-4. All calculations give the optimal $(k+1)$-th order, except that for the central flux, if we use the leap-frog scheme, the order of accuracy will be suboptimal when $k=1$.

Next, we investigate the numerical energy behaviors with $N=400$ grid points. The results are listed in Figure 2. Since $\theta=0$ and $\tau=\infty$, following the proof in Theorem 3, we obtain that the discrete energy $\mathcal{E}_{h}^{n}$ with alternating and central fluxes satisfies

[TABLE]

where

[TABLE]

for leap-frog scheme, and

[TABLE]

for fully implicit scheme. Therefore, the schemes are energy-conserving. Figure 2 shows that the numerical results are consistent with our analysis: the leap-frog scheme conserve discrete energy up to the machine error, while the fully implicit scheme has larger errors, which is caused by larger time steps and the error from the Newton solver (we set the tolerance as $\epsilon=10^{-10}$).

On the other hand, the upwind flux is dissipative. When employing the upwind flux and leap-frog scheme, we have

[TABLE]

where the discrete energy is

[TABLE]

For the fully implicit scheme with upwind flux, the discrete energy (103) will satisfy

[TABLE]

We observe the predicted behavior numerically (Figure 2). Note that for $k=3$, the initial energy increase is caused by error from the Newton solver.

6.2 Soliton propagation

In this example, we will consider the soliton propagation in the full Maxwell model (98), similar to the setup in [19]. The computational domain is $x\in[0,45].$ The coefficients in this example are chosen as

[TABLE]

Initially, all fields are zero. The left boundary is injected with an incoming solitary wave, for which the electric field is prescribed as

[TABLE]

where $f(t)=M\ \mathrm{sech}(t-20)$. $M$ is a constant to be specified later. Similar to [19], the boundary condition of $H$ can be approximated from the linearized dispersion relation. Assuming a space-time harmonic variation $e^{i(\omega t-kx)}$ of all fields, the exact dispersion relation associated with the linear parts of the system (98) is

[TABLE]

The solution corresponding to the wave propagating to the right is

[TABLE]

Then we take the approximate value of $H$ as

[TABLE]

where $c.c.$ denotes the complex conjugate of the first term, $f^{(m)}(t)$ is the $m$-th derivative of $f(t)$, and $(\frac{1}{Z})^{(m)}$ is the $m$-th derivative of $Z=-\omega/k$ with respect to $\omega$.

We treat the right boundary as an absorbing wall corresponding to the linearized system, similar to the procedure performed in [25]. Neglecting the nonlinear effects and the delayed response in (98), we have

[TABLE]

Because only waves that propagate to the right are allowed, the left going characteristic variable $H+\sqrt{\epsilon_{\infty}}E$ is set to be zero at the right boundary $x_{R}=x_{N+1/2}$. Therefore, for semi-discrete scheme, we require

[TABLE]

This corresponds to rewriting the central flux as

[TABLE]

and rewriting the upwind flux as

[TABLE]

To guarantee better stability results for the outflowing edge, when using alternating fluxes, we employ the central flux (108) at the right boundary instead. With this boundary condition, the energy relation such as those in Theorem 1 should be adjusted accordingly. For example, we can verify that the semi-discrete scheme with alternating and central fluxes satisfy

[TABLE]

with energy

[TABLE]

and the contribution from the inflow boundary

[TABLE]

The scheme with the upwind flux satisfies

[TABLE]

with the same energy definition as in (111) and

[TABLE]

Therefore,

[TABLE]

implying energy stability.

For the fully discrete schemes, there is no ambiguity defining the fluxes (108), (109) for implicit scheme. While for the leap-frog formulations with the upwind flux (109), we use

[TABLE]

While for the other fluxes (108),

[TABLE]

Implementation-wise, with (118), at the rightmost cell $I_{N}=[x_{N-1/2},x_{N+1/2}]$, we need to solve the nonlinear system to obtain $H_{h}$ by Newton’s method. The energy relation for the resulting fully discrete scheme is summarized in the appendix.

We take $N=6400$, and the time step $\Delta t=CFL\times h.$ CFL numbers, listed in Table 5, are chosen to ensure the convergence of Newton’s method to the correct solution.

We simulate the transient fundamental ($M=1$) and second-order ($M=2$) temporal soliton evolution using various schemes with different orders. The plots of the electric field at $t=40,80$ are provided in Figures 3-6. As shown in [25, 19], a daughter pulse travels ahead the soliton-like pulse, resulting from the third-harmonic generation. This daughter pulse is much smaller in amplitude than the soliton pulse, and the frequency is about 3 times as that of the soliton pulse. The daughter pulse is evident in all simulations except with the upwind flux and $k=1$, where the numerical dissipation damps its magnitude significantly. Some reflections from the right boundary is present for the central flux. This is also observed in [25] for the finite difference scheme due to the approximate boundary conditions. As a consequence, there will be spurious oscillation near the right boundary, especially for higher order scheme. On the other hand, such oscillations are not observed for alternating fluxes or the upwind flux.

In Figures 7 - 10, we plot the transient evolution of the total energy and pulse area. Here, the pulse area is obtained by the composite trapezoidal rule between two extrema points of $E$. To distinguish the soliton pulse and the daughter pulse effectively, we only consider the soliton pulse area when $|E|\geq 0.01$. Numerical results represent high agreement between the leap-frog scheme and the fully implicit scheme. Notice that we employ the approximation boundary condition $H(x=0,t)$, and the two alternating fluxes need different inflow information, which means one uses $E(x=0,t)$ and the other one uses $H(x=0,t)$. Hence, there is a slight discrepancy between the total energy with those two fluxes, as well as the pulse area. When using the central flux, both $E(x=0,t)$ and $H(x=0,t)$ are required, therefore the energy and pulse area calculated by the central flux stay in between the two alternating fluxes, which is consistent with our analysis in (115). In Figures 7 and 8, it is observed that the total energy decreases after the entire wave entering the domain, demonstrating the energy stability of the schemes. In particular, the energy calculated from upwind flux displays slightly more damping especially when $t$ is large and $k=1$.

7 Conclusions

In this paper, we propose fully discrete energy stable schemes for 1D Maxwell’s equations in nonlinear optics. The schemes use novel treatments in temporal discretizations and discontinous Galerkin schemes in space with various choices of fluxes. We prove semi-discrete and fully discrete energy stability of the proposed methods, and provide error estimates of the semi-discrete schemes with conditions on the strength of the nonlinearity of the system. Numerical results validate the theoretical predictions, which show that the fully implicit scheme allow larger CFL numbers than the leap-frog schemes. The upwind flux exhibits more dissipation, which can damp the spurious oscillations from the boundary treatment, but also in the mean time affect the effective capturing of the daughter pulse for the soliton propagation example for low order polynomial spaces. From our experience, the alternating fluxes outperform the central and upwind fluxes in the numerical examples studied in terms of accuracy and resolution of the wave profiles. Future work includes extensions to higher dimensions, to finite difference schemes, and Fourier analysis of the semi-discrete and fully discrete DG methods for linearized Maxwell systems in dispersive media.

Appendix A Energy relation for the fully discrete schemes with non-periodic boundary conditions in Section 6.2

Here, we list the energy relation for the fully discrete schemes with boundary conditions as discussed in Section 6.2.

The results with fully implicit time discretizations are very similar to the semi-discrete case, i.e. we have that the fully implicit scheme with alternating and central fluxes satisfies

[TABLE]

and that with the upwind flux satisfies

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

Moreover, for the upwind flux, we have

[TABLE]

Thus,

[TABLE]

On the other hand, the leap-frog scheme with alternating and central fluxes satisfies

[TABLE]

where

[TABLE]

[TABLE]

Unlike the previous cases, (122) cannot be shown as non-negative, which means some energy may be injected at the right boundary in this case.

The leap-frog scheme with the upwind flux satisfies

[TABLE]

where $\Theta_{in}$, $\Theta_{out}$ and the discrete energy $\mathcal{E}_{h}^{n}$ are

[TABLE]

Note that at fully discrete level, we can only prove energy stability for fully implicit scheme with upwind flux.