Numerical simulation of a coupled system of Maxwell equations and a gas dynamic model

TL;DR

This paper presents a numerical simulation method for plasmon dynamics in metallic nanostructures, combining Maxwell equations with a gas dynamic model including quantum effects, revealing insights into harmonic generation and resonance phenomena.

Contribution

It introduces a coupled Maxwell-gas dynamic simulation framework with quantum pressure effects for plasmon behavior in nanostructures, advancing computational modeling capabilities.

Findings

Quantum pressure term causes bulk plasmon resonance.

Bound electrons influence high order harmonic generation.

The numerical method accurately simulates nonlocal effects.

Abstract

It is known that both linear and nonlinear optical phenomena can be produced when the plasmon in metallic nanostructures are excited by the external electromagnetic waves. In this work, a coupled system of Maxwell equations and a gas dynamic model including a quantum pressure term is employed to simulate the plasmon dynamics of free electron fluid in different metallic nanostructures using a discontinuous Galerkin method in two dimensions. Numerical benchmarks demonstrate that the proposed numerical method can simulate both the high order harmonic generation and the nonlocal effect from metallic nanostructures. Based on the switch-on-and-off investigation, we can conclude that the quantum pressure term in gas dynamics is responsible for the bulk plasmon resonance. In addition, for the dielectric-filled nano-cavity, a coupled effective polarization model is further adopted to investigate the optical behavior of bound electrons. Concerning the numerical setting in this work, a strengthened influence of bound electrons on the generation of high order harmonic waves has been observed.

Tables

Table 1. Table 1: Unit system

Physical quantity	Reference scale	Redefined quantity
$\tilde{L}$	$L_0=1.0\times {10}^{-9}m$	$L=\tilde{L}/L_0$
$\tilde{t}$	$t_0=L_0/c,c=1/\sqrt{{\mu }_0{ϵ}_0}$	$t=\tilde{t}/t_0$
$\tilde{\nabla }$		$\nabla =L_0\tilde{\nabla }$
$\partial /\partial \tilde{t}$		$\partial /\partial t=t_0\partial /\partial \tilde{t}$
$\widetilde{𝐄}$	$E_0=1.0\times {10}^{7}V/m$	$𝐄=\widetilde{𝐄}/E_0$
$\widetilde{𝐇}$	$H_0=E_0/Z,Z=\sqrt{{\mu }_0/{ϵ}_0}$	$𝐇=\widetilde{𝐇}/H_0$
$\widetilde{𝐉}$	$J_0=E_0/\left(Z{L}_0\right)$	$𝐉=\widetilde{𝐉}/J_0$
${\widetilde{𝐉}}_b$	$J_0$	${𝐉}_b={\widetilde{𝐉}}_b/J_0$
$\widetilde{𝐏}$	$t_0{J}_0$	$𝐏=\widetilde{𝐏}/\left(t_0{J}_0\right)$
$\tilde{\rho }$	${\rho }_0={ϵ}_0/L_0$	$\rho =\tilde{\rho }/{\rho }_0$
$\tilde{\omega }$	${\omega }_0=c/L_0$	$\omega =\tilde{\omega }/{\omega }_0$

Table 2. Table 2: L 2 superscript 𝐿 2 L^{2} error and convergence order for ρ 𝜌 \rho , H z subscript 𝐻 𝑧 H_{z} , E x subscript 𝐸 𝑥 E_{x} and E y subscript 𝐸 𝑦 E_{y} approximated by P 1 superscript 𝑃 1 P^{1} and P 2 superscript 𝑃 2 P^{2} element at T = 0.2 𝑇 0.2 T=0.2 .

Mesh		$\rho$			$H_z$			$E_x$			$E_y$
		$L^2$ error	Order		$L^2$ error	Order		$L^2$ error	Order		$L^2$ error	Order
$P^1$
$40\times 40$		0.106E-02	-		0.121E-02	-		0.754E-03	-		0.713E-03	-
$80\times 80$		0.262E-03	2.02		0.299E-03	2.02		0.202E-03	1.90		0.190E-03	1.91
$160\times 160$		0.654E-04	2.00		0.747E-04	2.00		0.525E-04	1.94		0.491E-04	1.95
$320\times 320$		0.163E-04	2.00		0.187E-04	2.00		0.134E-04	2.00		0.125E-04	1.98
$640\times 640$		0.408E-05	2.00		0.467E-05	2.00		0.337E-05	1.99		0.315E-05	1.99
$P^2$
$40\times 40$		0.299E-04	-		0.302E-04	-		0.300E-04	-		0.297E-04	-
$80\times 80$		0.374E-05	3.00		0.379E-05	2.99		0.383E-05	2.97		0.378E-05	2.97
$160\times 160$		0.467E-06	3.00		0.474E-06	3.00		0.484E-06	2.99		0.477E-06	2.99
$320\times 320$		0.584E-07	3.00		0.593E-07	3.00		0.608E-07	2.99		0.599E-07	2.99
$640\times 640$		0.730E-08	3.00		0.742E-08	3.00		0.762E-08	3.00		0.751E-08	3.00

Equations

$${\mu }_0\frac{\partial \tilde{\mathbf{H}}}{\partial \tilde{t}}+\tilde{\nabla }\times \tilde{\mathbf{E}}$$

$${ϵ}_0\frac{\partial \tilde{\mathbf{E}}}{\partial \tilde{t}}-\tilde{\nabla }\times \tilde{\mathbf{H}}$$

$$\underset{|\tilde{\mathbf{x}}|\to \infty }{\lim }({\tilde{\mathbf{H}}}^s\times \tilde{\mathbf{x}}-|\tilde{\mathbf{x}}|{\tilde{\mathbf{E}}}^s)=0,(\text{or}\underset{|\tilde{\mathbf{x}}|\to \infty }{\lim }({\tilde{\mathbf{E}}}^s\times \tilde{\mathbf{x}}+|\tilde{\mathbf{x}}|{\tilde{\mathbf{H}}}^s)=0),$$

