---
author:
- 'Vrushali A. Bokil [^1]'
- 'Yingda Cheng [^2]'
- 'Yan Jiang [^3]'
- 'Fengyan Li [^4]'
bibliography:
- 'bibfile.bib'
- 'bib\_Li.bib'
- 'Bokil.bib'
- 'cheng.bib'
- 'cheng\_papers.bib'
- 'bokil\_papers.bib'
- 'mimetic.bib'
title: 'Energy Stable Discontinuous Galerkin Methods for Maxwell’s Equations in Nonlinear Optical Media'
---

Maxwell’s equations, nonlinear dispersion, discontinuous Galerkin method, energy stability, error estimates.

Energy relation for the fully discrete schemes with non-periodic boundary conditions in Section \[sec:num2\]
============================================================================================================

Here, we list the energy relation for the fully discrete schemes with boundary conditions as discussed in Section \[sec:num2\].

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 $$\mathcal{E}_h^{n+1}- \mathcal{E}_h^{n}=-\frac{\Delta t}{4 \omega_p^2 \tau} \int_\Omega (J_h^{n+1} +J_h^{n})^2 dx-\frac{a\theta \Delta t}{8\omega_v^2 \tau_v}  \int_\Omega  (\sigma_h^{n+1} +\sigma_h^{n})^2  dx -\Delta t\Theta_{in}^{n} -\Delta t \Theta^{n}_{out}\le -\Delta t\Theta^{n}_{in},$$ and that with the upwind flux satisfies $$\begin{aligned}
 \mathcal{E}_h^{n+1}- \mathcal{E}_h^{n} &=&
