One- and multi-dimensional stochastic Maxwell equations with additive noise are considered in this paper.
It is known that such system can be written in the multi-symplectic structure, and the stochastic energy increases linearly in time.
High order discontinuous Galerkin methods are designed for the stochastic
Maxwell equations with additive noise, and we show that the proposed methods satisfy the discrete form of the
stochastic energy linear growth property and preserve the multi-symplectic structure on the discrete level.
Optimal error estimate of the semi-discrete DG method is also analyzed. The fully discrete methods are obtained
by coupling with symplectic temporal discretizations. One- and two-dimensional numerical results are provided to demonstrate the performance of the proposed methods, and optimal error estimates and linear growth of the discrete energy can be observed for all cases.

Key words: Discontinuous Galerkin methods, Stochastic Maxwell equations, Additive noise, Multi-symplectic method, Optimal error estimate.

§ INTRODUCTION

In this paper we develop and analyze high order discontinuous Galerkin (DG) methods for one- and two-dimensional
stochastic Maxwell equations with additive noise. Maxwell equations play an important role in many physical applications,
and have been widely used in electromagnetism, electronic biology, optical imaging, etc. The general formulation of
Maxwell equations is

\begin{equation}\label{Max}
\begin{cases}
\partial_t\vect{D}=\nabla\times\vect{H}-{\vect{J}}_e,~~ \nabla\cdot\vect{D}=\rho,\\
\partial_t\vect{B}=-\nabla\times\vect{E},~~\nabla\cdot  \vect{B} = 0,
\end{cases}
\end{equation}

where $\bf{H}$ represents the magnetic field, $\bf{E}$ stands for the electric
field, $\bf{D}$ and $\bf{B}$ are the electric and magnetic flux density
respectively.  $\bf{J}_e$ is the electric current density, and $\rho$ is the
electric charge density.

Stochastic Maxwell equations are the generalized version of the deterministic Maxwell equations, which are often described as a random perturbation of the electric current density or the magnet current density by noise. The noises are commonly regarded as Brownian motion, Poisson process, etc. In [26], Rytov et al. introduced fluctuations of an electromagnetic field to obtain stochastic  Maxwell equations. Ord showed in [24] that the random walk model due to Mark Kac can be modified to produce Maxwell's field equations in 1+1 dimensions. In [18],
Liaskos et al. studied the stochastic integrodifferential equations in Hilbert spaces, and they examined the well posedness for the
Cauchy problem of the integrodifferential equations describing Maxwell equations.  The random electromagnetic fields using the spectral representation
is explored in [17], and the electromagnetic fields were coupled by Maxwell equations with a random source term.
Finite element approximations of a class of nonlinear stochastic wave equations with multiplicative noise were recently investigated in [16].
The semilinear stochastic Maxwell equations with additive noise in the following form:

\begin{align}\label{MUL}
\begin{cases}
\epsilon d\vect{E} -\nabla\times
\vect{H}dt=-\vect{J}_e(t,\vect{x},\vect{E},\vect{H})dt-\vect{J}_e^r(t,\vect{x})\circ dW, \\
\mu d\vect{H}+\nabla\times\vect{E}dt=-\vect{J}_m(t,\vect{x},\vect{E},\vect{H})dt-\vect{J}_m^r(t,\vect{x})\circ
\end{cases}
\end{align}

were studied by Chen et al. in [2], where $dW$ is a space-time mixed color noise, often driven by Brownian motion, and $\vect{J}_e$ and $\vect{J}_m$ are described as electric current and magnetic current, along with continuous bounded functions. Theoretical properties of the stochastic system (<ref>) such as regularity, energy and divergence evolution law, and symplecticity have been presented in that paper. In addition, a stochastic Runge-Kutta semidiscretization scheme was proposed for (<ref>) and proven to possess first order of mean accuracy.
In [8], David et al. developed an exponential integrator for a more generalized formulation of (<ref>), when $\bf{J}_e^r$ and $\bf{J}_m^r$ depend on $\bm{E}$ and $\bm{H}$. In a recent review article, Zhang et al. [27] presented different types of stochastic Maxwell equations with additive or multiplicative noises.

Stochastic Maxwell equations can be viewed as a type of stochastic Hamiltonian PDEs. In [12], Jiang et al. considered stochastic Hamiltonian PDEs in the form

\begin{equation}\label{defSPDE}
\end{equation}

where $M$ and $K$ are anti-symmetric matrices, and $S_1$ and $S_2$ are smooth functions of $z$.
It can be shown that the system (<ref>) satisfies the following stochastic multi-symplectic conservation law:

\[d\omega+\partial_x \kappa dt=0,\qquad \omega=MU\cdot V,\qquad \kappa=KU\cdot V,\]

which $U$ and $V$ are a pair of solutions to the variational equation

\[Md(\partial z)+K(\partial z)_x dt=\nabla^2_zS_1(z)\partial zdt+\nabla^2_zS_2(z)\partial zdW_t.\]

In [1, 10], multi-symplectic finite difference methods have been studied for the stochastic Maxwell equations with additive noise of the form

\begin{align}\label{ADD}
\begin{cases}
\epsilon d\vect{E} =\nabla \times \vect{H}dt- \lambda_1\vect{1}^T dW,\\
\mu d\vect{H}=-\nabla\times\vect{E}+\lambda_2\vect{1}^TdW,
\end{cases}
\end{align}

which can be reformulated in two slightly different formulations of (<ref>) in [1, 10] with different set of auxiliary variables introduced. Numerical methods based on these reformation have been presented and studied, and it was shown that these methods preserve stochastic multi-symplecticity on the discrete level. In addition, the linear growth property of stochastic energy was also preserved by the proposed methods. In [11], the extension to stochastic Maxwell equations with multiplicative noise was investigated.

In this work, we investigate the high-order schemes for stochastic Maxwell equations, following the recent work in [14] on DG methods for stochastic conservation laws. The DG method is a class of finite element methods that uses discontinuous piecewise polynomials as the basis functions. This method has been shown to adopt many advantages from both finite element and finite volume methods, which includes hp-adaptivity flexibility, efficient parallel implementation, the ability of handling complicated boundary conditions, etc. The DG methods were first introduced in [25] by Reed and Hill to solve transport equations, and later they were extended to solve hyperbolic conservation laws by Cockburn et al. in [3, 5, 6, 7].
There have been some recent studies in extending DG method for stochastic partial differential equations. In [14] Li et al. applied the DG method to the stochastic conservation laws with multiplicative noise

\[du+f(u)_xdt=g(x,t,u)dW_t. \]

In [15], they also proposed an ultra-weak DG method for the stochastic Korteweg-De Vries equations in the form

\[du=-(u_{xxx}+f(u)_x)dt+g(x,t,u)dW_t.\]

Optimal error estimate was proven for the semilinear equations, and numerically optimal convergence rate was also observed for many nonlinear cases.

In this paper we apply high order DG methods to stochastic Maxwell equations with additive noise (<ref>) in
one and two dimensions. The stochastic energy of the exact solutions are shown to satisfy the linear growth property, and we will demonstrate that the numerical solutions of the semi-discrete DG methods satisfy the similar energy law on the discrete level. When the standard Brownian motion $W_t$ is considered, the exactly same discrete energy law can be obtained. Following the error estimate for the deterministic Maxwell equations, we provide the optimal error estimate of the semi-discrete DG methods for the one- and two-dimensional stochastic Maxwell equations on cartesian meshes.  Furthermore, multi-symplectic property of certain DG methods for the one-dimensional deterministic multi-symplectic Hamiltonian partial differential equations (HPDEs) of the form

\[Mz_t+Kz_x=\nabla_z S(z)\]

has been recently investigated in [22].
Following the idea, we will start by establishing the multi-symplectic structure of the stochastic Maxwell equations, and then prove that the DG methods with suitable numerical fluxes can preserve the multi-symplecticity. The resulting semi-discrete methods are combined with symplectic temporal discretization. Both first order symplectic Euler method and second order symplectic partitioned Runge-Kutta (PRK) method will be introduced and analyzed in this paper.

The structure of this paper is as follows.  Section <ref> provides
theoretical results for one-dimensional stochastic Maxwell equations. We discuss the
conservation of multi-symplecticity of our DG scheme, demonstrate the energy
law of numerical solutions, and present the result on the optimal error estimate of
the proposed methods.  Section <ref> provides the same theoretical
results on DG methods for two-dimensional stochastic Maxwell equations. Temporal
discretization is provided in section <ref>. We consider both symplectic
Euler method and the second order symplectic PRK method. Numerical
results in both one and two dimensions are provided in section <ref> to validate the convergence rate and the linear energy growth property. Section <ref> contains some conclusion remarks.

Throughout this paper, $L^2$ norm is denoted by $\norm{\cdot}$, and $C$ represents a generic
positive constant independent of the spatial and temporal step size $h$ and $\Delta t$,
which can take different values in different cases. In general, $W_t$ represents a
Q-Wiener process defined on a given probability space
$(\Omega,\mathcal{F},\mathbb{P})$, which can be characterized as

\begin{equation}\label{defW}
t\ge 0,
\vect{x}\in \mathbb{R}~ \text{or}~ \mathbb{R}^2,
\end{equation}

where $\{e_m\}$ is an orthonormal basis of $L^2(\mathbb{D})$, with $\mathbb{D}\subset \mathbb{R}^d,~d=1,2$. $Q$ is a symmetric, non-negative operator such that $Tr(Q)<\infty$ and $Qe_m=\gamma_m e_m$. Furthermore, $\{\mathcal{B}_m\}$ is a sequence of independent standard Brownian motions.

§ ONE-DIMENSIONAL MAXWELL EQUATION WITH ADDITIVE NOISE

In this section, the one-dimensional computational domain is denoted by $I$, which is partitioned into subintervals $I_j=[x_{j-\frac12},x_{j+\frac12}]$
where $j=1,2,\cdots,N$. We also denote $x_j=\frac{1}{2}(x_{j-\frac12}+x_{j+\frac12})$ to be the center of each cell, and $h_j=x_{j+\frac12}-x_{j-\frac12}$ to be the mesh size. Let $h=\max_j h_j$ be the maximum mesh size. We further assume that $h/h_j$
is bounded over all $j$.  Assume $P^k(I_j)$ to be the space of polynomials of degree up to $k$ on $I_j$, and the piecewise polynomial space $V_h^k$ is defined as follows:

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

Note that functions in $V_h^k$ can be discontinuous at cell interfaces. Let $v_{j+\frac12}^+$ and $v_{j+\frac12}^-$ be the right and left limit of $v$ at the interface $x_{j+\frac12}$, and we denote by $\{ v\}_{j+\frac12} = \frac12(v_{j+\frac12}^+ + v_{j+\frac12}^-)$ and $[v]_{j+\frac12} = v_{j+\frac12}^+ - v_{j+\frac12}^-$ the average and jump of $v$ at $x_{j+\frac12}$.

We start by considering the one-dimensional (1D) stochastic Maxwell equation with additive noise of the form

\begin{equation}\label{1}
\begin{cases}
d\eta=-u_xdt-\lambda_1dW_t, &\\
\end{cases}
\end{equation}

which can be viewed as a 1D version of the model (<ref>).
The DG scheme for (<ref>) is formulated as:  for $x\in I,~ (\omega,t)\in \Omega\times [0,T]$, find $\eta_h(\omega,x,t), u_h(\omega,x,t) \in V^k_h$, such that for any test functions $ \varphi, \widetilde{\varphi}\in V_h^k$, it holds that