$$\nabla \times u$$

$$\frac{\partial ({\tilde{n}}_e{\tilde{q}}_e)}{\partial \tilde{t}}+\tilde{\nabla }\cdot ({\tilde{n}}_e{\tilde{q}}_e{\tilde{\mathbf{u}}}_e)=0,$$

$${\tilde{m}}_e\left[\frac{\partial {\tilde{\mathbf{u}}}_e}{\partial \tilde{t}}+({\tilde{\mathbf{u}}}_e\cdot \tilde{\nabla }){\tilde{\mathbf{u}}}_e\right]={\tilde{q}}_e(\tilde{\mathbf{E}}+{\mu }_0{\tilde{\mathbf{u}}}_e\times \tilde{\mathbf{H}})-\tilde{\gamma }{\tilde{m}}_e{\tilde{\mathbf{u}}}_e-\frac{\tilde{\nabla }\tilde{p}}{{\tilde{n}}_e},$$

$$\frac{\partial \tilde{\mathbf{P}}}{\partial \tilde{t}}={\tilde{\mathbf{J}}}_{\mathbf{b}},$$

$$\frac{\partial {\tilde{\mathbf{J}}}_{\mathbf{b}}}{\partial \tilde{t}}+{\tilde{\gamma }}_b{\tilde{\mathbf{J}}}_{\mathbf{b}}+{\tilde{\omega }}_b^2\tilde{\mathbf{P}}$$

$$\tilde{\rho }={\tilde{n}}_e{\tilde{q}}_e,\tilde{\mathbf{J}}=\tilde{\rho }{\tilde{\mathbf{u}}}_e,$$

$$\frac{\partial U_1}{\partial t}+\nabla \cdot {\mathbf{F}}_1(U_1)=S_1(U_2,U_3),\text{in}{\Omega }_1\times [0,T],$$

$$\frac{\partial U_2}{\partial t}+\nabla \cdot {\mathbf{F}}_2(U_2)=S_2(U_1,U_2),\text{in}{\Omega }_2\times [0,T],$$

$$\frac{d{U}_3}{dt}=S_3(U_1,U_3),\text{in}{\Omega }_3\times [0,T].$$

$$U_{1}=\left(\begin{array}[]{l}H_{x}\\ H_{y}\\ H_{z}\\ E_{x}\\ E_{y}\\ E_{z}\end{array}\right),F_{1}(U_{1})=\left(\begin{array}[]{l}~{}~{}~{}0{}{}\\ -E_{z}{}\\ ~{}~{}E_{y}\\ ~{}~{}~{}0{}\\ ~{}~{}H_{z}\\ -H_{y}{}\end{array}\right),G_{1}(U_{1})=\left(\begin{array}[]{l}~{}~{}E_{z}{}{}\\ ~{}~{}~{}0{}{}{}\\ -E_{x}{}{}{}\\ -H_{z}\\ ~{}~{}~{}0{}{}{}\\ ~{}~{}H_{x}\end{array}\right),S_{1}=-\left(\begin{array}[]{l}0\\ 0\\ 0\\ \rho u_{x}+J_{bx}\\ \rho u_{y}+J_{by}\\ \rho u_{z}+J_{bz}\\ \end{array}\right)$$

$$U_{2}=\left(\begin{array}[]{l}\rho\\ \rho u_{x}\\ \rho u_{y}\\ \rho u_{z}\end{array}\right),F_{2}(U_{2})=\left(\begin{array}[]{l}\rho u_{x}\\ \rho u_{x}u_{x}+k\rho^{5/3}\\ \rho u_{x}u_{y}\\ \rho u_{x}u_{z}\end{array}\right),G_{2}(U_{2})=\left(\begin{array}[]{l}\rho u_{y}\\ \rho u_{x}u_{y}\\ \rho u_{y}u_{y}+k\rho^{5/3}\\ \rho u_{x}u_{z}\end{array}\right),$$

$$k=\frac{1}{5}{\left(\frac{\hbar}{m_{e}}\right)}^{2}{\left(\frac{3\pi^{2}}{q_{e}}\right)}^{2/3}~{},S_{2}=\left(\begin{array}[]{l}0\\ \frac{\rho q_{e}}{m_{e}}(E_{x}+u_{y}H_{z}-u_{z}H_{y})-\gamma\rho u_{x}\\ \frac{\rho q_{e}}{m_{e}}(E_{y}+u_{z}H_{x}-u_{x}H_{z})-\gamma\rho u_{y}\\ \frac{\rho q_{e}}{m_{e}}(E_{z}+u_{x}H_{y}-u_{y}H_{x})-\gamma\rho u_{z}\\ \end{array}\right)~{},$$

$$U_{3}=\left(\begin{array}[]{l}P_{x}\\ P_{y}\\ P_{z}\\ J_{bx}\\ J_{by}\\ J_{bz}\end{array}\right)~{},S_{3}=-\left(\begin{array}[]{l}J_{bx}\\ J_{by}\\ J_{bz}\\ \frac{n_{b}q_{e}^{2}}{m_{e}}E_{x}+\frac{q_{e}}{m_{e}}(J_{by}H_{z}-J_{bz}H_{y})-\gamma_{b}J_{bx}-\omega_{b}^{2}P_{x}\\ \frac{n_{b}q_{e}^{2}}{m_{e}}E_{y}+\frac{q_{e}}{m_{e}}(J_{bz}H_{x}-J_{bx}H_{z})-\gamma_{b}J_{by}-\omega_{b}^{2}P_{y}\\ \frac{n_{b}q_{e}^{2}}{m_{e}}E_{z}+\frac{q_{e}}{m_{e}}(J_{bx}H_{y}-J_{by}H_{x})-\gamma_{b}J_{bz}-\omega_{b}^{2}P_{z}\end{array}\right)~{},$$