-\frac{\Delta t}{4 \omega_p^2 \tau} \int_\Omega (J_h^{n+1} +J_h^{n})^2 dx-\frac{a\theta \Delta t}{8\omega_v^2 \tau_v}  \int_\Omega  (\sigma_h^{n+1}+\sigma_h^{n})^2  dx 
-\frac{\Delta t}{8\sqrt{\epsilon_\infty}} \sum_{j=1}^{N-1}[H_h^{n}+H_h^{n+1}]_{j+1/2}^2\\
&&-\frac{\Delta t\sqrt{ \epsilon_\infty}}{8}\sum_{j=1}^{N-1}[E_h^n+E_h^{n+1}]_{j+1/2}^2 -\Delta t\Theta_{in}^{n} -\Delta t \Theta^{n}_{out}\le -\Delta t\Theta^{n}_{in},\notag\end{aligned}$$ where $$\begin{aligned}
\mathcal{E}_h^n&=&\int_{\Omega}  \frac{1}{2} (H_h^{n})^2 + \frac{ \epsilon_\infty}{2} (E_h^n)^2 +  \frac{1}{2\omega_p^2} (J_h^n)^2 + \frac{\omega_0^2}{2 \omega_p^2}   (P_h^n)^2+ \frac{a\theta}{4\omega_v^2} (\sigma_h^n)^2 + \frac{a\theta}{2}  Q_h^n (E_h^n)^2  \notag\\
&& + \frac{3 a (1-\theta)}{4} (E_h^n)^4+\frac{a\theta}{4}(Q_h^n)^2  dx.\notag\end{aligned}$$ $$\begin{aligned}
0 \le\Theta^{n}_{out}&=&\left\{\begin{array}{ll}
\frac{1}{16\sqrt{\epsilon_{\infty}}} ((H^{n+1}_{h}+H^{n}_{h})^{-}_{N+1/2}
-\sqrt{\epsilon_{\infty}}(E^{n+1}_{h}+E^{n}_{h})^{-}_{N+1/2})^2, & \mbox{ for central and alternating fluxes},\\
\frac{1}{8\sqrt{\epsilon_{\infty}}}((H^{n+1}_{h}+H^{n}_{h})^{-}_{N+1/2})^2
+\frac{\sqrt{\epsilon_{\infty}}}{8}((E^{n+1}_{h}+E^{n}_{h})^{-}_{N+1/2})^2, & \mbox{ for upwind flux},\\
\end{array}
\right.\end{aligned}$$

$$\Theta^{n}_{in}=\left\{\begin{array}{ll}
    \frac{1}{8} \left(E(0,t^{n+1})+E(0,t^{n})\right) \left(H^{n+1}_{h}+H^{n}_{h}\right)^{+}_{1/2}
    +\frac{1}{8} \left(H(0,t^{n+1})+H(0,t^{n})\right)\\ \left(E^{n+1}_{h}+E^{n}_{h}\right)^{+}_{1/2}, 
    & \mbox{ for central flux},\\
    \frac{1}{4} \left(H(0,t^{n+1})+H(0,t^{n})\right) \left(E^{n+1}_{h}+E^{n}_{h}\right)^{+}_{1/2}, & \mbox{ for alternating flux I},\\
    \frac{1}{4} \left(E(0,t^{n+1})+E(0,t^{n})\right) \left(H^{n+1}_{h}+H^{n}_{h}\right)^{+}_{1/2}, & \mbox{ for alternating flux II},\\
    \frac{1}{8} \left(E(0,t^{n+1})+E(0,t^{n})\right) \left(H^{n+1}_{h}+H^{n}_{h}\right)^{+}_{1/2}
    +\frac{1}{8} \left(H(0,t^{n+1})+H(0,t^{n})\right)\\ \left(E^{n+1}_{h}+E^{n}_{h}\right)^{+}_{1/2}
    +\frac{1}{8\sqrt{\epsilon_\infty}}(H^{n}_{h}+H^{n+1}_{h})^{+}_{1/2}[H^{n}_{h}+H^{n+1}_{h}]_{1/2}\\
    +\frac{\sqrt{\epsilon_\infty}}{8}(E^{n+1}_{h}+E^{n}_{h})^{+}_{1/2}[E^{n+1}_{h}+E^{n}_{h}]_{1/2},
    & \mbox{ for upwind flux}.\\
\end{array}
\right.$$

Moreover, for the upwind flux, we have $$\begin{aligned}
\Theta_{in}=& \frac{1}{16\sqrt{\epsilon_\infty}}[H^{n+1}_{h}+H^{n}_{h}]^{2}_{1/2} +\frac{\sqrt{\epsilon_\infty}}{16}[E^{n+1}_{h}+E^{n}_{h}]^{2}_{1/2}\\ 
& +\frac{1}{16\sqrt{\epsilon_\infty}}\left( (H^{n+1}_{h}+H^{n}_{h})^{+}_{1/2} +\sqrt{\epsilon_\infty} \left( E(0,t^{n+1})+E(0,t^{n}) \right) \right)^2\\ &+\frac{1}{16\sqrt{\epsilon_\infty}}\left( \left( H(t^{n+1},0)+H(t^{n},0) \right) +\sqrt{\epsilon_\infty} \left( E^{n+1}_{h}+E^{n}_{h} \right) \right)^2\\
& -\frac{1}{8\sqrt{\epsilon_\infty}} \left( H(t^{n+1},0)+H(t^{n},0) \right)^2 -\frac{\sqrt{\epsilon_\infty}}{8} \left( E(t^{n+1},0)+E(t^{n},0) \right)^2.\end{aligned}$$ Thus, $$\begin{aligned}
 \mathcal{E}_h^{n+1}- \mathcal{E}_h^{n}\leq
\frac{1}{8\sqrt{\epsilon_\infty}}\Delta t \left( H(t^{n+1},0)+H(t^{n},0) \right)^2 +\frac{\sqrt{\epsilon_\infty}}{8} \Delta t \left( E(t^{n+1},0)+E(t^{n},0) \right)^2.\end{aligned}$$ On the other hand, the leap-frog scheme with alternating and central fluxes satisfies $$\mathcal{E}_h^{n+1}- \mathcal{E}_h^{n}=-\frac{\Delta t}{4 \omega_p^2 \tau} \int_\Omega (J_h^{n+1} +J_h^{n})^2 dx-\frac{a\theta \Delta t}{8\omega_v^2 \tau_v}  \int_\Omega  (\sigma_h^{n+1} +\sigma_h^{n})^2  dx -\Delta t\Theta^{n}_{in} -\Delta t \Theta^{n}_{out} 
\le  -\Delta t\Theta^{n}_{in} -\Delta t \Theta^{n}_{out},$$ where $$\begin{aligned}
\mathcal{E}_h^n&=&\int_{\Omega}  \frac{1}{2} H_h^{n+1/2} H_h^{n-1/2} + \frac{\epsilon_\infty}{2} (E_h^n)^2 +  \frac{1}{2\omega_p^2} (J_h^n)^2 + \frac{\omega_0^2}{2 \omega_p^2}   (P_h^n)^2 + \frac{a\theta}{4\omega_v^2} (\sigma_h^n)^2 + \frac{a\theta}{2}  Q_h^n (E_h^n)^2,  \notag \\
&& + \frac{3a (1-\theta)}{4} (E_h^n)^4+\frac{a\theta}{4}(Q_h^n)^2  dx, \notag\\
\Theta^{n}_{out}&=& \frac{\sqrt{\epsilon_\infty}}{16} \left( \left(E^{n}_{h}+E^{n+1}_{h}\right)^{-}_{N+1/2} - \frac{2}{\sqrt{\epsilon_\infty}}\left(H^{n+1/2}_{h}\right)^{-}_{N+1/2}\right)^2 +\frac{1}{16\sqrt{\epsilon_\infty}}\left(H^{n+1/2}_{h}\right)^{-}_{N+1/2}\notag\\ &&\left( H^{n-1/2}_{h}-2H^{n+1/2}_{h}+H^{n+3/2}_{h}\right)^{-}_{N+1/2}, \label{eq:nosign}\end{aligned}$$ $$\begin{aligned}
\Theta^{n}_{in}&=&\left\{\begin{array}{ll}
\frac{1}{4}\left(E(0,t^{n})+E(0,t^{n+1})\right) \left(H^{n+1/2}_{h}\right)^{+}_{1/2} + \frac{1}{4}H(0,t^{n+1/2}) \left(E^{n}_{h}+E^{n+1}_{h}\right)^{+}_{1/2}, & \mbox{ for central flux},\\
\frac{1}{2}H(0,t^{n+1/2}) \left(E^{n}_{h}+E^{n+1}_{h}\right)^{+}_{1/2},& \mbox{ for alternating flux I},\\
 \frac{1}{2}\left(E(0,t^{n})+E(0,t^{n+1})\right) \left(H^{n+1/2}_{h}\right)^{+}_{1/2}, & \mbox{ for alternating flux II}.\\
\end{array}
\right.\end{aligned}$$ Unlike the previous cases, 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 $$\begin{aligned}
 \mathcal{E}_h^{n+1}- \mathcal{E}_h^{n} &=&
-\frac{\Delta t}{4 \omega_p^2 \tau} \int_\Omega (J_h^{n+1} +J_h^{n})^2 dx-\frac{a\theta \Delta t}{8\omega_v^2 \tau_v}  \int_\Omega  (\sigma_h^{n+1}+\sigma_h^{n})^2  dx \\
&& -\frac{\Delta t}{8\sqrt{\epsilon_\infty}} \sum_{j=1}^{N-1}[H_h^{n-1/2}+H_h^{n+1/2}]_{j+1/2}^2
-\frac{\Delta t\sqrt{ \epsilon_\infty}}{8} \sum_{j=1}^{N-1}[E_h^n+E_h^{n+1}]_{j+1/2}^2 \notag\\
&&-\Delta t\Theta_{out}^{n} -\Delta t\Theta_{in}^{n},\notag
$$ where $\Theta_{in}$, $\Theta_{out}$ and the discrete energy $\mathcal{E}_{h}^{n}$ are $$\begin{aligned}
\mathcal{E}_h^n&=&\int_{\Omega}  \frac{1}{2} H_h^{n+1/2} H_h^{n-1/2} + \frac{\epsilon_\infty}{2} (E_h^n)^2 +  \frac{1}{2\omega_p^2} (J_h^n)^2 + \frac{\omega_0^2}{2 \omega_p^2}   (P_h^n)^2+ \frac{a\theta}{4\omega_v^2} (\sigma_h^n)^2 \notag \\
&& + \frac{a\theta}{2}  Q_h^n (E_h^n)^2  + \frac{3 a (1-\theta)}{4} (E_h^n)^4+\frac{a\theta}{4}(Q_h^n)^2  dx\notag\\
&&+\frac{\Delta t}{8\sqrt{\epsilon_\infty}}\sum_{j=1}^{N-1} ([H_h^{n-1/2}] [H_h^{n-1/2}+H_h^{n+1/2}])_{j+1/2} \notag \\
&&+\frac{\Delta t}{8\sqrt{\epsilon_\infty}} (H_h^{n-1/2})^{-}_{N+1/2} (H_h^{n-1/2}+H_h^{n+1/2})|^{-}_{N+1/2}, \notag \\
\Theta_{out}^{n}&=& \frac{1}{8\sqrt{\epsilon_\infty}}\left( (H_h^{n-1/2}+H_h^{n+1/2})^{-}_{N+1/2} \right)^2
+\frac{\sqrt{ \epsilon_\infty}}{8} \left( (E_h^n+E_h^{n+1})^{-}_{N+1/2} \right)^2,\notag\\
\Theta_{in}^{n}&=& \frac{1}{4}\left(E(0,t^{n})+E(0,t^{n+1})\right) \left(H^{n+1/2}_{h}\right)^{+}_{1/2} + \frac{1}{4}H(0,t^{n+1/2}) \left(E^{n}_{h}+E^{n+1}_{h}\right)^{+}_{1/2}\notag\\
&&+\frac{1}{8\sqrt{\epsilon_\infty}}\left(H^{n+1/2}_{h}\right)^{+}_{1/2}[H^{n-1/2}_{h}+2H^{n+1/2}_{h}+H^{n+3/2}_{h}]_{1/2}
+\frac{\sqrt{\epsilon_\infty}}{8}(E^{n}_{h}+E^{n+1}_{h})^{+}_{1/2}[E^{n}_{h}+E^{n+1}_{h}]_{1/2}\notag\\
&=&\frac{1}{4\sqrt{\epsilon_\infty}}[H^{n+1/2}_{h}]^{2}_{1/2} +\frac{\sqrt{\epsilon_\infty}}{16}[E^{n+1}_{h}+E^{n}_{h}]^{2}_{1/2} +\frac{1}{4\sqrt{\epsilon_\infty}}\left( (H^{n+1/2}_{h})^{+}_{1/2} +\sqrt{\epsilon_\infty} \frac{E(0,t^{n+1})+E(0,t^{n})}{2} \right)^2\notag\\ &&+\frac{1}{4\sqrt{\epsilon_\infty}}\left( H(t^{n+1/2},0) +\sqrt{\epsilon_\infty} \frac{ E^{n+1}_{h}+E^{n}_{h}}{2}  \right)^2
 -\frac{1}{2\sqrt{\epsilon_\infty}} \left( H(t^{n+1/2},0) \right)^2 -\frac{\sqrt{\epsilon_\infty}}{8} \left( E(t^{n+1},0)+E(t^{n},0) \right)^2\notag\\
 &&+\frac{1}{8\sqrt{\epsilon_\infty}}\left(H^{n+1/2}_{h}\right)^{+}_{1/2}[H^{n-1/2}_{h}-2H^{n+1/2}_{h}+H^{n+3/2}_{h}]_{1/2}.\label{eq:nosign2}\end{aligned}$$ Note that at fully discrete level, we can only prove energy stability for fully implicit scheme with upwind flux.

[^1]: Department of Mathematics, Oregon State University, Corvallis, OR 97331 U.S.A. [bokilv@math.oregonstate.edu]{}.

[^2]: Department of Mathematics, Michigan State University, East Lansing, MI 48824 U.S.A. [ycheng@math.msu.edu]{}. Research is supported by NSF grant DMS-1453661.

[^3]: Department of Mathematics, Michigan State University, East Lansing, MI 48824 U.S.A. [jiangyan@math.msu.edu]{}.

[^4]: Department of Mathematical Sciences, Rensselaer Polytechnic Institute, Troy, NY 12180 U.S.A. [lif@rpi.edu]{}. Research is supported by NSF grant DMS-1318409.
