Abstract: In this work, we develop energy stable numerical methods to simulate electromagnetic waves propagating in optical media where the media responses include the linear Lorentz dispersion, the instantaneous nonlinear cubic Kerr response, and the nonlinear delayed Raman molecular vibrational response. Unlike the first-order PDE-ODE governing equations considered previously in Bokil et al. (J Comput Phys 350: 420-452, 2017) and Lyu et al. (J Sci Comput 89: 1-42, 2021), a model of mixed-order form is adopted here that consists of the first-order PDE part for Maxwell’s equations coupled with the second-order ODE part (i.e., the auxiliary differential equations) modeling the linear and nonlinear dispersion in the material. The main contribution is a new numerical strategy to treat the Kerr and Raman nonlinearities to achieve provable energy stability property within a second-order temporal discretization. A nodal discontinuous Galerkin (DG) method is further applied in space for efficiently handling nonlinear terms at the algebraic level, while preserving the energy stability and achieving high-order accuracy. Indeed with d_E as the number of the components of the electric field, only a d_E×d_E nonlinear algebraic system needs to be solved at each interpolation node, and more importantly, all these small nonlinear systems are completely decoupled over one time step, rendering very high parallel efficiency. We evaluate the proposed schemes by comparing them with the methods in Bokil et al. (2017) and Lyu et al. (2021) (implemented in nodal form) regarding the accuracy, computational efficiency, and energy stability, by a parallel scalability study, and also through the simulations of the soliton-like wave propagation in one dimension, as well as the spatial-soliton propagation and two-beam interactions modeled by the two-dimensional transverse electric (TE) mode of the equations.

Key words: Maxwell’s equations, Kerr and Raman, Discontinuous Galerkin method, Energy stability