$$V_{h,{\Omega }_1}^p:$$

$$W_{h,{\Omega }_2}^q:$$

$$V_{h,{\Omega }_3}^p:$$

$${\int }_K{U}_{1h}^{n+1}\cdot {\Phi }_{h}d\mathbf{x}$$

$$\mathcal{H}_{1}(U_{1h}^{n,int},U_{1h}^{n,ext})=\left(\begin{array}[]{l}-\bm{n}_{K}\times\frac{(Z\bm{H}+\bm{n}_{K}\times\bm{E})_{h}^{n,int}+(Z\bm{H}-\bm{n}_{K}\times\mathbf{E})_{h}^{n,ext}}{Z_{h}^{n,int}+Z_{h}^{n,ext}}\\ \\ \bm{n}_{K}\times\frac{(Y\bm{E}-\bm{n}_{K}\times\bm{H})_{h}^{n,int}+(Y\bm{E}+\bm{n}_{K}\times\mathbf{H})^{n,ext}_{h}}{Y_{h}^{n,int}+Y_{h}^{n,ext}}\end{array}\right)~{}.$$

$${\int }_K{U}_{2h}^{n+1}\cdot {\Psi }_{h}d\mathbf{x}$$

$${\mathcal{H}}_2(U_{2h}^{n,int},U_{2h}^{n,ext})=\frac{1}{2}\left[{\mathbf{F}}_2(U_{2h}^{n,int})\cdot {\mathbf{n}}_K+{\mathbf{F}}_2(U_{2h}^{n,ext})\cdot {\mathbf{n}}_K-{\alpha }^n\left(U_{2h}^{n,ext}-U_{2h}^{n,int}\right)\right],{\alpha }^n=\underset{K}{\max }|{\mathbf{F}}_2^{\prime }(U_{2h}^n)\cdot {\mathbf{n}}_K|.$$

$${\int }_K{U}_{3h}^{n+1}\cdot {\xi }_{h}d\mathbf{x}$$

$$\mathbf{E}(\mathbf{x},0)=\mathbf{H}(\mathbf{x},0)=0,\text{in}{\Omega }_1,$$

$$\rho (\mathbf{x},0)=q_e{n}_0(\mathbf{x}),\mathbf{u}(\mathbf{x},0)=0,\text{in}{\Omega }_2,$$

$$\mathbf{P}(\mathbf{x},0)={\mathbf{J}}_b(\mathbf{x},0)=0,\text{in}{\Omega }_3.$$

$$\frac{\partial U_1}{\partial t}+\nabla \cdot {\mathbf{F}}_1(U_1)=S_p(U_1,U_7),\text{in}{\Omega }_p\times [0,T],$$

$$\frac{d{U}_7}{dt}=S_7,\text{in}{\Omega }_p\times [0,T],$$

$$S_{p}=\left(\begin{array}[]{l}Q_{x}+(\sigma_{x}-\sigma_{y})H_{x}\\ Q_{y}+(\sigma_{y}-\sigma_{x})H_{y}\\ Q_{z}-(\sigma_{x}+\sigma_{y})H_{z}\\ P_{x}+(\sigma_{x}-\sigma_{y})E_{x}\\ P_{y}+(\sigma_{y}-\sigma_{x})E_{y}\\ P_{z}-(\sigma_{x}+\sigma_{y})E_{z}\\ \end{array}\right),~{}U_{7}=\left(\begin{array}[]{l}Q_{x}\\ Q_{y}\\ Q_{z}\\ P_{x}\\ P_{y}\\ P_{z}\\ \end{array}\right),~{}\\ S_{7}=\left(\begin{array}[]{l}-\sigma_{x}Q_{x}-\sigma_{x}(\sigma_{x}-\sigma_{y})H_{x}\\ -\sigma_{y}Q_{y}-\sigma_{y}(\sigma_{y}-\sigma_{x})H_{y}\\ -\sigma_{x}\sigma_{y}H_{z}\\ -\sigma_{x}P_{x}-\sigma_{x}(\sigma_{x}-\sigma_{y})E_{x}\\ -\sigma_{y}P_{y}-\sigma_{y}(\sigma_{y}-\sigma_{x})E_{y}\\ -\sigma_{x}\sigma_{y}E_{z}\\ \end{array}\right).$$

$${\sigma }_i={\sigma }_m{\left(\frac{d_i}{\delta }\right)}^n,i=x,y,$$

$$R(0)=e^{-2{\sigma }_m\delta /(n+1)},$$

$$E_z=E_{A}sin({\omega }_{m}t)e^{(-\frac{4\pi (t-t_d{)}^2}{t_b^2})},{\omega }_m=\frac{2\pi }{{\lambda }_0},$$

Full text

††journal: arXiv

1 Introduction