\begin{eqnarray}
\int_{I_j}d\eta_h\varphi(x)dx=\Big(\int_{I_j}u_h\varphi_xdx
-\int_{I_j}\lambda_1\varphi dW_tdx,   \label{phi}\\ \, \notag \\
\int_{I_j}du_h\widetilde{\varphi}(x)dx=\Big(\int_{I_j}\eta_h\widetilde{\varphi}_xdx
+\int_{I_j}\lambda_2\widetilde{\varphi} dW_tdx ,
\label{phii}
\end{eqnarray}

where the generalized alternating numerical fluxes are chosen to be

$$\widecheck{u}_h=\{u_h\}+\alpha[u_h], \qquad \widecheck{\eta}_h=\{\eta_h\}-\alpha[\eta_h],$$

for some non-zero constant $\alpha\in[-1,1]$.
Below, we will explore some theoretical properties of this DG method, including the discrete energy law, optimal error estimate and the multi-symplectic structure.

§.§ Discrete energy law

The exact solutions of the three dimensional stochastic Maxwell equation (<ref>) satisfy the linear energy growth property, as studied in [10]. Below, we start by presenting similar result for the one-dimensional model (<ref>).

Let $u$ and $\eta$ be the solutions to the model (<ref>) under periodic boundary condition. Denote the energy by $\mathcal{E}(t)=\int_I u^2(x,t)+\eta^2(x,t)dx$, then for any $t$, the global stochastic energy satisfies the following energy law

\begin{eqnarray}\label{thm2.9}
\mathcal{E}(t) = \mathcal{E}(0)  +2 \int_0^t\int_{I}(\lambda_2\eta-\lambda_1 u) dW_tdx
+(\lambda_1^2+\lambda_2^2)Tr(Q) t.
\end{eqnarray}

and, after taking the expectation,

\begin{eqnarray}\label{thm2.10}
\mathbb{E}\big(\mathcal{E}(t)\big)=\mathcal{E}(0)+(\lambda_1^2+\lambda_2^2)Tr(Q)t.
\end{eqnarray}

By utilizing the Itô's lemma and the equations (<ref>), we have

\begin{align}
&d\mathcal{E}(t)=\int_I \big(2udu+ (du)^2+2\eta d\eta +(d\eta)^2\big) dx \notag\\
&~~~~~~~=-\int_I (2u\eta_x+2\eta u_x)dtdx+2\int_I (\lambda_2u-\lambda_1 \eta) dW_tdx+(\lambda_1^2+\lambda_2^2)\int_I (dW_t)^2 dx.	\label{1denergyeq1}
\end{align}

It follows from the periodic boundary condition and the definition of $W_t$ that

\[\int_I (2u\eta_x+2\eta u_x)dtdx=0, \qquad\qquad \int_I (dW_t)^2 dx=Tr(Q)dt.\]

Therefore, (<ref>) reduces to

\begin{eqnarray}\label{beforeexp}
d\mathcal{E}(t)=2\int_I (\lambda_2u-\lambda_1 \eta)dW_tdx+(\lambda_1^2+\lambda_2^2)Tr(Q)dt.
\end{eqnarray}

Integrating (<ref>) over $t$ yields (<ref>), and taking an expectation leads to (<ref>).

Next, we show that the following discrete energy law is satisfied by the numerical solutions of the DG methods.

Let $u_h(\omega,x,t)$ and $\eta_h(\omega,x,t)$ be the numerical solutions to the DG methods (<ref>) and (<ref>).

(a): For any $t\in[0,T]$, the numerical solutions satisfy the discrete energy law

\begin{align}
&\norm{u_h(\omega,x,t)}^2+\norm{\eta_h(\omega,x,t)}^2 \notag\\
+\norm{u_h(x,0)}^2+\norm{\eta_h(x,0)}^2+(\lambda_1^2+\lambda_2^2)Kt. \label{1denergy11}
\end{align}

\begin{equation}
\mathbb{E}\Big(\norm{u_h(\omega,x,t)}^2+\norm{\eta_h(\omega,x,t)}^2\Big)=\norm{u_h(x,0)}^2+\norm{\eta_h(x,0)}^2+(\lambda_1^2+\lambda_2^2)Kt, \label{1denergy11_1}
\end{equation}

\begin{equation}
K=\sum_{j=1}^N\sum_{l=0}^k \mu_j^l\sum^\infty_{m=1} \Big(\int_{I_j}\phi_j^l\sqrt{\gamma_m}e_mdx\Big)^2, \label{added2}
\end{equation}

where  $\{\phi^l_j, j=0,\cdots,k\}$ represents the set of Legendre basis over cell $I_j$, and $\mu_j^l = (\int_{I_j}(\phi_j^l)^2dx)^{-1}$.

(b) The constant $K$ is bounded by $(k+1)Tr(Q)$ with $k$ being the polynomial degree of the DG methods.
Moreover, if there exists some constant $\alpha>0$ such that the series $\sum^\infty_{m=1}\gamma_m (K_m)^\alpha<\infty$ with $K_m=\norm{(e_m)_x}_{{\infty}}$,
we can show that $K=Tr(Q)+O(h^\alpha)$, i.e., $K$ is an approximation of $Tr(Q)$ appearing in the continuous energy law (<ref>)-(<ref>).

(a): From Itô's formula, we have

\begin{equation}\label{Ito}
d(u_h)^2=2u_hdu_h+(du_h)^2,\qquad d(\eta_h)^2=2\eta_hd\eta_h+(d\eta_h)^2.
\end{equation}

By taking the test functions $\varphi=\eta_h$ and $\widetilde{\varphi}=u_h$ in the DG methods (<ref>)-(<ref>), and summing the resulting equations up, we have

\begin{eqnarray}\label{qqq1}
\int_{I_j}((du_h)u_h+(d\eta_h)\eta_h) dx&=&\int_{I_j}(\lambda_2u_h-\lambda_1\eta_h)dW_tdx+\Big(\int_{I_j}u_h(\eta_h)_xdx
+\eta_h(u_h)_xdx\Big)dt \notag \\
&&+\Big( -\big((\{u_h\}+\alpha[u_h])\eta_h^-\big)_{j+\frac12} + \big((\{u_h\}+\alpha[u_h])\eta_h^+\big)_{j-\frac12}\Big)dt
\notag \\
&&+\Big( -\big((\{\eta_h\}-\alpha[\eta_h])u_h^-\big)_{j+\frac12} + \big((\{\eta_h\}-\alpha[\eta_h])u_h^+\big)_{j-\frac12}\Big)dt\notag \\
&=&-\Theta_{j+\frac12}+\Theta_{j-\frac12} +\int_{I_j}(\lambda_2u_h-\lambda_1\eta_h)dW_tdx,
\end{eqnarray}

\[\Theta=\left(\frac12+\alpha\right)\eta^-_hu_h^++\left(\frac12-\alpha\right)u_h^-\eta_h^+.\]

Let us represent the numerical solutions $u_h$ in the cell $I_j$ as

\[u_h(\omega,x,t)=\sum^k_{l=0}u_j^l(\omega,t)\phi_j^l(x),\]

which leads to

\[du_h=\sum^k_{l=0}du_j^l\phi_j^l.\]

Taking the test function $\widetilde{\varphi}=\phi_j^m$, $m=0,\cdots,k$, in (<ref>), we obtain

\begin{equation}\label{system}
\sum^k_{l=0}\Big(\int_{I_j}\phi_j^m\phi_j^ldx\Big)du_j^l=
\mathcal{A}_j(\eta_h;\phi_j^m)dt+\int_{I_j}\lambda_2\phi_j^mdW_tdx
\end{equation}

where the operator

\[\mathcal{A}_j(f;g)=\int_{I_j}fg_xdx-\big(\widecheck{f}g^-\big)_{j+\frac12}+
\big(\widecheck{f}g^+\big)_{j-\frac12},\]

is introduced for ease of presentation.
Denote the mass matrix by $L_j$ with the $(m,l)$ entry being
$\int_{I_j}\phi_j^l\phi_j^mdx$. Note that the orthogonal Legendre basis are chosen, therefore $L_j$ is a diagonal matrix,
so is the inverse matrix $L_j^{-1}$. Let us denote $L_j^{-1}$=diag$(\mu_j^0,\mu_j^1,\cdots,\mu_j^k)$ with $\mu_j^l=(\int_{I_j}(\phi_j^l)^2dx)^{-1}$. In addition, introduce the vectors

\[\textbf{u}_j=\bigg(u_j^0,u_j^1\cdots,u_j^k\bigg)^T,
\textbf{A}_j=\bigg(\mathcal{A}_j(\eta_h;\phi_j^0),\mathcal{A}_j(\eta_h;\phi_j^1)

and the equation (<ref>) can be rewritten as a linear system

\begin{equation}\label{linearsystem}
\end{equation}

which leads to

\begin{equation}\label{du}
\end{equation}

Therefore, we have

\begin{eqnarray}
\int_{I_j}(du_h)^2dx
=\sum^k_{l=0}\bigg(\int_{I_j}\phi_j^l\phi_j^ldx\bigg)(du_j^l)^2 \notag \\
d\textbf{u}_j=\Big(\int_{I_j}\lambda_2\Phi_jdW_tdx\Big)\cdot L_j^{-1}\Big(\int_{I_j}\lambda_2\Phi_jdW_tdx\Big)\notag\\
&=&\lambda_2^2\sum_{l=0}^k \mu_j^l\Big(\int_{I_j}\phi_j^ldW_tdx\Big)^2. \label{718663}
\end{eqnarray}

The definition of $W_t$ yields

\begin{eqnarray*}
\Big(\int_{I_j}\phi_j^ldW_tdx\Big)^2=\Big(\int_{I_j}\phi_j^l\sum^\infty_{m=1}\sqrt{\gamma_m}e_m(x)d\mathcal{B}_mdx\Big)^2
=\sum^\infty_{m=1} \Big(\int_{I_j}\phi_j^l\sqrt{\gamma_m}e_mdx\Big)^2dt.
\end{eqnarray*}

Similarly, we have

\begin{equation}
\int_{I_j} (d\eta_h)^2dx=\lambda_1^2\sum_{l=0}^k \mu_j^l\Big(\int_{I_j}\phi_j^ldW_tdx\Big)^2.
\label{718663_1}
\end{equation}

The combination of (<ref>), (<ref>), (<ref>) and (<ref>) leads to

\begin{align}
\int_{I_j} (d(u_h)^2 + d(\eta_h)^2) dx =& -2\Theta_{j+\frac12}+2\Theta_{j-\frac12} \notag\\
&+2\int_{I_j}(\lambda_2u_h-\lambda_1\eta_h)dW_tdx + (\lambda_1^2+\lambda_2^2)\sum_{l=0}^k \mu_j^l\sum^\infty_{m=1} \Big(\int_{I_j}\phi_j^l\sqrt{\gamma_m}e_mdx\Big)^2dt.
\label{added1}
\end{align}

Summing over all the cells and integrating in time from $0$ to $t$, we can obtain (<ref>). Note that $\int_I(\lambda_2u_h-\lambda_1\eta_h)dW_sdx$ is an Itô integral, thus
$\mathbb{E}\Big(\int^t_0\int_I(\lambda_2u_h-\lambda_1\eta_h)dW_sdx\Big)=0$, and taking the expectation (<ref>) leads to (<ref>).

(b): For the constant $K$ defined in (<ref>), we can bound it by

\begin{align*}
K\le & \sum_{j=1}^N\sum_{l=0}^k \mu_j^l\sum^\infty_{m=1}
\int_{I_j}(\phi_j^l)^2dx\int_{I_j}{\gamma_m}e_m^2dx
\int_{I_j}{\gamma_m}e_m^2dx\\
\int_{I}{\gamma_m}e_m^2dx=(k+1)Tr(Q),
\end{align*}

where $\mu_j^l=(\int_{I_j}(\phi_j^l)^2dx)^{-1}$ is used.

Next, we show that $K$ is an approximation of $Tr(Q)$. Note that $\phi_j^l=1$, $\mu_j^1 =1/h_j$, and $\int_{I_j}\phi_j^l dx=0$ for $l\ge 1$.
Let us define $\overline{e}_m=\frac{1}{h_j} \int_{I_j} e_m dx$, and rewrite $K$ as

\begin{eqnarray}\label{esC}
K=\sum^N_{j=1}\sum^\infty_{m=1}\frac{\gamma_m}{h_j}\Big(\int_{I_j}e_mdx\Big)^2 +
\sum^N_{j=1}\sum^\infty_{m=1}\sum_{l=1}^k\mu_j^l \gamma_m \Big(\int_{I_j} \phi_j^l
(e_m-\overline{e}_m) dx\Big)^2.
\end{eqnarray}

Using the fact that $Tr(Q)=\sum^\infty_{m=1}\gamma_m$ and $\sum^N_{j=1} \int_{I_j}e_m^2dx = 1$, we have

\begin{align}\notag
+  \sum^N_{j=1}\sum^\infty_{m=1}\sum_{l=1}^k\mu_j^l \gamma_m \Big(\int_{I_j} \phi_j^l
(e_m-\overline{e}_m) dx\Big)^2 \\
& =: \RNum{1}+\RNum{2}. \label{CandTr}
\end{align}

Notice that

\begin{eqnarray*}
\int_{I_j}e_m^2 dx=\int_{I_j}(e_m-\overline{e}_m+\overline{e}_m)^2dx
=\int_{I_j}(e_m-\overline{e}_m)^2dx+\int_{I_j}\overline{e}_m^2 dx=
\int_{I_j}(e_m-\overline{e}_m)^2dx+\frac1h\Big(\int_{I_j}e_mdx\Big)^2,
\end{eqnarray*}

which leads to

\begin{eqnarray*}
\RNum{1}=-\sum^N_{j=1}\sum^\infty_{m=1}\gamma_m\int_{I_j}(e_m-\overline{e}_m)^2dx.
\end{eqnarray*}

Since $e_m$ is uniformly bounded, there exists some constant $M$ such that $|e_m-\overline{e}_m|^{2-\alpha}\le M$, and we have the following estimate

\begin{align*}
\RNum{1}&=-\sum^N_{j=1}\sum^\infty_{m=1}\gamma_m\int_{I_j}(e_m-\overline{e}_m)^2dx \\
&\le M\sum^N_{j=1}\sum_{m=1}^\infty\gamma_m\int_{I_j}|e_m-\overline{e}_m|^\alpha dx
\le M\sum^N_{j=1}\sum_{m=1}^\infty\gamma_m (K_m)^\alpha h^{1+\alpha}	  \\
&= M \Big(\sum^N_{j=1} h^{1+\alpha} \Big) \Big(\sum_{m=1}^\infty\gamma_m (K_m)^\alpha\Big) = O(h^{\alpha}),
\end{align*}

where $K_m=\norm{(e_m)_x}_{{\infty}}$ and the assumption $\sum_{m=1}^\infty\gamma_m (K_m)^\alpha<\infty$ is used.
For the other term $\RNum{2}$, we can apply Young's inequality and follow the similar analysis as above to obtain

\begin{align*}
\RNum{2}&=\sum^N_{j=1}\sum^\infty_{m=1}\sum_{l=1}^k\mu_j^l \gamma_m \Big(\int_{I_j} \phi_j^l
(e_m-\overline{e}_m) dx\Big)^2
\le \sum^N_{j=1}\sum^\infty_{m=1}\sum_{l=1}^k \mu_j^l \gamma_m \int_{I_j} (\phi_j^l )^2dx \int_{I_j}(e_m-\overline{e}_m)^2dx\\
&= \sum^N_{j=1}\sum^\infty_{m=1}\sum_{l=1}^k\gamma_m\int_{I_j}(e_m-\overline{e}_m)^2dx
= k\sum^N_{j=1}\sum^\infty_{m=1}\gamma_m\int_{I_j}(e_m-\overline{e}_m)^2dx =  O(h^{\alpha}).
\end{align*}

The combination of these results lead to the conclusion that
$K-Tr(Q)=O(h^\alpha)$, which finishes the proof.

If $W_t$ is the standard Brownian motion, we have $e_1=1$ and $e_m=0$ for $m>1$, therefore it can be easily shown that $\RNum{1}=\RNum{2}=0$ which leads to $K= Tr(Q)$. This means that the continuous energy law (<ref>) is exactly preserved by the proposed method.

§.§ Optimal error estimates

The optimal error estimate analysis of the proposed semi-discrete DG method will be provided in this subsection.

We firstly define the generalized Radau $\mathcal{P}^\alpha$ projection operators which will be used in the analysis. On any cell $I_j$ and for any function $g(x)$, its projection $\mathcal{P}^\alpha g$ into the space $V_h^k$ is given by

\[\int_{I_j}(\mathcal{P}^\alpha g-g(x))v(x)dx=0,~\forall v(x)\in P^{k-1}(I_j),~~\text{and} ~
(\{\mathcal{P}^\alpha g\}+\alpha[\mathcal{P}^\alpha g])_{j+\frac12}=g(x_{j+\frac12}).\]

The following property on the projection error is studied in [21, 23] and will be used throughout this section.

Let $\mathcal{P}^\alpha$ (with $\alpha\neq 0$) be the generalized Radau projection defined above. There exists some constant $C$ which is independent of $h$, such that

\[\norm{\mathcal{P}^{\alpha}g-g}\le Ch^{k+1}.\]

Let $u$, $\eta\in L^2(\Omega\times[0,T];H^{k+2})$ be the strong solutions to the one-dimensional stochastic Maxwell equations with additive noise (<ref>), and $u_h,~\eta_h\in{V}^k_h$ be the numerical solutions given by the DG scheme (<ref>) and (<ref>). With the initial conditions chosen as

$$\eta_h(x,0)=\mathcal{P}^{-\alpha}\eta(x,0), ~u_h(x,0)=\mathcal{P}^{\alpha} u(x,0), $$

there holds the following error estimates

\begin{eqnarray}\label{thm2}
\norm{u-u_h}^2+\norm{\eta-\eta_h}^2 \le Ch^{2k+2}.
\end{eqnarray}

Note that both the exact solution $\eta$ and the numerical solution $\eta_h$ satisfy the equation (<ref>), and both $u$ and $u_h$
satisfy (<ref>), therefore we have the following error equations

\begin{eqnarray}
\int_{I_j}d(\eta-\eta_h)\varphi(x)dx=\Big(\int_{I_j}(u-u_h)\varphi_xdx
\int_{I_j}d(u-u_h)\widetilde{\varphi}(x)dx=\Big(\int_{I_j}(\eta-\eta_h)\widetilde{\varphi}_xdx
\end{eqnarray}

\begin{eqnarray*}
\xi^\eta=\mathcal{P}^{-\alpha}\eta-\eta_h,~~\epsilon^\eta=\mathcal{P}^{-\alpha}\eta-\eta,~~
\xi^u=\mathcal{P}^\alpha u-u_h,~~\epsilon^u=\mathcal{P}^\alpha u-u,
\end{eqnarray*}

so that we can decompose the numerical error into two terms

\begin{equation}\label{ooo}
\eta-\eta_h=\xi^\eta-\epsilon^\eta,\qquad~u-u_h=\xi^u-\epsilon^u.
\end{equation}

By choosing the test functions $\varphi=\xi^\eta$, $\widetilde{\varphi}=\xi^u$ in (<ref>) and (<ref>), and summing up the resulting equations, we obtain

\begin{align}
\int_{I_j}(d\xi^\eta\xi^\eta+d\xi^u\xi^u) dx&=\int_{I_j}(d\epsilon^\eta\xi^\eta+d\epsilon^u\xi^u )dx\notag \\
+\big((\{\xi^u\}+\alpha[\xi^u])(\xi^\eta)^+\big)_{j-\frac12}\Big)dt \notag \\
+\big((\{\epsilon^u\}+\alpha[\epsilon^u])(\xi^\eta)^+\big)_{j-\frac12}\Big)dt\notag \\
+\big((\{\xi^\eta\}-\alpha[\xi^\eta])(\xi^u)^+\big)_{j-\frac12}\Big)dt \notag \\
+\big((\{\epsilon^\eta\}-\alpha[\epsilon^\eta])(\xi^u)^+\big)_{j-\frac12}\Big)dt\notag \\
&= \int_{I_j}(d\epsilon^\eta\xi^\eta+d\epsilon^u\xi^u) dx+ \Big(\widetilde{\Theta}_{j-\frac12}-\widetilde{\Theta}_{j+\frac12}\Big)dt,
\label{err3}
\end{align}

and the last equality follows from an integration by parts and the definition of $\epsilon^u$ and $\epsilon^\eta$:

\[\int_{I_j}\epsilon^u\xi^\eta_xdx=\int_{I_j}\epsilon^\eta\xi^u_xdx=(\{\epsilon^u\}+\alpha[\epsilon^u])_{j\pm\frac12}=

By Itô's lemma, we have

\begin{eqnarray}\label{ito2}
\end{eqnarray}

Note that

\[d(\mathcal{P}^\alpha u)=\mathcal{P}^\alpha(du)=\mathcal{P}^\alpha(-\eta_xdt+\lambda_2dW_t)=
\mathcal{P}^\alpha(-\eta_xdt)+\lambda_2\mathcal{P}^\alpha(dW_t).\]

For any test function $\widetilde{\varphi}$, we have

\begin{eqnarray}\label{proj}
\int_{I_j}(d\mathcal{P}^\alpha u)\widetilde{\varphi}dx
\mathcal{P}^\alpha(dW_t)dx.
\end{eqnarray}

Subtracting (<ref>) from (<ref>), we obtain

\begin{eqnarray*}
\int_{I_j}d\xi^u\widetilde{\varphi}dx=\Big(\int_{I_j}(-\eta_h\widetilde{\varphi_x}+\mathcal{P}^\alpha
\Big)dt+\lambda_2\int_{I_j}(\mathcal{P}^\alpha(dW_t)-dW_t)\widetilde{\varphi}dx.
\end{eqnarray*}

By choosing $\widetilde{\varphi}=d\xi^u$, the first term on the right-hand side becomes zero due to Itô's calculus, and we have

\[\int_{I_j}(d\xi^u)^2 dx=\lambda_2\int_{I_j}(\mathcal{P}^\alpha(dW_t)-dW_t)d\xi^udx\le \frac12\int_{I_j}(d\xi^u)^2dx+C
\int_{I_j}(\mathcal{P}^\alpha(dW_t)-dW_t)^2dx,\]

which leads to

\begin{equation} \label{added10}
\int_{I_j}(d\xi^u)^2dx\le C\int_{I_j}(\mathcal{P}^\alpha(dW_t)-dW_t)^2dx.
\end{equation}

In the same fashion, we can have

\begin{equation}\label{added11}
\int_{I_j}(d\xi^\eta)^2 dx\le C\int_{I_j}(\mathcal{P}^{-\alpha}(dW_t)-dW_t)^2dx.
\end{equation}

The combination of (<ref>), (<ref>) with (<ref>), (<ref>) leads to

\begin{eqnarray}\label{err4}
\begin{aligned}
&\frac12\int_{I_j} \big( d(\xi^\eta)^2+d(\xi^u)^2 \big)dx\\
&= \int_{I_j}(d\epsilon^\eta\xi^\eta + d\epsilon^u\xi^u)dx+\Big(\widetilde{\Theta}_{j-\frac12}-\widetilde{\Theta}_{j+\frac12}\Big)dt+\frac12\int_{I_j}\big( (d\xi^u)^2+(d\xi^\eta)^2 \big)dx\\
&\le \int_{I_j}(d\epsilon^\eta\xi^\eta + d\epsilon^u\xi^u) dx+\Big(\widetilde{\Theta}_{j-\frac12}-\widetilde{\Theta}_{j+\frac12}\Big)dt+C\int_{I_j}(\mathcal{P}^{\pm\alpha}(dW_t)-dW_t)^2dx.
\end{aligned}
\end{eqnarray}

Summing over $I_j$ and utilizing the periodic boundary conditions yields

\begin{align}\label{err5}
&\frac12\int_{I} \big( d(\xi^\eta)^2+d(\xi^u)^2 \big)dx
\le \int_{I}(d\epsilon^\eta\xi^\eta + d\epsilon^u\xi^u) dx+C\norm{\mathcal{P}^{\pm\alpha}(dW_s)-dW_s}^2.
\end{align}

By integrating (<ref>) from $0$ to $t$ and noting that $\norm{\xi^\eta(x,0)}^2=\norm{\xi^u(x,0)}^2=0$, we have

\begin{align*}
\norm{\xi^\eta(x,t)}^2+\norm{\xi^u(x,t)}^2
& \le 2\Big(\int_0^t\int_{I} (d\epsilon^\eta\xi^\eta
+ d\epsilon^u\xi^u)dxds\Big)+C\int^t_0\norm{\mathcal{P}^{\pm\alpha}(dW_s)-dW_s}^2\\
& \le \Big(\int^t_0\norm{\xi^\eta(x,s)}^2+\norm{\xi^u(x,s)}^2ds\Big)
\end{align*}

where the projection error is used in the last inequality.
Applying the Gronwall's inequality and combining with the optimal projection error yields the desired optimal error estimate (<ref>).

§.§ Multi-symplectic structure

The stochastic Maxwell equations (<ref>) have a multi-symplectic structure, and we will show that the proposed DG method (<ref>)-(<ref>) preserves the multi-symplecticity.

Following the idea in [10], we introduce new variables

\[v_t=u,~~\zeta_t=\eta,~~P=u+\frac{1}{2}\zeta_x,~~Q=\eta+\frac{1}{2}v_x,\]

and rewrite (<ref>) into the following system:

\begin{equation}\label{2}
\begin{cases}
\frac12\zeta_x=P-u, &\\
\frac12v_x=Q-\eta,&\\
d\zeta=\eta dt.&
\end{cases}
\end{equation}

Introduce the notation $z=(u,~\eta,~v,~\zeta,~P,~Q)^T$, and the system (<ref>) can be rewritten as the multi-symplectic

\begin{equation}\label{multi}
\end{equation}

\[M=\left(
\begin{array}{cccccc}
\end{array}
\right),\qquad K=\left(\begin{array}{cccccc}
\end{array} \right), \]

\[S_1(z)=Pu+Q\eta-\frac{u^2}{2}-\frac{\eta^2}{2},\qquad S_2(z)=\lambda_1\zeta-\lambda_2v.\]

Applying the exterior derivative to the system (<ref>) yields the following variation equation for the one-form $Z$:

\begin{equation}{\label{variation}}
\end{equation}

Let $U, V\in\mathbb{R}^d$ be any solutions to the variation equation (<ref>), and define

\begin{equation*}
\omega(U,V)=MU\cdot V, \qquad \kappa(U,V)=KU\cdot V,
\end{equation*}

then the system (<ref>) can be shown to satisfy the multi-symplectic conservation law given by

\begin{equation}\label{consmulti}
\end{equation}

Multi-symplectic numerical method refers to the method that satisfies a consistent discrete version of this conservation law.

Since the new system (<ref>) is equivalent to the original model (<ref>), we can rewrite the proposed DG methods (<ref>) and (<ref>) into a consistent formulation for the system (<ref>).
We start by defining the variables $v_h$ and $\zeta_h$ as follows:  for $u_h$ and $\eta_h$
defined in (<ref>) and (<ref>),
find $v_h, ~\zeta_h\in V_h^k$, such that for all test functions $\psi, \widetilde{\psi}\in V_h^k$, it holds that

\begin{eqnarray}
&\int_{I_j}dv_h\psi dx=\int_{I_j}u_h\psi dxdt, \qquad\qquad
\int_{I_j}d\zeta_h\widetilde{\psi}dx=\int_{I_j}\eta_h\widetilde{\psi}dxdt. \label{v}
\end{eqnarray}

Next, find $P_h,~Q_h\in V_h^k$, such that for all $\varphi, \widetilde{\varphi}\in V_h^k$, it holds that

\begin{eqnarray}
&&\int_{I_j}(P_h-u_h)\varphi dx=-\frac12\Big(\int_{I_j}
\zeta_h\varphi_xdx-(\widehat{\zeta}_h\varphi^-)_{j+\frac12}
+(\widehat{\zeta}_h\varphi^+)_{j-\frac12}\Big), \label{Ph} \\
+(\widehat{v}_h\widetilde{\varphi}^+)_{j-\frac12}\Big) \label{Qh},
\end{eqnarray}

where $\widehat{v}_h=\{v_h\}+2n[v_h],~
\widehat{\zeta}_h=[\zeta_h]+2m[\zeta_h]$ for some $m,~n\in\mathbb{R}$ with $n-m=\alpha$.
By combining the derivative of (<ref>) and (<ref>) with the equations (<ref>) and (<ref>), we obtain

\begin{eqnarray}
&&\int_{I_j}dP_h\phi dx=\frac12\Big(\int_{I_j}\eta_h\phi_xdx-
+(\widehat{\eta}_h\phi^+)_{j-\frac12}\Big)dt+\int_{I_j}\lambda_2\phi dW_tdx
\label{PPh},\\
\widetilde{\phi}dW_tdx\label{QQh},
\end{eqnarray}

where $\widehat{u}_h=\{u_h\}-2m[u_h],~\widehat{\eta}_h=\{\eta_h\}-2n[\eta_h]$.
Combining (<ref>)-(<ref>), we have derived the following expression:
find $P_h,~Q_h,~u_h,~\eta_h,~v_h,~\zeta_h~\in V_h^k$ such that

\begin{eqnarray}
\begin{aligned}
&\int_{I_j}(P_h-u_h)\varphi dx=-\frac12\Big(\int_{I_j}
\zeta_h\varphi_xdx-(\widehat{\zeta}_h\varphi^-)_{j+\frac12}
+(\widehat{\zeta}_h\varphi^+)_{j-\frac12}\Big), \label{DG1} \\
+(\widehat{v}_h\widetilde{\varphi}^+)_{j-\frac12}\Big), \label{DG2}\\
&\int_{I_j}dP_h\phi dx=\frac12\Big(\int_{I_j}\eta_h\phi_xdx-
+(\widehat{\eta}_h\phi^+)_{j-\frac12}\Big)dt+\int_{I_j}\lambda_2\phi dW_tdx,
\label{DG3}\\
\widetilde{\phi}dW_tdx,\label{DG4}\\
&\int_{I_j}dv_h\psi dx=\int_{I_j}u_h\psi dxdt,\label{DG5}\\
&\int_{I_j}d\zeta_h\widetilde{\psi}=\int_{I_j}\eta_h\widetilde{\psi}dxdt , \label{DG6}
\end{aligned}
\end{eqnarray}

hold for any $\varphi,~\widetilde{\varphi},~\phi,~\widetilde{\phi},~\psi,~\widetilde{\psi}~\in
The numerical fluxes take the form

\[\widehat{u}_h=\{u_h\}-2m[u_h],~~\widehat{\eta}_h=\{\eta_h\}-2n[\eta_h],

Note that the method (<ref>) can be rewritten into the corresponding DG scheme for the new system (<ref>): find $z_h\in (V_h^k)^6$, such that for all  $\bm{\varphi}\in (V^k_h)^6$, it holds that

\begin{eqnarray}\label{DG}
\int_{I_j}Mdz_h\cdot \bm{\varphi} dx-\Big(\int_{I_j}Kz_h\cdot\bm{\varphi}_xdx-
\big(\widehat{Kz_h}\cdot\bm{\varphi}^-\big)_{j+\frac12}+
\big(\widehat{Kz_h}\cdot\bm{\varphi}^+\big)_{j-\frac12}\Big)dt \notag \\
=\int_{I_j}\nabla S_1(z_h)\cdot \bm{\varphi} dxdt+\int_{I_j}\nabla S_2(z_h)\cdot
\bm{\varphi} dW_tdx,
\end{eqnarray}

where $\widehat{Kz_h}$ is the numerical flux in the form

\[\widehat{Kz_h}=K\{z_h\}+A[z_h],
\qquad
A= \left(\begin{array}{cccccc}
\end{array} \right). \]

Applying the exterior derivative to this scheme (<ref>) leads to the variational equation

\begin{eqnarray}\label{varia}
\int_{I_j}MdZ_h\cdot \varphi dx-\Big(\int_{I_j}KZ_h\cdot\varphi_xdx-
\big(\widehat{KZ_h}\cdot\varphi^-\big)_{j+\frac12}+
\big(\widehat{KZ_h}\cdot\varphi^+\big)_{j-\frac12}\Big)dt %\notag \\
=\int_{I_j}\nabla^2 S_1(z)Z_h\cdot \varphi dxdt,
\end{eqnarray}

Following the proof of multi-symplecticity of DG methods in [22], we have the following results.

For any $U_h, V_h\in {V}^k_h$, we have

\begin{eqnarray}\label{lemma1-1}
KU^-_h\cdot V_h^--\widehat{KU_h}\cdot V_h^-+
\widehat{KV_h}\cdot U_h^-=
KU^+_h\cdot V_h^+-\widehat{KU_h}\cdot V_h^++
\widehat{KV_h}\cdot U_h^+=\mathcal{F}_K(U_h,V_h),
\end{eqnarray}

\[\mathcal{F}_K(U_h,V_h)=\{KU_h\cdot V_h\}-\widehat{KU_h}\cdot \{V_h\}+
\widehat{KV_h}\cdot \{U_h\}.\]

Let  $U_h, V_h$ be any solutions to the variational equation (<ref>), we have the semi-discrete version of
the multi-symplectic conservation laws

\begin{equation}\label{dismul}
\int_{I_j} d(MU_h\cdot V_h)dx-\Big(\mathcal{F}_K(U_h,V_h)_{j+\frac12}
\end{equation}

By Itô's lemma, we have

\begin{equation}\label{added12}
d(MU_h\cdot V_h)=MdU_h\cdot V_h+MU_h\cdot dV_h+MdU_h\cdot dV_h.
\end{equation}

Note that $U_h(\omega,x,t), V_h(\omega,x,t)\in{V}^k_h$, which can be rewritten as

\[U_h(\omega,x,t)=\sum^k_{l=0}\varphi_j^l(x)U^l_j(\omega,t),\qquad V_h(\omega,x,t)=\sum^k_{l=0}\varphi_j^l(x)V^l_j(\omega,t),\]

where  $\{\varphi_j^l(x),~l=0,\cdots, k\}$ is a set of basis.
Therefore one has

\[dU_h=\sum^k_{l=0}\varphi_j^l(x)dU^l_j,\qquad dV_h=\sum^k_{l=0}\varphi_j^l(x)dV^l_j.\]

In the variational equation (<ref>), we set $dZ_h$ to be $dU_h$, and take the test function $\varphi$ to be $dV_h$. As the second and third terms in (<ref>) are both drift terms, we can conclude that $\int_{I_j}MdU_h\cdot dV_hdx=0$.
By combining (<ref>) and (<ref>), and utilizing the fact that $M$ is anti-symmetric and $\nabla^2S_1$ is symmetric, we obtain

\begin{eqnarray*}
&&\int_{I_j}d(MU_h\cdot V_h)dx=\int_{I_j}(MdU_h\cdot V_h+MU_h\cdot dV_h ) dx
=\int_{I_j} (MdU_h\cdot V_h-MdV_h\cdot U_h ) dx\\
&&=\Bigg(\int_{I_j}KU_h\cdot(V_h)_xdx-\big(\widehat{KU_h}\cdot V_h^-\big)_{j+\frac12}
V_h^+\big)_{j-\frac12}\Bigg)dt+\int_{I_j}\nabla^2S_1U_h\cdot V_hdxdt\\
&&-\Bigg(\int_{I_j}KV_h\cdot(U_h)_xdx-\big(\widehat{KV_h}\cdot U_h^-\big)_{j+\frac12}
+\big(\widehat{KV_h}\cdot U_h^+\big)_{j-\frac12}\Bigg)dt-\int_{I_j}\nabla^2S_1V_h\cdot
U_hdxdt \\
&&=\Big(\big(KU_h^-\cdot V_h^--\widehat{KU_h}\cdot V_h^-+\widehat{KV_h}\cdot U_h^- \big)_{j+\frac12}
-\big(KU_h^+\cdot V_h^+-\widehat{KU_h}\cdot V_h^++\widehat{KV_h}\cdot U_h^+\big)_{j-\frac12} \Big)dt\\
&&= \Big(\mathcal{F}_K(U_h,V_h)_{j+\frac12}
\end{eqnarray*}

where the last equality follows from Lemma <ref>. This finishes the proof.

§ TWO-DIMENSIONAL STOCHASTIC MAXWELL EQUATIONS WITH ADDITIVE NOISE

In this section, we will present the discontinuous Galerkin methods for two-dimensional stochastic Maxwell equations with additive noise on cartesian meshes, and study the stability, error estimate and multi-symplecticity of the proposed methods.

The two-dimensional rectangular computational domain is set to be $I\times J$, and we consider the rectangular partition with the cells denoted by $I_i\times J_j=[x_{i-\frac12},x_{i+\frac12}]\times [y_{j-\frac12},y_{j+\frac12}]$ for $i=1,2,\cdots, N_x$ and $j=1,2\cdots, N_y$.
Let $x_i=\frac12(x_{i-\frac12}+x_{i+\frac12})$, and $y_j=\frac12(y_{j-\frac12}+y_{j+\frac12})$ . Furthermore, we define the mesh size
in both directions as $h_{x,i}=x_{i+\frac12}-x_{i-\frac12}$, $h_{y,j}=y_{j+\frac12}-y_{j-\frac12}$, with $h_x=\max_i h_{x,i}$, $h_y=\max_j
h_{y,j}$ and $h=\max(h_x,h_y)$ being the maximum mesh size. Similar to the
one-dimensional case, we define the two dimensional piecewise polynomial space $\mathbb{V}_h^k$ as follows:

\[\mathbb{V}_h^k=\{v(x,y): v|_{I_i\times J_j}\in Q^k(I_i\times J_j)=P^k(I_i)\otimes P^k(J_j),~~i=1,2,\cdots,N_x; j=1,2,\cdots,N_y\}.\]

The two-dimensional stochastic Maxwell equations with additive noise take the form

\begin{equation}\label{2d}
\begin{cases}
\end{cases}
\end{equation}

The DG scheme for (<ref>) is formulated as follows: find $E_h,~S_h,~T_h\in
\mathbb{V}_h^k$, such that for all test functions $\varphi,~\psi,~\phi\in\mathbb{V}_h^k$, it holds that

\begin{eqnarray}
\int_{J_j}\int_{I_i}dE_h\varphi dxdy&=&-\int_{J_j}\Bigg(
\int_{I_i}T_h\varphi_xdx-\big(\widetilde{(T_h)_{\alpha_1}}\varphi^-\big)_{i+\frac12,y}+
\big(\widetilde{(T_h)_{\alpha_1}}\varphi^+\big)_{i-\frac12,y}\Bigg)
dydt \notag \\
+\big(\widetilde{(S_h)_{-\alpha_2}}\varphi^+\big)_{x,j-\frac12}\Bigg)dxdt\notag \\
\int_{J_j}\int_{I_i}\lambda_1\varphi dW_tdxdy, \label{2ddgg1}\\
\int_{J_j}\int_{I_i}dS_h\psi dxdy &=&
\int_{I_i}\Bigg(\int_{J_j}E_h\psi_ydy-\big(\widetilde{(E_h)_{\alpha_2}}\psi^-\big)_{x,j+\frac12}
+\big(\widetilde{(E_h)_{\alpha_2}}\psi^+\big)_{x,j-\frac12}\Bigg)dxdt \notag \\
\lambda_2\psi dW_tdxdy , \label{2ddgg2}\\
\int_{J_j}\int_{I_i} dT_h \phi dxdy
&=& -\int_{J_j}\Bigg(\int_{I_i}E_h\phi_x dx -\big(\widetilde{(E_h)_{-\alpha_1}}\phi^-\big)_{i+\frac12,y}
\phi^+\big)_{i-\frac12,y}\Bigg)dydt\notag \\
\phi dW_tdxdy, \label{2ddgg3}
\end{eqnarray}

where the generalized alternating numerical fluxes are defined as follows:

\[\widetilde{g_\alpha}=\{g\}+\alpha[g], \qquad \text{for}~ g\in\mathbb{V}_h^k~\text{and}~ \alpha\in\{\pm\alpha_1,\pm\alpha_2\}\subset\mathbb{R}.\]

For the ease of presentation, we also introduce the following operators: for $\alpha\in\mathbb{R}$,

\[\mathcal{A}_{I_i}(f,g;\alpha)=\int_{I_i}fg_xdx-(\widetilde{f_\alpha}g^-)_{i+\frac12,y}+(\widetilde{f_\alpha}g^+)_{i-\frac12,y},\]

\[\mathcal{A}_{J_j}(f,g;\alpha)=\int_{J_j}fg_ydy-(\widetilde{f_\alpha}g^-)_{x,j+\frac12}+(\widetilde{f_\alpha}g^+)_{x,j-\frac12}.\]

§.§ Discrete energy law

Similar to the 1D case, we start by presenting the linear energy growth property of the exact solutions.

Let $E,~S,~T$ be the solutions to the equation (<ref>) under the periodic boundary condition, and define the two-dimensional energy as $\mathcal{E}(t)=\int_J\int_I E(x,y,t)^2+S(x,y,t)^2+T(x,y,t)^2
dxdy$, then for any $t$, the global stochastic energy satisfies the following
energy law

\begin{eqnarray}\label{2dEL}
\mathcal{E}(t)=\mathcal{E}(0)+2\int_0^t\int_J\int_I \lambda_2(S+T)-\lambda_1 E dW_sdxdy+(\lambda_1^2+2\lambda_2^2)Tr(Q)t,
\end{eqnarray}

and, after taking the expectation,

\begin{eqnarray}\label{linearpdelevel}
\mathbb{E}\Big(\mathcal{E}(t)\Big)=\mathcal{E}(0)+(\lambda_1^2+2\lambda_2^2)Tr(Q)t.
\end{eqnarray}

The proof is skipped here as it follows the same analysis as that of Theorem <ref>.

Let $E_h(\omega,x,y,t)$, $S_h(\omega,x,y,t)$ and $T_h(\omega,x,y,t)$ be the numerical solutions to the DG methods (<ref>) - (<ref>).

(a): For any $t\in[0,T]$, the numerical solutions satisfy the discrete energy law

\begin{eqnarray}\label{ccccccccxian}
&&\norm{E_h(\omega,x,y,t)}^2+\norm{S_h(\omega,x,y,t)}^2+\norm{T_h(\omega,x,y,t)}^2 \\
+\norm{T_h(x,y,0)}^2+(\lambda_1^2+2\lambda_2^2)Kt ,
\notag
\end{eqnarray}

\begin{align}\notag
&~~~~~~~~~~~~~~~~~~~~ =\norm{E_h(x,y,0)}^2+\norm{S_h(x,y,0)}^2+\norm{T_h(x,y,0)}^2+(\lambda_1^2+2\lambda_2^2)Kt.\label{2denergyy}
\end{align}

\begin{equation}
K=\sum_{i,j}\sum_{l=0}^{k^2+2k} \mu_{i,j}^l\sum^\infty_{m=1} \Big(\int_{J_j}\int_{I_i}\phi_{i,j}^l\sqrt{\gamma_m}e_mdxdy\Big)^2, \label{added2d2}
\end{equation}

where  $\{\phi^l_{i,j}, i,j=0,\cdots,k\}$ represents the set of Legendre basis over cell $J_j\times I_i$, and
$\mu_{i,j}^l=(\int_{J_j}\int_{I_i}(\phi_{i,j}^l )^2dxdy)^{-1}$.

(b) The constant $K$ is bounded by $(k^2+2k)Tr(Q)$ with $k$ being the polynomial degree of the DG methods.
Moreover, if there exists some constant $\alpha>0$ such that the series
$\sum^\infty_{m=1}\gamma_m (K_m)^\alpha<\infty$ with $K_m=\norm{\nabla e_m}_{{\infty}}$,
we can show that $K=Tr(Q)+O(h^\alpha)$, i.e., $K$ is an approximation of $Tr(Q)$ appearing in the continuous energy law (<ref>)-(<ref>).

(a):   By Itô's lemma, we know that

\begin{equation}\label{ito2d}
\end{equation}

By taking the test functions $\varphi=E_h$, $\psi=S_h$ and  $\phi=T_h$ in (<ref>)-(<ref>), and adding these equations, we have

\begin{eqnarray}\label{qqqq}
+\int_{I_i}(\widetilde{\Theta})_{x,j-\frac12}-(\widetilde{\Theta})_{x,j+\frac12}dxdt\notag \\
\end{eqnarray}

\[\Theta=\big(\frac12+\alpha_1\big)T_h^+E_h^-+\big(\frac12-\alpha_1\big)E_h^+T_h^-,
~~ \widetilde{\Theta}=\big(\frac12+\alpha_2\big)E_h^+S_h^-+\big(\frac12-\alpha_2\big)E_h^-S_h^+.\]

Following the same analysis as in Theorem <ref>, we can evaluate the terms
we have

\begin{align}
&\int_{J_j}\int_{I_i}(dS_h)^2dxdy= \int_{J_j}\int_{I_i}(dT_h)^2dxdy=\lambda_2^2 \sum_{l=0}^{k^2+2k} \mu_{i,j}^l\sum^\infty_{m=1} \Big(\int_{J_j}\int_{I_i}\phi_{i,j}^l\sqrt{\gamma_m}e_mdxdy\Big)^2dt, \notag \\
&\int_{J_j}\int_{I_i}(dE_h)^2dxdy =\lambda_1^2 \sum_{l=0}^{k^2+2k} \mu_{i,j}^l\sum^\infty_{m=1} \Big(\int_{J_j}\int_{I_i}\phi_{i,j}^l\sqrt{\gamma_m}e_mdxdy\Big)^2dt,   \label{added13}
\end{align}

By combining the equations (<ref>)-(<ref>), summing over all the cells and integrating in time from $0$ to $t$,
we can obtain

\begin{eqnarray}\label{cccccccc}
&&\norm{E_h(x,y,t)}^2+\norm{S_h(x,y,t)}^2+\norm{T_h(x,y,t)}^2 \\
2\int_0^t\int_{J}\int_{I}\lambda_2(T_h+S_h)-\lambda_1E_hdW_sdxdy+\norm{E_h(x,y,0)}^2+\norm{S_h(x,y,0)}^2+\norm{T_h(x,y,0)}^2+(\lambda_1^2+2\lambda_2^2)Kt. \notag
\end{eqnarray}

Note that  $\mathbb{E}\Big(\int_0^t\int_{J}\int_{I}\lambda_2(T_h+S_h)-\lambda_1E_hdW_sdxdy\Big)=0$ since it is an
Itô integral, therefore taking expectation of (<ref>) leads to the discrete energy law (<ref>).

(b): The estimate of the constant $K$ follows an exact same analysis as in the proof of Theorem <ref> (b), that is,

\[K\le (k^2+2k)Tr(Q), \]

and if  $\sum^\infty_{m=1}\gamma_m (K_m)^\alpha <\infty$, we can show that

\[K-Tr(Q)=O(h^\alpha).\]

If $W_t$ is the standard Brownian motion, it can be shown that $K= Tr(Q)$ (see Remark <ref>),
which means that the two-dimensional continuous energy law (<ref>) is exactly preserved by the proposed method.

§.§ Optimal error estimate

In this section we study the convergence rate of the DG scheme (<ref>)-(<ref>).
We start with defining the generalized Radau projection in $\mathbb{R}^2$ as

\begin{eqnarray}\label{2dproj}
\mathbb{P}^\alpha_x=\mathcal{P}^\alpha_x \otimes \mathcal{P}_y,~\mathbb{P}^\alpha_y=\mathcal{P}_x \otimes \mathcal{P}^\alpha_y,~
\mathbb{P}^{\alpha,\beta}= \mathcal{P}^\alpha_x \otimes \mathcal{P}^\beta_y,
\end{eqnarray}

where $\mathcal{P} $ is the $L^2$ projection.
The following lemmas are provided  in [21] and will be useful in our analysis.

For any function $w$, denote $\epsilon=\mathbb{P}^{\alpha,\beta}w-w$, with $\mathbb{P}^{\alpha,\beta}$ being the projection defined in (<ref>). For any $v\in \mathbb{V}_h^k$, there exists some constant $C$ such that

\[\Bigg|\sum_{i,j}\int_{J_j}\mathcal{A}_{I_i}(\epsilon,v,\alpha)dy\Bigg|\le
\Bigg|\sum_{i,j}\int_{I_i}\mathcal{A}_{J_j}(\epsilon,v,\alpha)dx\Bigg|\le Ch^{k+1}\norm{v}.\]

Let $\Pi$ be any projection defined in (<ref>). For any function $w(x,y)$, there exists some constant $C$, such that

\[\norm{\Pi w-w}\le Ch^{k+1}.\]

Now we turn to the main result on the error estimate of the DG methods.

Let $E_h, S_h, T_h\in\mathbb{V}_h^k$ be the numerical solutions given by the DG scheme (<ref>) - (<ref>),
and $E, T, S\in L^2(\Omega\times[0,T];H^{k+2})$ are strong solutions to (<ref>). With the initial conditions chosen as

\[E_h(x,y,0)=\mathbb{P}^{-\alpha_1,\alpha_2}E(x,y,0);~ S_h(x,y,0)=\mathbb{P}^{-\alpha_2}_yS(x,y,0);~ T_h(x,y,0)=\mathbb{P}_x^{\alpha_1}T(x,y,0),\]

there holds the following error estimates

\begin{equation}\label{2derroreq}
\norm{E-E_h}^2+\norm{S-S_h}^2+\norm{T-T_h}^2\le Ch^{2k+2},
\end{equation}

where the constant C denotes a generic positive constant independent of the spatial step sizes h.

Since both the numerical and exact solutions satisfy the equations (<ref>) - (<ref>),  we have the error equations

\begin{eqnarray}
\int_{J_j}\int_{I_i}d(E-E_h)\varphi dxdy&=&-\int_{J_j}\mathcal{A}_{I_i}(T-T_h,\varphi;\alpha_1)
\label{2ddgge1}\\
\int_{J_j}\int_{I_i}d(S-S_h)\psi dxdy &=&
\int_{I_i}\mathcal{A}_{J_j}(E-E_h,\psi;\alpha_2)dxdt ,
\label{2ddgge2}\\
\int_{J_j}\int_{I_i} d(T-T_h) \phi dxdy
&=& -\int_{J_j}\mathcal{A}_{I_i}(E-E_h,\phi;-\alpha_1)dydt.
\label{2ddgge3}
\end{eqnarray}

\begin{align*}
&\xi^E=\mathbb{P}^{-\alpha_1,\alpha_2}E-E_h, \quad \xi^S=\mathbb{P}_y^{-\alpha_2}S-S_h, \quad
\xi^T=\mathbb{P}_x^{\alpha_1}T-T_h,  \\
&\epsilon^E=\mathbb{P}^{-\alpha_1,\alpha_2}E-E, \quad
\epsilon^S=\mathbb{P}_y^{-\alpha_2}S-S, \quad \epsilon^T=\mathbb{P}_x^{\alpha_1}T-T,
\end{align*}

which leads to the decomposition of the error into two terms as

\[E-E_h=\xi^E-\epsilon^E,~~S-S_h=\xi^S-\epsilon^S,~~T-T_h=\xi^T-\epsilon^T.\]

By choosing the test functions as $\varphi=\xi^E,~\psi=\xi^S,~\phi=\xi^T$ in (<ref>)-(<ref>), and noting that

\[\int_{J_j}\mathcal{A}_{I_i}(\epsilon^T,\xi^E;\alpha_1)dydt=\int_{I_i}\mathcal{A}_{J_j}(\epsilon^S,\xi^E;-\alpha_2)dxdt=0\]

by the definition of the projections, we have

\begin{eqnarray}
&&\int_{J_j}\int_{I_i}d\xi^E\xi^E - d\epsilon^E\xi^Edxdy
+\int_{I_i}\mathcal{A}_{J_j}(\xi^S,\xi^E;-\alpha_2) dxdt, \label{2derror1}\\
&& \int_{J_j}\int_{I_i}d\xi^S\xi^S - d\epsilon^S\xi^Sdxdy
-\int_{I_i}\mathcal{A}_{J_j}(\epsilon^E,\xi^S;\alpha_2)dxdt,   \label{2derror2}\\
&& \int_{J_j}\int_{I_i}d\xi^T\xi^T - d\epsilon^T\xi^T
-\int_{J_j} \mathcal{A}_{I_i}(\xi^E,\xi^T;-\alpha_1)dydt
+\int_{J_j}\mathcal{A}_{I_i}(\epsilon^E,\xi^T;-\alpha_1)dydt.  \label{2derror3}
\end{eqnarray}

Summing up (<ref>)-(<ref>) and using integration by parts, we

\begin{align}
&\int_{J_j}\int_{I_i}d\xi^E\xi^E+d\xi^S\xi^S+d\xi^T\xi^T dxdy~~~~~~~~~~~~~~~~~~~~~~~~~~  \notag \\
&=    \int_{J_j}\int_{I_i}d\epsilon^E\xi^E+d\epsilon^S\xi^S+d\epsilon^T\xi^T    dxdy
-\bar{\Pi}_{x,j+\frac12}\Bigg)dxdt	\notag \\
+\int_{J_j}\mathcal{A}_{I_i}(\epsilon^E,\xi^T;-\alpha_1)dydt, \label{2dqqqq}
\end{align}

\[\Pi=\big(\frac12+\alpha_1\big)(\xi^T)^+(\xi^E)^-+\big(\frac12-\alpha_1\big)(\xi^T)^-(\xi^E)^+,~~
\bar{\Pi}=\big(\frac12+\alpha_2\big)(\xi^S)^-(\xi^E)^++ \big(\frac12-\alpha_2\big)(\xi^S)^+(\xi^E)^-.\]

By Itô's lemma, we have

\[d(\xi^E)^2=2d\xi^E\xi^E+d(\xi^E)^2,~~d(\xi^S)^2=2d\xi^S\xi^S+d(\xi^S)^2,~~d(\xi^T)^2=2d\xi^T\xi^T+d(\xi^T)^2.\]

Following the exact same analysis as shown in the proof of Theorem <ref>, we have

\[\int_{J_j}\int_{I_i}d(\xi^E)^2dxdy\le C\int_{J_j}\int_{I_i}(\mathbb{P}^{-\alpha_1,\alpha_2}(dW_t)-dW_t)^2dxdy,\]

\[\int_{J_j}\int_{I_i}d(\xi^S)^2dxdy\le C\int_{J_j}\int_{I_i}(\mathbb{P}^{-\alpha_2}_y(dW_t)-dW_t)^2dxdy,\]

\[\int_{J_j}\int_{I_i}d(\xi^T)^2dxdy\le C\int_{J_j}\int_{I_i}(\mathbb{P}^{\alpha_1}_x(dW_t)-dW_t)^2dxdy.\]

Therefore, by summing over all the cells $I_i\times J_j$ in (<ref>) and applying the superconvergence property in Lemma <ref>, we obtain

\begin{align}\label{2dwwwww}
\begin{aligned}
+ \int_{J}\int_{I}d\xi^E\xi^E+d\xi^S\xi^S+d\xi^T\xi^T dxdy\\
&  \le  \int_{J}\int_{I}d\epsilon^E\xi^E+d\epsilon^S\xi^S+d\epsilon^T\xi^T
\end{aligned}
\end{align}

Integrating (<ref>) over $t$,  applying the projection error property and Young's inequality, we can show that

\[\norm{\xi^E}^2+\norm{\xi^S}^2+\norm{\xi^T}^2\le C\int_0^t \norm{\xi^E(x,y,s)}^2+\norm{\xi^S(x,y,s)}^2+\norm{\xi^T(x,y,s)}^2ds+ Ch^{2k+2}.\]

Applying the Gronwall's inequality and combining with the optimal projection error yields the desired optimal error estimate (<ref>).

§.§ Multi-symplectic structure

Similar to the 1D case, to rewrite the two-dimensional stochastic Maxwell equations (<ref>) in the multi-symplectic form, we introduce the new variables such that

\[u_t=E,~v_t=S,~w_t=T,~P=T-\frac12u_x,~Q=S+\frac12u_y,~R=E-\frac12w_x+\frac12v_y.\]

and the system (<ref>) becomes

\begin{equation}\label{2dre}
\begin{cases}
\frac12u_y=Q-S,\\
\end{cases}
\end{equation}

We set $z=(T,S,E,w,v,u,P,Q,R)^T$, and define

\[M=\left(
\begin{array}{ccccccccc}
\end{array}
\right),\]

\[ K_1= \left(
\begin{array}{ccccccccc}
\frac12&0&0&0&0&0&0&0&0\\
\end{array}
\right), \qquad K_2=\left(
\begin{array}{ccccccccc}
\end{array}
\right),\]

\[S_1=PT+QS+RE-\frac{T^2+S^2+E^2}{2},~S_2=\lambda_1u-\lambda_2(w+v), \]

therefore, the two-dimensional stochastic Maxwell equations (<ref>) can be rewritten in the following multi-symplectic system

\begin{equation}\label{2dmulti}
Mdz+K_1z_xdt+K_2z_ydt=\nabla S_1(z)dt+\nabla S_2(z)dW_t.
\end{equation}

Its variational equation takes the form

\begin{equation}\label{2dvar}
\end{equation}

Let $U,~V$ be any solution to the variational equation (<ref>), and define $\omega=MU\cdot V, ~\kappa_x=K_1U\cdot V,~\kappa_y=K_2U\cdot V$, then we can derive the multi-symplectic conservation law given by

\[d\omega +\kappa_xdt+\kappa_ydt=0.\]

Since the new system (<ref>) is equivalent to the original model (<ref>), we can rewrite the proposed DG methods (<ref>)-(<ref>) into a consistent formulation for the system (<ref>).
For $E_h,~S_h,~T_h$ defined in (<ref>) - (<ref>), find $w_h,~v_h,~u_h\in\mathbb{V}_h^k$, such that for all
$\psi,~\bar{\psi},~\widetilde{\psi}\in\mathbb{V}_h^k$, it holds that

\begin{align}
&\int_{J_j}\int_{I_i}dw_h\psi dxdy=\int_{J_j}\int_{I_i}T_h\psi dxdy\, dt,\label{wh}\\
&\int_{J_j}\int_{I_i}dv_h\bar{\psi }dxdy=\int_{J_j}\int_{I_i}S_h\bar{\psi }dxdy\, dt,\label{vh}\\
&\int_{J_j}\int_{I_i}du_h\widetilde{\psi }dxdy=\int_{J_j}\int_{I_i}E_h\widetilde{\psi }dxdy\, dt.\label{uh}
\end{align}

Next, find $P_h,~Q_h,~R_h\in\mathbb{V}_h^k$ such that for any $\varphi,~\widetilde{\varphi},~\bar{\varphi}\in\mathbb{V}_h^k$, it holds that

\begin{eqnarray}
&&\int_{J_j}\int_{I_i}(P_h-T_h)\varphi dxdy=\frac12\int_{J_j}\Bigg(\int_{I_i}u_h\varphi_xdx-(\widehat{u}_h\varphi^-)_{i+\frac12,y}
+(\widehat{u}_h\varphi^+)_{i-\frac12,y}\Bigg)dy, \label{P}\\
+(\widehat{u}_h{\bar{\varphi}}^+)_{x,j-\frac12}\Bigg)dx, \label{Q}\\
&&\int_{J_j}\int_{I_i}(R_h-E_h)\widetilde{\varphi} dxdy=\frac12\int_{J_j}\Bigg(\int_{I_i}w_h\widetilde{\varphi}_xdx-(\widehat{w}_h\widetilde{\varphi}^-)_{i+\frac12,y}
&& ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~-\frac12\int_{I_i}\Bigg(\int_{J_j}v_h\widetilde{\varphi}_ydy-(\widehat{v}_h\widetilde{\varphi}^-)_{x,j+\frac12}
\end{eqnarray}

where the numerical fluxes are chosen as

\[(\widehat{u}_h)_{x_0,y}=(\{u_h\}-2m_1[u_h])_{x_0,y},~(\widehat{w}_h)_{x_0,y}=(\{w_h\}-2n_1[w_h])_{x_0,y},~~m_n,~n_2\in\mathbb{R} ~\text{with}~m_1-n_1=\alpha_1,\]

\[(\widehat{u}_h)_{x,y_0}=(\{u_h\}+2m_2[u_h])_{x,y_0},~(\widehat{v}_h)_{x,y_0}=(\{v_h\}+2n_2[u_h])_{x,y_0},~~m_n,~n_2\in\mathbb{R} ~\text{with}~m_2-n_2=\alpha_2.\]

By combining the derivative of (<ref>)-(<ref>) with the equations (<ref>)-(<ref>), we obtain

\begin{align} \notag
&\int_{J_j}\int_{I_i}dP_h\phi dxdy=-\frac12\int_{J_j}\Bigg(\int_{I_i}E_h\phi_xdx-(\widehat{E}_h\phi^-)_{i+\frac12,y}
+(\widehat{E}_h\phi^+)_{i-\frac12,y}\Bigg)dy dt
\int_{J_j}\int_{I_i}\phi dW_tdxdy, \\
\int_{J_j}\int_{I_i}\bar{\phi} dW_tdxdy, \notag \\
+(\widehat{T}_h\widetilde{\phi}^+)_{i-\frac12,y}\Bigg)dy \label{PQR3}\\
& ~~~~~~~~~~~~~~~~~~~~~~~+\frac12\int_{I_i}\Bigg(\int_{J_j}S_h\widetilde{\phi}_ydy-(\widehat{S}_h\widetilde{\phi}^-)_{x,j+\frac12}
\int_{J_j}\int_{I_i}\widetilde{\phi} dW_tdxdy,  \notag
\end{align}

where the numerical fluxes are

\[(\widehat{T}_h)_{x_0,y}=(\{T_h\}+2m_1[T_h])_{x_0,y},~~(\widehat{E}_h)_{x_0,y}=(\{E_h\}+2n_1[E_h])_{x_0,y},\]

\[(\widehat{E}_h)_{x,y_0}=(\{E_h\}-2n_2[E_h])_{x,y_0},~~(\widehat{S}_h)_{x,y_0}=(\{S_h\}-2m_2[S_h])_{x,y_0}.\]

Combining (<ref>)-(<ref>), we have derived the equivalent formulation of the DG scheme (<ref>)-(<ref>) for the new system (<ref>):

find $z_h\in
(\mathbb{V}^k_h)^9$, such that for all $\bm{\varphi}\in(\mathbb{V}^k_h)^9$, we have

\begin{eqnarray}
&&\int_{J_j}\int_{I_i}Mdz_h\cdot \bm{\varphi} dxdy-\int_{J_j}\Big(\int_{I_i}K_1z_h\cdot\bm{\varphi}_xdx
\big(\widehat{K_1z_h}\cdot\bm{\varphi}^+\big)_{i-\frac12,y}\Big)dydt \notag \\
\big(\widehat{K_2z_h}\cdot\bm{\varphi}^+\big)_{x,j-\frac12}\Big)dxdt\notag \\
&&=\int_{J_j}\int_{I_i}\nabla S_1(z_h)\cdot\bm{\varphi} dxdydt+\int_{J_j}\int_{I_i}\nabla S_2(z_h)\cdot\bm{\varphi} dW_tdxdy, \label{2ddg}
\end{eqnarray}

where  $(\widehat{K_1z_h})_{x_0,y}=(K_1\{z_h\}+A_1[z_h])_{x_0,y}$ , and

\[A_1=\left(
\begin{array}{ccccccccc}
\end{array}\right), \qquad A_2=\left(
\begin{array}{ccccccccc}
\end{array}
\right),\]

Applying the exterior derivative to this scheme (<ref>) leads to the variation equation

\begin{eqnarray}
&&\int_{J_j}\int_{I_i}MdZ_h\cdot \bm{\varphi} dxdy-\int_{J_j}\Big(\int_{I_i}K_1Z_h\cdot\bm{\varphi} _xdx
-\big(\widehat{K_1Z_h}\cdot\bm{\varphi} ^-\big)_{i+\frac12,y}+
\big(\widehat{K_1Z_h}\cdot\bm{\varphi} ^+\big)_{i-\frac12,y}\Big)dydt \label{2ddgvar} \\
&&-\int_{I_i}\Big(\int_{J_j}K_2z_h\cdot\bm{\varphi} _ydy
-\big(\widehat{K_2Z_h}\cdot\bm{\varphi} ^-\big)_{x,j+\frac12}+
\big(\widehat{K_2Z_h}\cdot\bm{\varphi} ^+\big)_{x,j-\frac12}\Big)dxdt=\int_{J_j}\int_{I_i}\nabla^2S_1Z_h\cdot\bm{\varphi}  dxdydt, \notag
\end{eqnarray}

Let  $U_h, V_h$ be any solutions to the variational equation (<ref>), we have the semi-discrete version of
the multi-symplectic conservation laws

\begin{eqnarray}
\int_{J_j}\int_{I_i} d(MU_h\cdot V_h)dx-\int_{J_j}\Big(\mathcal{F}_{K_1}(U_h,V_h)_{i+\frac12,y}
-\mathcal{F}_{K_1}(U_h,V_h)_{i-\frac12,y}\Big)dydt \notag \\
-\mathcal{F}_{K_2}(U_h,V_h)_{x,j-\frac12}\Big)dxdt=0. \label{2ddismul}
\end{eqnarray}