In recent years, there are increasing interests in complex optical phenomena associated with metallic nanostructures. One of them is the classical local optical response that features the macroscopic properties of materials. However, for small metallic nanoparticles and metallic clusters, the experimental studies [5] have retrieved a size-dependent surface resonance shift and a multiple bulk resonance at frequencies above the plasma frequency in the extinction cross sections (ECS). These novel phenomena, which are not observed in the classical local response, are due to the instantaneous response to the excitation in a nonlocal manner: $\mathbb{D}(\bm{x},\omega)=\epsilon_{0}\int\bm{\epsilon}(\bm{x},\bm{x}^{{}^{\prime}},\omega)\cdot\mathbb{E}(\bm{x}^{{}^{\prime}},\omega)d\bm{x}^{{}^{\prime}}$ , and thus are termed as the so-called nonlocal effect. Here, $\mathbb{D}$ is the electric displacement, $\epsilon_{0}$ is the electric permittivity in vacuum, $\bm{\epsilon}$ is the relative permittivity tensor, $\omega$ is the frequency, and $\mathbb{E}$ denotes the electric field. These plasmonic responses have found their wide applications in biosensing [1], plasmonic waveguiding [3] and cancer therapy [15]. Another attractive optical phenomenon is the high order harmonic generation, in particular, the second harmonic generation (SHG). Physically, the SHG is an optical process in which an electromagnetic wave at the fundamental frequency interacts with metallic nanoparticles to generate a new wave with twice of the fundamental frequency. Experiments [12, 13, 19, 24] have shown that the second harmonic wave can be generated from various metamaterials. The SHG is of great importance with broad applications such as in development of the laser sources [21], the optical parametric amplifiers [28], and imaging and microscopy technology [4].

There have been many numerical methods and models proposed to simulate the above optical phenomena in nanostructures. In [16], the authors generalized a finite-difference time domain (FDTD) method to simulate the SHG from metallic nanostructures by using a fully coupled fluid-Maxwell system (called the nonlinear hydrodynamic Drude model) that was derived from the cold-plasma wave equations [29] and the Maxwell equations. In the coupled system, the charge density $\rho$ depends on the divergence of the electric field $\mathbf{E}$, namely, $\rho=\epsilon_{0}\nabla\cdot\mathbf{E}$. Since the normal electric field is discontinuous at the dielectric-metal interface, the computation of the charge density is challenging. The authors introduced a smooth transition layer between the metal and dielectric materials so that the ion density varies from its bulk value to zero smoothly, leading to an efficient computation of $\nabla\cdot\mathbf{E}$. In [14], a fully second order hydrodynamic model has been employed to explore the mechanism of nonlocal feature on the nonlinear high order harmonic generation. In [2], an energy stable discontinuous Galerkin method has been designed for the Maxwell equations in Kerr-Raman-type nonlinear optical media for the simulation of third harmonic wave. Numerical investigations on the high order harmonic generation using a perturbation hydrodynamic model [29] and a fully hydrodynamic model [7, 9] are also reported in literatures. Concerning the nonlocal effect, a mixed finite element method (FEM) adopting the Nédélec element has been developed in [10] for simulating the nonlocal effect of a groove and a nanowire by using a nonlocal hydrodynamic Drude model in the two-dimensional frequency domain. In [23], the authors presented a DGTD method to solve a linearized nonlocal dispersion model for studying the nonlocal dispersion effect from the interaction of light with nanometer scale metallic structures. In [27], the authors applied a hybridizable discontinuous Galerkin (HDG) method to solve the Maxwell equations coupled with the nonlocal hydrodynamic Drude model in the frequency domain for computing the nonlocal electromagnetic effect from a two-dimensional gold nanowire and a three-dimensional periodic annular nanogap structure.

In this paper, we employ a coupled system of the gas dynamic equations including the pressure term and the Maxwell equations (termed as the modified nonlinear hydrodynamic Drude (MNHD) model in the current work) to simulate both the high-order harmonic generation and the nonlocal effect in two-dimensional metallic nanostructures using the high order Runge-Kutta discontinuous Galerkin (RKDG) method [6]. The gas dynamic model is essentially the Euler equations governing the motion of the electron fluid in metallic nanostructures. Being slightly different from the fully coupled fluid-Maxwell system in [16] where the high order harmonic generation has been successfully observed in numerical results, the MNHD model introduces the Thomas-Fermi pressure which characterizes the electron fluid equation of state. Meanwhile, the classical linear Drude model which can be derived from the MNHD model has been applied successfully to study the nonlocal effect from small metallic nano-particles by introducing a current diffusion term [18, 26]. These two facts indicate the potential of the MNHD model for the simulation of both the nonlocal effect and the high order harmonic generation under an uniform framework, and this is indeed one of the motivations of current work. To our knowledge, there have not been literatures addressing these two important optical phenomena in an uniform numerical model and method yet. Moreover, in order to simulate the SHG from a metallic nanostructure adjacent to a nanostructure with some kind of non-metallic material, we couple the MNHD model with an effective polarization model [22] to numerically investigate the influence of bound electrons on the generation of high order harmonics. A strengthened influence of bound electrons on the generation of high order harmonic waves associated to the numerical setting in this work has been successfully observed. Finally, we point out that, due to the ability of the discontinuous Galerkin (DG) method on dealing with discontinuous physical quantities in simulations, we will compute directly the charge density based on the mass conservation equation other than Gauss’s law.

This paper is organized as follows. In Section 2, we introduce the mathematical governing equations, and then present in details the numerical scheme for solving the coupled system in Section 3. In Section 4, numerical tests on the nonlocal effect and the second harmonic generation are presented to show the efficient performance of the model and numerical method. Conclusions are finally presented in Section 5.

2 Mathematical model

2.1 Maxwell equations

The governing equations for the propagation of electromagnetic fields are the Maxwell equations

[TABLE]

Here, ${\epsilon}_{0},{\mu}_{0}$ are the permittivity and the permeability in free space, respectively, $\tilde{\mathbf{E}}=\tilde{\mathbf{E}}^{i}+\tilde{\mathbf{E}}^{s}$ and $\tilde{\mathbf{H}}=\tilde{\mathbf{H}}^{i}+\tilde{\mathbf{H}}^{s}$ are the total fields, with $\tilde{\mathbf{E}}^{i}$ and $\tilde{\mathbf{H}}^{i}$ being the incident fields, and $\tilde{\mathbf{E}}^{s}$ and $\tilde{\mathbf{H}}^{s}$ being the scattered fields satisfying the Silver-M$\ddot{u}$ller radiation condition

[TABLE]

where $\tilde{\mathbf{J}}$ denotes the current density generated by the motion of electrons in a metallic nanostructure and thus identically equals to zero outside the metallic nanostructure, and $\mathbf{\tilde{J}}_{b}$ is the bound current density which is trivial if the effect of bound electrons is neglected. In this paper, we take an assumption of $z$-invariance, and therefore the following notation for the curl operator is used for the vector field $\mathbb{u}=(u_{x},u_{y},u_{z})$:

[TABLE]

2.2 Gas dynamic equations

As being shown in Figure 1, $\Omega_{2}\subset\mathbb{R}^{2}$ is a bounded domain occupied by a metallic nanostructure, which is excited by an external electromagnetic field with the electric field $\tilde{\mathbf{E}}(\tilde{\mathbf{x}},\tilde{t})$ and the magnetic field $\tilde{\mathbf{H}}(\tilde{\mathbf{x}},\tilde{t})$. On the other hand, the infinite mass of the ions in metal is further assumed, and it implies that the ions density $\tilde{n}_{0}(\bm{x})$ is a time independent function. Therefore, the positively charged ions merely play a role in providing background for the motion of electrons without making any contribution to the current density.

As a result, in terms of continuum mechanics, the motion of free electrons in the metallic nanostructure $\Omega_{2}$ satisfies the following Euler equations

[TABLE]

[TABLE]

where $\tilde{q}_{e},\tilde{m}_{e},\tilde{n}_{e},\tilde{\mathbf{u}}_{e}$ denote the electron charge, mass, number density and velocity field, respectively. $\tilde{\gamma}=1/\tilde{\tau}$ denotes the relaxation time, the average time of collisions between the electrons and the ions. $\tilde{p}={(3\pi^{2})}^{2/3}(\hbar^{2}/5\tilde{m}_{e})\tilde{n}_{e}^{5/3}$, with $\hbar$ being the Planck’s constant, is the electron gas pressure evaluated by the Thomas-Fermi theory [20]. In this work, we term the equations (1)-(5) as the modified nonlinear hydrodynamic Drude (MNHD) model. Since the coupling MNHD system is a self-consistent model with strong nonlinearity and diffusion effects, one expects to be able to observe both the high-order harmonic generation and the nonlocal optical response in numerics using this model.

2.3 Polarization model

Let $\Omega_{3}\subset\mathbb{R}^{2}$ (see Figure 1) be a bounded domain filled with some kind of non-metallic material, being adjacent to $\Omega_{2}$. In the domain $\Omega_{3}$, the electrons are subject to the linear and nonlinear restoring forces and can not leave far away from their atomic nucleus. We call these electrons as the bound electrons. Let $\mathbf{\tilde{P}}$ be the polarization and $\mathbf{\tilde{J}}_{b}$ be the bound current density, respectively, an effective polarization model for bound charges can be obtained by using Newton’s second law and neglecting the nonlinear restoring forces [22]

[TABLE]

Here, $\tilde{n}_{b}$, $\tilde{\gamma}_{b}$, $\tilde{\omega}_{b}$ denote the constant ion density, damping coefficient and resonance frequency for the bound electrons, respectively, and ${c}$ denotes the speed of light in vacuum. We point out that this model defined on the domain $\Omega_{3}$ is coupled with the MNHD model for the numerical simulation only when we consider the numerical test on the influence of bound electrons (the last part in Section 4.2).

2.4 Compact nondimensionalized forms

Assume that $\Omega_{1}\subset\mathbb{R}^{2}$ is a finite truncation such that $(\Omega_{2}\cup\Omega_{3})\subset\Omega_{1}$. Defining the charge density $\tilde{\rho}$ and the current density $\tilde{\mathbf{J}}$ as follows

[TABLE]

we can write the equations (1)-(2) and (4)-(5) in the following compact hyperbolic system after a nondimensionalization procedure according to Table 1

[TABLE]

[TABLE]

where $\bm{F}_{1}(U_{1})=\left(F_{1}(U_{1}),G_{1}(U_{1})\right)$, $\bm{F}_{2}(U_{2})=\left(F_{2}(U_{2}),G_{2}(U_{2})\right)$. We could also formulate the polarization model into the following compact ordinary differential system according to Table 1

[TABLE]

In the above formulations,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

$T$ is the final simulation time, and the subscript $x$ ($y$ or $z$) denotes the $x$ ($y$ or $z$) component of corresponding unknowns. Finally, we indicate that the radiation condition (3) will be replaced with an approximate boundary condition on the boundary $\partial\Omega_{1}$, which is to be discussed in the next section.

3 Numerical schemes

In this section, we will present the numerical method solving the equations (8)-(10). Let $\mathscr{T}_{h}$ be a partition of $\Omega_{1}$. For each element $K\in\mathscr{T}_{h}$, we define the following finite dimensional discrete spaces consisting of piecewise polynomials with the degree at most $k$

[TABLE]

where $\tilde{\mathscr{T}}_{h}=\{K\in{\mathscr{T}}_{h}:K\subset\Omega_{2}\}$, and $\hat{\mathscr{T}}_{h}=\{K\in{\mathscr{T}}_{h}:K\subset\Omega_{3}\}$. We assume that the boundary of each subdomain $\Omega_{i}$, $i=1,2,3$, belongs to the set of boundary of $K$, or contains the vertex nodes of $K$.

3.1 Schemes with the forward Euler time discretization