The proof follows the same idea as that of Theorem <ref>, and is skipped here to save space.

§ SYMPLECTIC TIME DISCRETIZATION

In this section, we discuss the symplectic temporal discretization for the semi-discrete DG methods presented in the previous section.

In the one-dimensional case, we can set $\eta_h|_{I_j}=\sum^k_{l=1}\eta_j^l\varphi_j^l$,
$u_h|_{I_j}=\sum^k_{l=1}u_j^l\varphi_j^l$, with the set $\{\varphi_j^l\}$ being the basis of $V_h^k$, and introduce the notations

\[p=(\eta_1^0,\cdots,\eta_1^k,\eta_2^0,\cdots,\eta_2^k,\cdots,\eta_J^k)^T,
\quad q=(u_1^0,\cdots,u_1^k,u_2^0,\cdots,u_2^k,\cdots,u_J^k)^T.\]

In two-dimensional case, consider

\[E_h=\sum^k_{l=1} E_{i,j}^l\varphi_{i,j}^l,\quad  S_h=\sum^k_{l=1} S_{i,j}^l\varphi_{i,j}^l, \quad T_h=\sum^k_{l=1} T_{i,j}^l\varphi_{i,j}^l,\]

where $\{\varphi_{i,j}^l\}$ is the basis of $\mathbb{V}_h^k$.  Define