We start introducing the schemes with the first order forward Euler method for the time discretization, and the higher order time discretization will be discussed in Section 3.2. The proposed schemes evolve the numerical solutions $U_{1h}$, $U_{2h}$ and $U_{3h}$, which are assumed to be available at $t=t^{n}$, denoted by $U_{1h}^{n}\in V_{h,\Omega_{1}}^{6}$, $U_{2h}^{n}\in W_{h,\Omega_{2}}^{4}$ and $U_{3h}^{n}\in V_{h,\Omega_{3}}^{6}$, and will be computed at $t=t^{n+1}=t^{n}+\Delta t^{n}$, denoted by $U_{1h}^{n+1}$, $U_{2h}^{n+1}$ and $U_{3h}^{n+1}$.

3.1.1 Updating $U_{1h}^{n+1}$

To get $U_{1h}^{n+1}$, we apply to (8) with the DG method for the space discretization and the first order forward Euler method for the time discretization. That is, to look for $U_{1h}^{n+1}\in V_{h,\Omega_{1}}^{6}$, for $\forall\,\Phi_{h}\in V_{h,\Omega_{1}}^{6}$ and $\forall K\,\in\mathscr{T}_{h}$, such that

[TABLE]

where ${\mathcal{H}_{1}(\cdot,\cdot)}$ denotes the numerical flux evaluated on the interface between two adjacent elements, and $U_{1h}^{n,int},~{}U_{1h}^{n,ext}$ are the traces of $U_{1h}^{n}$ on $\partial K$ evaluated from the interior and exterior of element $K$. In this paper, we employ the upwind numerical flux [17] given by

[TABLE]

where $Z=\frac{1}{Y}=\sqrt{\frac{\mu}{\epsilon}}$ denotes the local impedance, and $\bm{n}_{K}$ denotes the unit outward normal of $K$.

3.1.2 Updating $U_{2h}^{n+1}$

To get $U_{2h}^{n+1}$, we apply to (9) with the DG method for the space discretization and the first order forward Euler method for the time discretization. That is, to look for $U_{2h}^{n+1}\in W_{h,\Omega_{2}}^{4}$, for $\forall\,\Psi_{h}\in W_{h,\Omega_{2}}^{4}$ and $\forall K\,\in\tilde{\mathscr{T}}_{h}$, such that

[TABLE]

where $\mathcal{H}_{2}(\cdot,\cdot)$ denotes the numerical flux evaluated on the interface between two adjacent elements, and $U_{2h}^{n,int},~{}U_{2h}^{n,ext}$ are the traces of $U_{2h}^{n}$ on $\partial K$ evaluated from the interior and exterior of element $K$. In this paper, we employ the Lax-Friedrichs numerical flux given by

[TABLE]

3.1.3 Updating $U_{3h}^{n+1}$

We look for $U_{3h}^{n+1}\in V_{h,\Omega_{3}}^{6}$, for $\forall\,\xi_{h}\in V_{h,\Omega_{3}}^{6}$ and $\forall\,K\in\hat{\mathscr{T}}_{h}$, such that

[TABLE]

3.2 Schemes with high order time discretizations

In the previous subsection, we have discussed the first order time discretization. To increase the accuracy in the time domain, strong stability preserving (SSP) high-order time discretizations [8] can be used, and we employ the third order TVD Runge-Kutta method [6] for the time discretization in this work.

3.3 Numerical ingredients

Numerical investigations on complex optical phenomena are significantly correlated with numerical settings. In this subsection, we will describe in details numerical ingredients, including the initial conditions, and the interface and absorbing boundary conditions, etc..

3.3.1 Initial conditions

Before excited by the electromagnetic fields, the electrons in the nanostructures are at rest if the thermal effect is ignored. Under this circumstance, the electron number density is equal to the ion number density $n_{0}(\bm{x})$ so that nanostructures are electrically neutral, where $n_{0}(\bm{x})$ could be evaluated via the plasma frequency $\omega_{p}=\sqrt{\frac{n_{0}q_{e}^{2}}{m_{e}}}$. Therefore, the initial conditions are set as follows

[TABLE]

3.3.2 Boundary conditions on $\Gamma_{2}=\partial\Omega_{2}$

In order to solve the MNHD model, two boundary conditions need to be prescribed, one of which is the boundary condition on the metal-vacuum interface $\Gamma_{2}=\partial\Omega_{2}$. At the microscopic level, the charge density $\rho$ varies continuously across the dielectric-metal interface, and there is actually a transition region with a scale of a few atomic diameters where the charge density changes gradually down to be trivial [26, 16]. However, the thickness of the transition layer is a negligible scale compared to the finest mesh that we are able to afford in domain discretizations. This fact makes it impossible for us to implement the ab-$initio$ boundary condition in a macroscopic model [11]. We use a natural boundary condition in this work, namely $\frac{\partial\rho}{\partial\bm{n}}=0$. For the current density $\mathbf{J}$, we employ the so-called $slip$ boundary condition, i.e. $\mathbf{n}\cdot\mathbf{J}=0$. It implies that the current density is prohibited to travel out of the nanoparticle surface in a normal direction with respect to the interface while a tangential current shift is allowed.

3.3.3 Artificial boundary conditions on $\Gamma_{1}=\partial\Omega_{1}$

As it is mentioned above, we need to employ an artificial boundary surrounding $\Omega_{2}\cup\Omega_{3}$, denoted by $\Gamma_{1}=\partial\Omega_{1}$, for the practical computation. In this paper, we use the uniaxial perfectly matched layer (PML) [17] to absorb the electromagnetic waves propagating through the boundary $\Gamma_{1}$ except for those tests with particular specifications. Let $\Omega_{p}$ be the PML region (see Figure 2) surrounding the finite truncation $\Omega_{1}$, and the modified formulations for (8) in $\Omega_{p}$ are written as

[TABLE]

and

[TABLE]

where

[TABLE]

Here, the parameters of the dissipative layer for absorbing the fields propagating in the $i$-th direction $\sigma_{i}$ are given by

[TABLE]

where $d_{i},\delta,n$ denote the distance from the PML-vacuum interface, the thickness of PML, and the degree of polynomials, respectively. $\sigma_{m}$ is the maximum electric conductivity which can be determined by

[TABLE]

where $R(0)$ denotes the theoretical reflection at normal incidence. The equations (14)-(15) are solved by the RKDG method as well.

3.3.4 TF/SF technique

In order to stimulate the plasmon resonances in the metallic nanostructures, the initial conditions presented in 3.3.1 are not sufficient, and an extra appropriate wave source should be added during the computation. In this paper, we take the incident wave as a $z$-polarized Gaussian pulse modulated by the sine function

[TABLE]

where ${\lambda}_{0}$ is the carrier center wavelength, ${E}_{A}$ is the peak amplitude, ${t}_{d}$ is the pulse duration, and ${t}_{b}$ is the optical bandwidth (see Figure 3). To implement this source injection in simulations, since the popular method of hard source in computational electromagnetics may produce backward-scattered waves in a longtime simulation, we apply in this work the total-field/scattered-field (TF/SF) technique [25] which requires to divide the computational domain into a total-field zone and a scattered-field zone (see Figure 4) through a virtual TF/SF boundary $F$ inside $\Gamma_{1}$. An incident wave is then introduced into the total-field zone from the virtual boundary without introducing any nonphysical effects. Meanwhile, a simple process on numerical fluxes along the virtual boundary allows us to realize this purpose efficiently and accurately (see Figure 5).

4 Numerical results

We present in this section numerical results which are computed based on the coupling system of (8)-(9) ($\bm{J}_{b}=0$ and $\Omega_{3}$ is empty) except for the last case in Section 4.2 where we compute the MNHD model (8)-(9) together with the polarization model (10). In addition, all the simulations are performed with both the piecewise linear ($P^{1}$) and quadratic ($P^{2}$) elements on structured grids, and if not specified, only $P^{2}$ results are presented for the illustration. The time step $\Delta t$ is dynamically determined by

[TABLE]

where $a_{x}=\textmd{max}(|u_{x}|+c^{f},1.0)$, $a_{y}=\textmd{max}(|u_{y}|+c^{f},1.0)$, $c^{f}=\sqrt{\frac{5}{3}k|\rho|^{\frac{1}{3}}}$, $C_{cfl}$ is the CFL number. If not specified, the vaule of CFL number is selected to be $0.3$ for all tests, and we also point out that no limiting procedures have been used in all numerical tests reported in this work.

4.1 Accuracy test

We start with a manufactured example to study the accuracy of our schemes for MNHD model. For simplicity, we choose $k=0,\frac{q_{e}}{m_{e}}=1,\gamma=1$ in the equations (8) and (9) and employ the following as the manufactured solutions

[TABLE]

where $\alpha=\frac{\sqrt{2}}{2}$. The computational domain is $\Omega_{2}=\Omega_{1}=[0,1]\times[0,1]$. Periodic boundary conditions are applied in both $x-$ and $y-$directions. In Table 2, we present the $L^{2}$ errors and the corresponding order of accuracy for $\rho$, $H_{z}$, $E_{x}$ and $E_{y}$ at $T=0.2$. It shows that the orders of accuracy are optimal for both $P^{k}$ solutions with $k=1,2$.

4.2 High-order harmonic generation

The high-order harmonic generation is a nonlinear optical process sensitive to the configurations of metallic nanostructures. In this test, we simulate the high-order harmonic generation from two typical metallic nanostructures, namely, an array of rectangular nanostructures (see the left of Figure 6) and an array of L-shaped nanostructures (see the middle of Figure 6). In addition, we also consider an array of L-shaped nanostructures with metallic materials and an array of rectangular nanostructures with non-metallic materials (see the right of Figure 6) for purpose of studying the influence of bound electrons on the generation of high order harmonic waves. In our simulations, the nanostructures are arranged periodically in $x-$direction with plasma frequency ${\omega}_{p}=4.560\times 10^{-2}$, and ${\gamma}=2.160\times 10^{-4}$. For simplicity, we only focus on one single structure by using periodic boundary condition in $x-$direction and PML in $y-$direction [16]. Measurements on the linear response and the SHG signal can be computed in terms of the integrals of the electric fields $E_{z}$ and $E_{x}$ along the probing line $S$ as follows

[TABLE]

We employ an uniform mesh $\Delta x=\Delta y=2.0$ in the simulation and the final time $T=1.499\times 10^{5}$.

Firstly, we study the high order harmonic generation from the rectangular metallic nanostructures whose dimension is shown in Figure 6 (Left). In Figure 7, we plot the time history of linear response and SHG. As being shown in Figure 7, there is no SHG from the rectangular metallic particle. This result is reasonable since the rectangular metallic nanostructure possesses a perfect centrosymmetric property which leads to a vanishing second-order nonlinear optical susceptibility tensor $\chi^{(2)}$ prohibiting the SHG. However, the rectangular metallic nanostructure does allow for the third harmonic generation (THG) since the symmetric property would not remove the third-order nonlinear optical susceptibility tensor $\chi^{(3)}$. In Figure 8 (Left), spectrums of these responses are presented. In order to show the spectrums apparently, a zoom-in plot of the spectrum is presented in Figure 8 (Right). It can be observed that, the third-order harmonic generation ($400$), the fifth-order harmonic generation ($240$) and the seventh-order harmonic generation ($171.4$) are captured in the simulation. These figures also indicate that, apart from the nonlinear optical responses, the MNHD model maintains the description of the linear optical interaction.