\begin{align*}
&E_{i,j}=(E_{i,j}^0, E_{i,j}^1, \cdots, E_{i,j}^{k^2+2k})^T, \\
&\bm{E}_h= \left(E_{1,1}, E_{2,1}, \cdots, E_{I,1}, E_{1,2}, \cdots, E_{I,J} \right)^T,
\end{align*}

and similarly for $\bm{S}_h$, $\bm{T}_h$, and then introduce the notations

\[p=\bm{E}_h, \qquad q=\left( \bm{S}_h , ~ \bm{T}_h\right)^T.\]

With these notations, either the one-dimensional scheme (<ref>)-(<ref>) or the one-dimensional method (<ref>)-(<ref>) can be simplified into the following stochastic differential equation:

\begin{align} \label{tq}
\begin{aligned}
dp=Aq+Ld\mathcal{B}_t, \\
\end{aligned}
\end{align}

where $A,~B,~L,~N$ are some constant matrices which may take different values in 1D or 2D setting, and $\mathcal{B}_t$ is a Brownian motion. The system (<ref>) will be considered in this section.

§.§ Symplectic Euler

We let $0=t_0\leq t_1 \leq \cdots \leq t_N=T$ be a partition of the time interval $[0,T]$. By setting $\tau =t_{k+1}-t_k$, and $\Delta \mathcal{B}=\mathcal{B}_{t_{k+1}}-\mathcal{B}_{t_k}$, the symplectic Euler methods for the system (<ref>) are given by