Then, we take an investigation on the high order harmonic generation from the L-shaped metallic nanostructures (see Figure 6 (middle) for the setup). In Figure 9, we plot the time history of the linear response and the SHG. In Figure 10 (Left), spectrums of these responses are presented. Since the centrosymmetry is broken in this case, the second order harmonic generation appears. As we can see from Figure 10 that the fundamental wave ($1200$) can be frequency-doubled after propagating through the nanostructure, and the higher order harmonic waves, such as the third-order harmonic generation ($400$), the fifth-order harmonic generation ($240$) and the seventh-order harmonic generation ($171.4$) in the transmitted $E_{z}$, are also captured in the simulation. The above results associated with both the rectangular and the L-shaped metallic nanostructures are in a good agreement with those presented in [16].

Finally, we investigate the high order harmonic generation from an L-shaped metallic nanostructure, with a rectangular non-metallic material occupied nanostructure being located at its corner (see the right of Figure 6), to study the influence of bound electrons on the generation of high order harmonic waves. The simulation setting is the same as the L-shaped case. In the simulations, we make a selection as ${n}_{b}=58.0$, ${\omega}_{b}=\sqrt{{n}_{b}{q}_{e}^{2}/{m}_{e}}$, ${\gamma}_{b}=0.8{\gamma}$. In Figure 11 (Left), spectrums of the linear and nonlinear optical responses are presented, in comparison with spectrums obtained from the L-shaped metallic nanostructure in the previous test. It can be observed that the high-order harmonics are enhanced when the interaction between the bound electrons and the external electromagnetic waves has been taken into account.

4.3 Nonlocal effect

For small optical particles, especially for particles with size down to subwavelength, apart from the surface plasmon, the bulk plasmon could be excited as well. In this test, we solve the MNHD model to investigate the nonlocal effect for the Ag nanowire. As in [10], we choose the nanowire with radius ${r}=2.0$, the plasma frequency ${\omega}_{p}=2.885\times 10^{-2}$, and the damping constant ${\gamma}=0.01{\omega}_{p}$. PML boundary conditions are employed in both $x-$ and $y-$directions. To resolve the ECS in the time domain, we introduce the $x-$polarized incident wave propagating in $y-$direction by the TF/SF technique and collect the Fourier-transformed total field and scattered field on the TF/SF boundary. The ECS is calculated as follows

[TABLE]

where $\mathbf{n}$ denotes the outward unit normal to the TF/SF boundary $F$, and $\mathbf{S}_{\vartriangle}$ denotes the time averaged Poynting vector

[TABLE]

We try to recover the normalized ECS from $0.4{\omega}_{p}$-$1.4{\omega}_{p}$ [10]. In the simulations, we adopt a short Gaussian pulse with ${t}_{b}=2.998,~{}{\lambda}_{0}=300$ carrying effective information between the range $0.4{\omega}_{p}$-$1.4{\omega}_{p}$ to perform a broad-band calculation. We simulate the optical interaction for a long time such that the scattered field decays adequately into a steady state and we set the final time ${T}=5.996\times 10^{4}$ for this purpose. In Figure 12, we present the ECS being normalized with respect to the diameter of nanowire as a function of the normalized angular frequency $\omega/\omega_{p}$. It can be seen from Figure 12 that the bulk plasmon resonances beyond the plasma frequency can be excited with the MNHD model. A slight blue-shift (from $\omega/\omega_{p}=0.703$ to $\omega/\omega_{p}=0.714$) of the surface plasmon resonance has been retrieved in the ECS as well. In order to make a further exploration on the source for the appearance of bulk plasmon resonances, we switch off the different nonlinear terms (quantum pressure term: $k\rho^{5/3}$, convection terms: $\rho u_{i}u_{j},i,j=x,y,z$, magnetic terms: $u_{i}H_{j}-u_{j}-H_{i},(i,j)=(x,y),(y,z),(z,x)$) in the equation (9), and make the calculation of the ECS, respectively. It can be observed from Figure 12 that, apart from the ECS calculated from the MNHD model without the quantum pressure term, all other ECSs are in perfect match with the calculated ECSs associated to the full MNHD model. This indicates that, among all three terms, the quantum pressure makes a unique contribution to appearance of the nonlocal effect. In Figure 13, we provide with the Fourier-transformed current density at $\omega/\omega_{p}=1.1835$. As it is shown, the resonances, namely the bulk plasmon, get excited in the bulk of nanowire.

5 Conclusions

In this paper, we develop a RKDG method for Maxwell equations nonlinearly coupled with gas dynamic models with both the quantum pressure and bound electrons being taken into account. Proper initial and boundary conditions, coupled with the DG method, have been designed for the efficient numerical simulation. Numerical results show that the high order harmonic waves can be produced from the L-shaped nanostructure and the rectangular nanostructure, and the bulk plasmon resonance can be excited for the metallic nanowire. The effect of bound electrons on the generation of high order harmonic waves have been confirmed in numerics. Meanwhile, a switch-off comparison confirms that the quantum pressure term in the MNHD model is essentially important for the production of nonlocal effects. Theoretical study on the quantum pressure terms, numerical investigations on the influence of spill-out electrons on these important and complex optical phenomena, and three dimensional simulations will be envisioned in our future work.

Acknowledgments

The work of M. Lyu is partially supported by the Chongqing University Graduate Student Research Innovation Project (Project No. CYS15018). The work of L. Jiang is partially supported by the Research Grants Council of Hong Kong GRF (Grant No. 17209918). The work of M. Li is partially supported by NSFC (Grant Nos. 11871139, 11701055). The work of L. Xu is partially supported by a Key Project of the Major Research Plan of NSFC (Grant No. 91630205) and a grant of NSFC (Grant No. 11771068).