\begin{align}\label{Eq}
& p^{k+1}=p^k+\tau Aq^k+L\Delta \mathcal{B},  \quad
& q^{k+1}=q^k+\tau Bp^{k+1}+N\Delta \mathcal{B}.
\end{align}

Following the studies in [20], we have the following result:

Symplectic Euler method (<ref>) preserves the symplectic structure, and has the first mean-square order of

§.§ Partitioned Runge-Kutta method

Consider the two-stage PRK methods [20] for the system (<ref>) of the form

\begin{eqnarray}
\begin{aligned}\label{prk}
&\mathcal{Q}_1=q^k+N\Big(J_k+\frac{1}{\sqrt{2}}\Delta \mathcal{B}\Big)\\
&\mathcal{P}_1=p^k+\frac14\tau A\mathcal{Q}_1
+L\Big(J_k+\frac{1}{2\sqrt3}\Delta \mathcal{B}\Big),\\
&\mathcal{Q}_2=q^k+\frac23\tau B\mathcal{P}_1+
N\Big(J_k-\frac{1}{3\sqrt2}\Delta \mathcal{B}\Big),\\
&\mathcal{P}_2=p^k+\tau \Big(\frac14A\mathcal{Q}_1+\frac34 A\mathcal{Q}_2\Big)
&p^{k+1}=p^k+L\Delta\mathcal{B}+\tau \Big(\frac14 A\mathcal{Q}_1+\frac34
\Big(\frac23B\mathcal{P}_1+\frac13B\mathcal{P}_2\Big),
\end{aligned}
\end{eqnarray}

$$\tau =t_{k+1}-t_k, \qquad \Delta \mathcal{B}=\mathcal{B}_{t_{k+1}}-\mathcal{B}_{t_k}, \qquad J_k=\frac{1}{\tau}\int_{t_k}^{t_{k+1}}\mathcal{B}_s-\mathcal{B}_{t_k} ds.$$

In order to analyze the convergence rate of the PRK methods (<ref>), we recall the following theorem in <cit.>.

Consider a general system of stochastic differential equation

\begin{equation}\label{SDE}
\end{equation}

where $X$, $a$ and $b_r$ are column vectors defined on $t\in [t_0, T]$, and $\mathcal{B}_r$ are independent
standard Brownian motions, and let $\overline{X}_{t,x}(t+\tau)$ be a one-step
approximation satisfying that for any $t_0\le t\le T-\tau$,

\begin{align}\label{th11}
&\Bigg|\mathbb{E}\Big(X_{t,x}(t+\tau)-\overline{X}_{t,x}(t+\tau)\Big)\Bigg|\le K(1+|x|^2)^{1/2}\tau^{p_1}, \\
K(1+|x|^2)^{1/2}\tau^{p_2}. \label{th22}
\end{align}

Let $p_2\ge\frac12,~~p_1\ge p_2+\frac12$, then for any $N$ and $k=0,1,\cdots,N$,
the following result holds:

\begin{eqnarray}\label{thcon}
\Bigg(\mathbb{E}\Big|X_{t_0,X_0}(t_k)-\overline{X}_{t_0,X_0}(t_k)\Big|^2\Bigg)^{1/2}\le
\end{eqnarray}

The following Lemma is a direct result of Theorem <ref>, and will be useful in our analysis.

Let the one-step approximation $\overline{X}_{t,x}(t+\tau)$ satisfies the
conditions of Theorem <ref>, and suppose that another one step
approximation $\widetilde{X}_{t,x}(t+\tau)$ satisfies

\begin{align}\label{lem11}
& \Bigg|\mathbb{E}\Big(\widetilde{X}_{t,x}(t+\tau)-\overline{X}_{t,x}(t+\tau)\Big)\Bigg|
=O(\tau^{p_1}), \\
& \Bigg(\mathbb{E}\Big|\widetilde{X}_{t,x}(t+\tau)-\overline{X}_{t,x}(t+\tau)\Big|^2\Bigg)^{1/2}
=O(\tau^{p_2}), \label{lem22}
\end{align}

with the same $p_1$ and $p_2$ as in Theorem <ref>, then the mean-square order of accuracy for $\widetilde{X}_{t,x}(t+\tau)$
is also $p_2-1/2$, same as $\overline{X}_{t,x}(t+\tau)$.

We have the following main result on the convergence rate of the PRK method (<ref>).

The explicit PRK method (<ref>) for the system (<ref>) preserves symplectic structure and has the mean-square order of 2.

The preservation of symplectic structure is shown in <cit.>, here we only prove the correct convergence order.
According to <cit.>, the following one-step approximation for the system (<ref>) has the mean-square order $2$:

\begin{eqnarray}
\overline{p}^{k+1}=p^k+\tau A{q}^k+\frac12
\tau^2AB{p}^k+\tau AN J_k+L\Delta\mathcal{B},\label{2nd1}\\
\overline{q}^{k+1}={q}^k+\tau Bp^k+\frac12 \tau^2BA q^k+\tau BL J_k+N\Delta\mathcal{B},\label{znd2}
\end{eqnarray}

which will be used as a reference method.

The PRK method (<ref>) can be rewritten as

\[p^{k+1}=p^k+\tau Aq^k+\frac12\tau^2 ABp^k+\tau AN J_k+L\Delta\mathcal{B}+\frac12\tau^2 ABL\Big(J_k+\frac{1}{2\sqrt3}\Delta \mathcal{B}\Big)

Notice that

\[\mathbb{E}J_k=\mathbb{E}{\Delta \mathcal{B}}=0,\qquad \mathbb{E}{(J_k)^2}=\frac{\tau}{3},\qquad \mathbb{E}(\Delta \mathcal{B})^2=\tau,\]

which, by Young's inequality, leads to

\[\mathbb{E}\Big(J_k+\frac{1}{2\sqrt3}\Delta \mathcal{B}\Big)^2\le C\mathbb{E}(J_k)^2+C\mathbb{E}(\Delta \mathcal{B})^2=O(\tau).\]

Therefore, we have

\[\Bigg|\mathbb{E}\Big(p^{k+1}-\overline{p}^{k+1}\Big)\Bigg|=O(\tau^3),\]

\begin{eqnarray*}
\mathbb{E}\Big(p^{k+1}-\overline{p}^{k+1}\Big)^2=\mathbb{E}\Big(\frac12\tau^2 AB\Big(J_k+\frac{1}{2\sqrt3}\Delta \mathcal{B}\Big)
\end{eqnarray*}

which leads to

\[\Bigg(\mathbb{E}\Big(p^{k+1}-\overline{p}^{k+1}\Big)^2\Bigg)^{1/2}=O(\tau^{5/2}).\]

Following a similar approach, we also have

\[\Bigg|\mathbb{E}\Big(q^{k+1}-\overline{q}^{k+1}\Big)\Bigg|=O(\tau^3), \qquad
\Bigg(\mathbb{E}\Big(q^{k+1}-\overline{q}^{k+1}\Big)^2\Bigg)^{1/2}=O(\tau^{5/2}).\]

Therefore, by applying the result of Lemma <ref>, the PRK method (<ref>) has mean-square order of $2$.

The PRK scheme (<ref>) was shown in [20] to have the mean-square order of 3/2 for general system. For this linear system (<ref>), it can be shown to have second order mean-square convergence rate.

§ NUMERICAL EXPERIMENTS

The numerical results of the DG scheme is presented in this section to demonstrate the performance of the proposed methods.
We consider the DG method with various polynomial degree $k$ as the spatial discretization, and utilize the PRK method (<ref>) for temporal discretization.  The accuracy tests are provided for both 1D case and 2D case to demonstrate the convergence rate of the methods, and the linear growth of discrete energy is also studied for both cases.

§.§ Accuracy test

For the one-dimensional system (<ref>) with periodic boundary conditions, we set $\lambda_1=\lambda_2=1$, and choose $W_t=\mathcal{B}_t$, which is the standard Brownian motion, so that one exact solution to (<ref>) takes the form

\begin{equation}\label{exact}
\begin{cases}
\eta=\sin(x-t)+\cos(x+t)-\mathcal{B}_t&\\
\end{cases}
\end{equation}

The computational domain is $[0,2\pi]$ and the final time is set to $T=3$. Initial conditions for $\eta(x,0)$ and $u(x,0)$ are
obtained by letting $t=0$ in (<ref>).  We use Nx and Nt to denote number of spatial steps and temporal steps respectively.

Table <ref> and Table <ref> demonstrates the order of convergence rate for the case $k=1$ and $k=2$ respectively, where $\Delta t$ is chosen to be small enough to ensure that the spatial error dominates. Under both cases, we can observe the optimal convergence orders, i.e., $(k+1)$-th order of accuracy, which is consistent with the result in Theorem <ref> for the semi-discrete method. Note that the second order temporal discretization is used for both cases, therefore, one would expect a second order convergence even coupled with third order DG spatial discretization. This second order convergence is observed in Table <ref>, when larger $\Delta t$ is used.

Numerical error and convergence rates of 1D case when $k=1$.

Nx   Nt   $\norm{e_u}$   order   $\norm{e_\eta}$   order

20   600   0.01295   0   0.01725   0

40   1200   3.201E-3   2.0162   4.281E-3   2.0106

80   2400   7.851E-4   2.0276   1.083E-3   1.9829

160   4800   1.958E-4   2.0036   2.714E-4   1.9966

Numerical error and convergence rates of 1D case when $k=2$.

Nx   Nt   $\norm{e_u}$   order   $\norm{e_\eta}$   order

20   600   3.176E-4   0   4.407E-4   0

40   1200   3.958E-5   3.0045   5.486E-5   3.0061

80   2400   4.920E-6   3.0081   6.874E-6   2.9964

160   4800   6.144E-7   3.0013   8.645E-7   2.9912

Numerical error and convergence rates of 1D case when $k=2$ with larger $\Delta t$.

Nx   Nt   $\norm{e_u}$   order   $\norm{e_\eta}$   order

20   60   5.6019E-4   0   5.3547E-4   0

40   120   1.2175E-4   2.2020   9.5595E-5   2.4858

80   240   2.9228E-5   2.0585   2.0804E-5   2.2000

160   480   7.2337E-6   2.0145   4.9803E-6   2.0626

Next, we consider the two-dimensional stochastic Maxwell equations (<ref>) with periodic boundary conditions.
Set $\lambda_1=\lambda_2=1$, and choose $W_t=\mathcal{B}_t$, then the exact solution to (<ref>) takes the form

\begin{equation}\label{2dexact}
\begin{cases}
\end{cases}
\end{equation}

The spatial domain is set to be $[0,2\pi]^2$, and the final stopping time is taken to be $T=1$. The initial conditions of $E, ~S,~T$ can be obtained by setting $t=0$ in the exact solutions (<ref>).  We use Nx and Nt to denote number of spatial steps and temporal steps respectively.
The numerical errors and the corresponding convergence rates are shown in Table <ref> for $k=1$ and in Table <ref> for $k=2$.
Under both cases, we can observe the optimal convergence orders, i.e., $(k+1)$-th order of accuracy, which matches
the analysis in Theorem <ref> for the semi-discrete method.

Numerical error and convergence rates of 2D case when $k=1$.

Nx   Ny   Nt   $\norm{e_E}$   order   $\norm{e_S}$   order   $\norm{e_T}$   order

20   20   200   0.02715   0   0.01749   0   0.01749   0

40   40   400   6.617E-3   2.0366   4.466E-3   1.9696   4.466E-3   1.9696

80   80   800   1.634E-3   2.0176   1.129E-3   1.9839   1.129E-3   1.9839

160   160   1600   4.061E-4   2.0086   2.839E-4   1.9918   2.839E-4   1.9918

Numerical error and convergence rates of 2D case when $k=2$.

Nx   Ny   Nt   $\norm{e_E}$   order   $\norm{e_S}$   order   $\norm{e_T}$   order

20   20   200   8.529E-4   0   5.899E-4   0   5.899E-4   0

40   40   400   9.980E-5   3.0952   7.061E-5   3.0628   7.061E-5   3.0628

80   80   800   1.253E-5   2.9939   9.119E-6   2.9529   9.119E-6   2.9529

160   160   1600   1.521E-6   3.0422   1.137E-6   3.0038   1.137E-6   3.0038

§.§ Averaged energy growth

The discrete energy law satisfied by the numerical solutions $u_h$ and $\eta_h$ was studied in
Theorem <ref> for the one-dimensional system. Under the same setup as in the previous example, we simulate the
system up to the final stopping time $T=3$, with the time step size $\Delta t=0.0075$ and 80 space meshes. In this test, we run the simulations with $1000$ samples, and take the average to approximate the expectation and compute the averaged energy. Figure <ref>
shows the time history of the averaged energy of our numerical solutions with different noise size $\lambda_{1,2}$,
from which we can observe that the averaged energy is almost linear with respect to time. Note that when $\lambda_1=\lambda_2=0$,
the global energy is preserved exactly on the discrete level. When the noise sizes $\lambda_{1,2}$ increase, the
linear growth rate of the discrete energy also increases. The slopes of the lines in Figure <ref> (b-d) can be computed via least square fitting, and are approximately 0.1228, 3.0154, 12.4137, respectively. They are proportional to $\lambda_1^2+\lambda_2^2$, which is consistent with the result in Theorem <ref>.

Averaged energy with different sizes of noise for 1D case.

Similarly, Theorem <ref> studied the discrete energy law satisfied by the numerical solutions $E_h$, $S_h$ and $T_h$ for the two-dimensional system.
For this example we simulate the
system up to the final stopping time $T=1$, with the time step size $\Delta t=0.0025$ and $80\times 80$ space meshes.
In this test, we run the simulations with $500$ samples, and take the average to approximate the expectation and compute the averaged energy. Figure <ref>
shows the time history of the averaged energy of our numerical solutions with different noise size $\lambda_{1,2}$,
from which we can observe that the averaged energy is almost linear with respect to time. Note that when $\lambda_1=\lambda_2=0$,
the global energy is preserved exactly on the discrete level. When the noise sizes $\lambda_{1,2}$ increase, the
linear growth rate of the discrete energy also increases.
The slopes of the lines in Figure <ref> (b-d) can be computed via least square fitting, and are approximately 1.2420, 29.9783, 117.2386, respectively. They are proportional to $\lambda_1^2+\lambda_2^2$, which is consistent with the result in Theorem <ref>.

Averaged energy with different sizes of noise for 2D case.

§.§ Test with noises of various sizes

In this example, we consider the two-dimensional system (<ref>), and use the initial conditions studied in [1]:

\begin{eqnarray*}
\begin{cases}
E(x,y,0)=\sin(3\pi x)\sin(4\pi y),\\
S(x,y,0)=-\frac45\cos(3\pi x)\cos(4\pi y),\\
T(x,y,0)=-\frac35\sin(3\pi x)\cos(4\pi y),
\end{cases}
\end{eqnarray*}

with $I\times J=[0,\frac23]\times[0,\frac12]$. The final stopping time is set as $T=1$.
Following the definition (<ref>), we construct the Wiener process as

\begin{eqnarray}\label{constW}
W_t=\sum^\infty_{m,n=1}2\sqrt{\frac{3}{m^3+n^3}}\sin\Big(\frac32m\pi x\Big)\sin\Big(2n\pi y\Big)\mathcal{B}_m(t),
\end{eqnarray}

and truncate the sum (<ref>) by taking the sum over $m,n$ from $1$ to $50$.

In this example, the mesh sizes $\Delta x=\Delta y=0.0083$ and time step size $\Delta t=0.00083$ are considered.
In order to show the effect of noise with various sizes on the numerical solution, we run the simulations with six sets of parameters: $\lambda_{1,2}=0.1, ~0.2, ~0.3, ~0.5, ~0.7,~1$. The contour plots of the numerical solution $S_h$ with these choices of noises are shown in Figure <ref>. The 3D plots of $T_h$ are also provided in Figure <ref>. It can be observed from these figures that the perturbation of the numerical solutions becomes more and more obvious as the size of the noise grows.

Contour plots of $S_h(x,y)$ with different sizes of noise for the test in Section <ref>. $8$ uniformly spaced contour lines  within the range $[-0.7, 0.7]$ are used.

3D plots of $T_h(x,y)$ with different sizes of noise for the test in Section <ref>.

§ CONCLUSION REMARK

In this paper we have developed and analyzed the DG scheme for the one- and two-dimensional
stochastic Maxwell equations with additive noise. The proposed methods are shown to satisfy the discrete form of the stochastic
energy linear growth property. The optimal error estimate of the semi-discrete methods is also proven analytically.
By introducing auxiliary variables, we also rewrite the stochastic Maxwell equations into the multi-symplectic structure and demonstrate that the proposed DG methods preserve the multi-symplectic structure. The semi-discrete methods are then combined with the symplectic Euler or PRK temporal discretization methods. Numerical experiments are provided to test the performance of the resulting methods, and optimal error estimates and linear growth of the discrete energy can be observed for all cases.

[1]
C. Chen, J. Hong, L. Zhang. Preservation of physical properties of stochastic Maxwell equations with additive noise via stochastic multi-symplectic methods,
Journal of Computational Physics, 306, 500-519 (2016).

[2]
C. Chen, J. Hong, L. Ji, Runge-Kutta semidiscretizations for stochastic Maxwell equations with additive
noise. SIAM Journal on Numerical Analysis, 57, 702-727 (2019).

[3]
B. Cockburn, S. Hou, C.-W. Shu.
The Runge-Kutta local projection discontinuous Galerkin finite
element method for conservation laws IV: the multidimensional case, Mathematics of Computation, 54, 545-581 (1990).

[4]
B. Cockburn, G. Kanschat, I. Perugia, D. Schötzau, Superconvergence of the local discontinuous Galerkin method for elliptic problems on Cartesian
grids, SIAM Journal on Numerical Analysis, 39, 264-285 (2001).

[5]
B. Cockburn, G. Karniadakis, C.-W. Shu.
The development of discontinuous galerkin methods.
In B. Cockburn, G. Karniadakis, and C.-W. Shu, editors, Discontinuous Galerkin Methods: Theory, Computation and Applications, pages
3–50. Lecture Notes in Computational Science and Engineering, Part I:
Overview, volume 11, Springer, 2000.

[6]
B. Cockburn, S.-Y. Lin, C.-W. Shu,
TVB Runge-Kutta local projection discontinuous Galerkin
finite element method for conservation laws III: one dimensional systems,
Journal of Computational Physics, 84, 90-113 (1989).

[7]
B. Cockburn, C.-W. Shu,
TVB Runge-Kutta local projection discontinuous Galerkin
finite element method for conservation laws II: general framework, Mathematics of Computation, 52, 411-435 (1989).

[8]
D. Cohen, J. Cui, J. Hong, L. Sun, Exponential integrators for stochastic Maxwell's equations driven by Itô
noise, Journal of Computational Physics, 410, 109382 (2020).

[10]
J. Hong, L. Ji, L. Zhang. A stochastic multi-symplectic scheme for stochastic Maxwell equations with additive noise,
Journal of Computational Physics, 268, 255-268 (2014).

[11]
J. Hong, L. Ji, L. Zhang, J. Cai. An energy-conserving method for stochastic Maxwell equations with multiplicative
noise, Journal of Computational Physics, 351, 216-229 (2017).

[12]
S. Jiang, L. Wang, J. Hong. Stochastic multi-symplectic integrator for stochastic nonlinear Schrödinger
equation,  Communications in Computational Physics, 14, 393-411 (2013).

[13]
P. Kloeden, E. Platen, Numerical Solution of Stochastic Differential
Equations, 3rd ed., Vol. 23, Applications in Mathematics, Stochastic Modelling
and Applied Probability, Springer-Verlag, Berlin, 1999.

[14]
Y. Li, C.-W. Shu, S. Tang. A discontinuous Galerkin method for stochastic conservation
laws, SIAM Journal on Scientific Computing, 42, A54-A86 (2020).

[15]
Y. Li, C.-W. Shu, S. Tang. An ultra-weak DG method with IMEX time-marching for generalized stochastic KdV
equations, Journal of Scientific Computing, 82, 61 (2020).

[16]
Y. Li, S. Wu, Y. Xing. Finite element approximations of a class of nonlinear stochastic wave equation
with multiplicative noise, submitted.

[17]
G. Liu, Stochastic wave propagation in Maxwell's equations, Journal of
Statistical Physics, 158, 1126-1146 (2015).

[18]
K.B. Liaskos, I.G. Stratis, A.N. Yannacopoulos. Stochastic integrodifferential equations in Hilbert spaces
with applications in electromagnetics,
Journal of Integral Equations and Applications, 22, 559-590 (2010).

[19]
G.N. Milstein, Numerical Integration of Stochastic Differential
Equations, Ural State University, Sverdlovsk. English translation by Kluwer
Academic Publishers, Mathematics and its Applications, v.313, 1995.

[20]
G.N. Milstein, Yu. M. Repin, M.V. Tretyakov. Symplectic integration of Hamiltonian systems with additive noise,
SIAM Journal on Numerical Analysis, 39, 2066-2088 (2002).

[21]
X. Meng, C.-W. Shu, B. Wu, Optimal error estimates for discontinuous Galerkin methods based on upwind-biased fluxes for linear hyperbolic
equations, Mathematics of Computation, 85, 1225-1261 (2016).

[22]
Z. Sun, Y. Xing. On structure-preserving discontinuous Galerkin methods for Hamiltonian partial
differential equations: Energy conservation and multi-symplecticity,
Journal of Computational Physics, 419, 109662 (2020).

[23]
Z. Sun, Y. Xing. Optimal error estimates of discontinuous Galerkin methods with generalized fluxes
for wave equations on unstructured meshes,
Mathematics of Computation, 90, 1741-1772 (2021).

[24]
G.N. Ord. A stochastic model of Maxwell's equations in 1+1 dimensions,
International Journal of Theoretical Physics,  35, 263-266 (1996).

[25]
W.H. Reed, T. Hill, Triangular mesh methods for the neutron transport
equation, Technical report, Los Alamos Scientific Lab., N. Mex. (USA), 1973.

[26]
S.M. Rytov,  Y.A. Kravtsov,  V.I. Tatarskii. Principles of Statistical Radiophysics: Elements and Random Fields
3, Springer, Berlin, 1989.

[27]
L. Zhang, C. Chen, J. Hong, L. Ji. A review on stochastic multi-symplectic methods for
stochastic Maxwell equations, Communications on Applied Mathematics and
Computation 1, 467-501 (2019).