{'fraction_non_alphanumeric': 0.04873985324873276, 'fraction_numerical': 0.01938959980149587, 'mean_word_length': 4.291072768192048, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Discontinuous Galerkin methods for stochastic Maxwell equations with multiplicative noise

In this paper we propose and analyze finite element discontinuous Galerkin methods for the one- and two-dimensional stochastic Maxwell equations with multiplicative noise. The discrete energy law of the semi-discrete DG methods were studied. Optimal error estimate of the semi-discrete method is obtained for the one-dimensional case, and the two-dimensional case on both rectangular meshes and triangular meshes under certain mesh assumptions. Strong Taylor 2.0 scheme is used as the temporal discretization. Both one- and two-dimensional numerical results are presented to validate the theoretical analysis results.


Introduction
Stochastic Maxwell equations are commonly used to model microscopic origins of randomness in electromagnetism. The concept of stochastic Maxwell equations was firstly introduced by Rytov et al. in [23] to describe the fluctuations of an electromagnetic field. In [21], Ord et al. studied the Mark Kac random walk model and modified it into the Max field equations in 1+1 dimensions. Such model describes most of the telegraph equations, and the author's modification of the Kac model constructed a strong connection between the telegraph and Maxwell equations. In [12] Horsin et al. applied an abstract approach and a constructive approach by generalizing the Hilbert uniqueness method to analyze the approximate controllability of the stochastic Maxwell equations. Furthermore in [19] the deterministic and stochastic integrodifferential equations in Hilbert spaces were studied and the well-posedness for the Cauchy problem of the integrodifferential equation were analyzed. Such results were motivated from mathematical modeling of electromagnetics fields in complex random media.
Numerical methods are often used to solve Maxwell equations with various forms of stochasticity. In [1], Benner et al. studied the time-harmonic Maxwell's equations with some uncertainty in material parameters. They compared stochastic collocation and Monte Carlo simulation, and computed a reduced model in order to lower the computational cost. In [13] Jung considered the wave and Maxwell equations with fluctuations by random change in media parameters. The author used polynomial chaos Galerkin projections to develop the evolution of probability distribution function. Later in [14] Jung et al. studied two-dimensional transverse magnetic Maxwell equations with multiple random interfaces. They applied the polynomial chaos projection method, and computed the stochastic and deterministic part separately. Furthermore, stochastic collocation methods for metamaterial Maxwell's equations are studied by Li et al. in [16]. They considered the equations with random coefficient and random initial conditions. They also developed regularity analysis for stochastic metamaterial Maxwell's equations.
In the recent years, there have been many studies on various numerical methods for the stochastic Maxwell equations with either additive or multiplicative noises. In [11], Hong et al. studied the the stochastic Maxwell equations driven by a multiplicative noise, which are given as follows: where H represents the magnetic field, E stands for the electric field, λ is the scale of the noise, • denotes the stochastic integral in Stratonovich sense, and W is a space time mixed Wiener process. They studied the multi-symplectic structure and the energy conservation law for (1.1), and developed a fully discrete numerical method which conserves both multi-symplecticity and energy in the discrete level. Furthermore, Cohen et al. in [9] analyzed the general form of stochastic Maxwell equations with multiplicative noise, and constructed an exponential integrator which has a general mean square convergence order 0.5 in time, and first order temporal rate under some assumptions. Chen et al. in [3] applied a semi-implicit scheme for the model under a general setting and showed that the proposed method has the mean-square order of 0.5 in time. In [26], Zhang et al. presented a nice review article to summarize numerical methods for different kinds of stochastic Maxwell equations with both additive noise and multiplicative noise, including (1.1) studied in [11], and the model dE − ∇ × Hdt = −J e (t, x, E, H)dt − J r e (t, x, E, H)dW, µdH + ∇ × Edt = −J m (t, x, E, H)dt − J r m t, x, E, H)dW, (1.2) which was studied in [3]. The properties of stochastic Maxwell equations are also provided in that paper. Hong et al. in [10] studied the stochastic wave equation and developed numerical schemes that preserve the averaged energy evolution law. Both the compact finite difference method and the interior penalty discontinuous Galerkin (DG) finite element methods were proposed. Finite element approximations of a class of nonlinear stochastic wave equations with multiplicative noise were recently investigated in [18].
The high order DG finite element methods are considered in this paper. The DG method is a class of finite element methods that uses discontinuous piecewise polynomials as the basis functions, and was first introduced by Reed and Hill in [22] to solve linear transport equation. In the early 1990s, Cockburn et al. studied the extension of DG methods for hyperbolic conservation laws in a series of papers [5,6,7,8]. DG methods adopt many advantages from both finite element and finite volume methods, including hp-adaptivity flexibility, efficient parallel implementation, the ability of handling complicated boundary conditions, etc., making them a popular choice for numerical methods of conservation laws.
In [4], Cheng et al. studied DG methods for the one-dimensional (1D) deterministic two-way wave equations, and investigated a family of L 2 stable high order DG methods defined through a general form of numerical fluxes. A systematic study of stability, error estimates, and dispersion analysis was carried out. In a recent paper [25], Sun and Xing extended the analysis on optimal error estimate to multi-dimensional wave equations. For the DG methods with generalized numerical fluxes, one key idea to the optimal error estimate was by constructing a global projection on unstructured meshes, which will also be useful for the analysis in this paper. Recently, DG methods have been extended for stochastic partial differential equations. Li et al. investigated DG methods [17] for stochastic conservation laws with multiplicative noise and provided optimal error estimate for the semilinear equations.
In this paper, we present the DG methods for both one-and two-dimensional (2D) stochastic Maxwell equations with multiplicative noise. This is an extension of our previous work [24], where we studied multi-symplectic DG methods for stochastic Maxwell equations with additive noise. We showed in [24] 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 was also analyzed in [24]. Another related work was studied by Chen in [2], where a symplectic full discretization of multi-dimensional stochastic Maxwell equations with additive noise is provided. DG methods were used for spatial discretization and midpoint method was used for temporal discretization. A first order convergence in time and 1.5th order convergence in space was discussed in [2]. In this work, we plan to apply DG spatial discretization with generalized numerical fluxes to stochastic Maxwell equations with multiplicative noise. We will present the discrete energy growth property of the numerical solutions obtained from our semi-discrete DG scheme, following the energy law of the exact solutions. With the help of a global projection constructed in [25], we present optimal error estimate of the semi-discrete DG methods for stochastic Maxwell equations with multiplicative noise for both one-and twodimensional cases. In the two-dimensional case, we study both rectangular meshes and triangular meshes. Strong Taylor 2.0 temporal discretization is combined with the semi-discrete DG method to derive a fully discrete method for numerical implementation. Both one-and two-dimensional numerical results are presented to validate the theoretical analysis results.
The structure of this paper is as follows. Section 2 describes the DG method for one-dimensional stochastic Maxwell equations with multiplicative noises. Energy growth and the optimal error estimate of the proposed method are provided. Section 3 studies the DG scheme for two-dimensional stochastic Maxwell equations. Both rectangular meshes and triangular meshes are considered. The corresponding discrete energy law and optimal error estimate are also studied for both cases. The temporal discretization is briefly discussed in Section 4. Section 5 presents the numerical results to validate the theoretical results. Conclusion remarks are provided in Section 6.
2 One-dimensional stochastic Maxwell equations with multiplicative noise In this section, we consider one-dimensional stochastic Maxwell equations with multiplicative noise: where W t is a standard one-dimensional Brownian motion on a given probability space (Ω, F, P) and f, g are functions that satisfy the following Lipschitz continuous and linear growth assumptions: The stochastic Maxwell equations with multiplicative noise (2.1) satisfy the following energy law. A similar result was discussed in [10] for stochastic wave equations. (2.5) Proof. By utilizing the Itô's lemma and equations (2.1), we have where the last equality follows from the integration by parts and periodic boundary condition. Integrating over time leads to Integrating the second equation of (2.1) over t leads to By Itô isometry, we obtain For the same reason, it follows that Note that the process t 0 I (gu + f v)dxdW t is an Itô integral, thus it has zero expectation. Taking the expectation of (2.7) yields the continuous energy law (2.5).
In [3], Chen et al. established the regularity properties of the solution of stochastic Maxwell equations with multiplicative noise (1.2). Let M be the differential operator defined as It was shown in [3] that, under certain assumptions, the solutions are uniformly bounded in the following way: (2.10) The one-dimensional computational domain I is partitioned into subintervals I j = [x j−1/2 , x j+1/2 ], j = 1, 2, · · · , N . Let us denote x j = (x j−1/2 +x j+1/2 )/2 to be the center of each subinterval. Denote h j = x j+1/2 − x j−1/2 to be the mesh size of each subinterval and h = max j h j to be the maximum mesh size. We further assume that the ratio h/h j is bounded over all j during mesh refinement. The piecewise polynomial solution and test function space V k h is defined as where P k (I j ) stands for the space of polynomials of degree up to k on the cell I j . Since the function in V k h can be discontinuous at cell interface, we use w + h and w − h to represent the function limit from the right and left respectively. We denote the average and jump of the functions at the cell interfaces by Throughout this paper, we denote L 2 norm by · , and C represents a generic positive constant which is independent of the spatial and temporal step size h, ∆t, and may take different values in different cases.
The DG scheme for the one-dimensional model (2.1) takes the following formulation: for where u h , v h are the numerical fluxes defined on the cell interfaces. In this paper, we follow the study in [4,25] and consider the following generalized numerical fluxes for some α ∈ R and β 1 , β 2 ≥ 0. A subgroup of such numerical fluxes, named "αβ" fluxes, were considered in [4] for one-dimensional deterministic two-way wave equations, and optimal error estimate were investigated based on a specially constructed projection operator. DG method with more general numerical fluxes were studied in [25] for the one-and multi-dimensional deterministic wave equations. To provide an optimal error estimate, the key ingredient was to construct an optimal global projection on one-dimensional meshes or multi-dimensional unstructured meshes, which will also be utilized in this paper. Note that similar DG method was studied in [2] for the stochastic Maxwell equation with additive noise, using a specific choice of numerical fluxes (upwind flux) and slightly different treatment of the noise term, and suboptimal error estimate was investigated.
Next, we start by showing the following semi-discrete energy law satisfied by the numerical solutions of the proposed DG methods.
Theorem 2.2 (Semi-discrete energy law). Let v h and u h be the numerical solutions obtained in (2.11) and (2.12) with β 1 , β 2 ≥ 0, and let P(f ) and P(g) be the L 2 projections of f and g onto V k h , then we have Moreover, the equality holds when β 1 = β 2 = 0 in the numerical fluxes (2.13).
Proof. Taking the test functions ϕ = v h in (2.11) and ϕ = u h in (2.12), and adding the resulting two equations together, we have Note that the equality in (2.14) holds when β 1 = β 2 = 0 in the numerical fluxes (2.13).
By Itô's formula, we have It follows from (2.11) that therefore, for any continuous semimartingale Y , where {φ l j } represents the set of orthogonal Legendre basis of V k h over cell I j . It can be shown that (2.16) and (2.17) lead to Repeating the same process as in (2.18), we can obtain Similarly, we have Integrating (2.14) over t, summing up over all cells I j and utilizing the results in (2.15), (2.20)-(2.21), we have with periodic boundary conditions. Note that the process t 0 I gu h + f v h dxdW τ is an Itô integral, thus it has zero expectation. Taking expectation of equation (2.22), and applying Itô isometry onto which finishes the proof.
Next, we will provide an optimal error estimate analysis of the proposed semi-discrete DG method. We start by defining a pair of global projection operators which will be used in the error estimate analysis: on any cell I j and for any pair of smooth functions (q(x), r(x)), define P α,β 1 and P −α,β 2 as This set of global projections was also introduced by Sun and Xing in [25], and the following property on the projection error is studied in [25, Lemma 2.1] and will be useful in the error estimate analysis.
are well defined. Furthermore, let q, r ∈ H k+1 be smooth functions, then there exists some constant C such that Proof. Since the exact solutions v and u also satisfy equations (2.11)-(2.12), taking the difference of the semi-discrete method and the equations satisfied by the exact solutions yields the following error equations: Decompose the numerical error into the following two terms For the initial condition, we choose which is an Itô integral, hence we have E E j q (w) = 0. By choosing the test functions ϕ = ξ v , ϕ = ξ u in (2.28)-(2.29) and summing up the equations, we obtain The last equality follows from the definition of the special projection which leads to (for the error term and an integration by parts which leads to By Itô's lemma, we have We then have, for any test function ϕ, (2.34) Therefore, after integrating over t, we have (2.35) For any semimartingale Y , we have (2.36) Let us write ξ u = k l=0 (ξ u ) l j φ l j . After repeating the same process as in (2.18) and applying the Cauchy inequality, we obtain (2.37) Introduce the notations Recall the assumptions (2.2) and (2.3), and by Itô isometry we have and E(Qu 1 ) = Therefore, after taking the expectation of equation (2.37), we conclude and for the same reason, After summing over all cells I j , utilizing the periodic boundary conditions and integrating equations (2.31), (2.32) from 0 to t, we can take the expectation of the resulting equations and obtain Following the same derivation of Eq. (2.34), we have Based on the assumptions, we know that I j (P α,β 1 (g) − g) ϕdxdW t is a martingale. Therefore after applying the error estimate of the projection in Lemma 2.1. Similarly, we have Note that the initial condition satisfies ξ v (x, 0) = ξ u (x, 0) = 0, and we can utilize the results in (2.39)-(2.43) to obtain The optimal error estimate (2.27) follows from applying the Gronwall's inequality and the optimal projection error.
Remark 2.1. In the proof, we assumed enough regularity of the exact solutions to study the "best" spatial convergence rate of the proposed method. Such convergence rate has also been observed on some numerical examples in Section 5.

Two-dimensional stochastic Maxwell equations with multiplicative noise
In this section, we study two-dimensional stochastic Maxwell equations in the following form where u = (E, S, T ), and f, g, r are smooth functions that satisfy the following Lipschitz continuous and linear growth assumptions: We refer to the beginning of Section 2 on the discussion of regularity properties of the solutions. For this model, we have the following energy law satisfied by the exact solutions. The proof follows almost the same analysis as that of Theorem 2.1 and is skipped here.

Triangular meshes
In this subsection we will consider DG methods for stochastic Maxwell equations with multiplicative noise on triangular discretization of the domain Ω, and carry out the corresponding analysis. Firstly we rewrite the equation (3.1) in the following compact form: Let Ω = ∪ K∈T h K be a quasi-uniform triangulation of the domain Ω, and Γ be the collection of triangle faces. For a triangle K, and a face F ∈ ∂K, let n be the outer normal vector on F . Given a face F , let n F denote the unit vector across F , whose direction is not essential for unspecified K, and we further assume there is a unit vector r such that |r · n F | ≥ κ > 0 for all n F . We define the finite element DG space V k h as where P k is the space of polynomial of degree at most k. V h is used to denote (V k h ) 2 , and denote (P k ) 2 as P k . For a given face F , let K ± be the two neighboring cells. For any function w or vector valued function v, define The vector n ± is the outer normal vector on F with respect to K ± . Finally, we define the following notations for averages and jumps across a face F : Note that over all the triangles and faces, we have the following equalities The DG scheme for the two-dimensional system (3.5) is: where the numerical fluxes are chosen to be (3.11) with α = α sgn(r · n F )n F for some number α. Note that these generalized numerical fluxes can be viewed as two-dimensional extension of the one-dimensional numerical fluxes (2.13).
Next, we start by showing the following semi-discrete energy law satisfied by numerical solutions of the proposed DG methods.
Proof. By taking the test function ϕ = E h in (3.9) and ψ = U h in (3.10), and summing the resulting equations over all cells K, we obtain (3.13) By an integration by parts and applying the equality (3.8), we can have (3.14) Therefore, equation (3.13) becomes By Itô's lemma, we have Following an exact same process as in the derivation of (2.20), and applying the Itô isometry, we have By plugging these into (3.15), integrating over time t, summing over all K, and taking expectation, we obtain where we use the fact that the right-hand side of (3.13) Ω F E h +G·U h dxdW t is a martingale, hence its expectation equals to zero. It is easy to observe that the equality holds when β 1 = β 2 = 0.
Similar to the one-dimensional case, to provide the optimal error estimate, we need to introduce the following pair of projections P U U and P E E [25]: for any K ∈ T h , (3.20) The following lemma on the projection is studied in [25, Lemma 3.1] and will be useful in the analysis of error estimate: Lemma 3.1. Suppose β 1 ≥ 0, β 2 > 0, and |α| 2 + β 1 β 2 = 0, then the projection pair defined in (3.17) -(3.20) is well defined, and there exists some constant C independent of mesh size h, such that Proof. For simplicity, for f, ϕ ∈ V k h , and g, ψ ∈ V k h , introduce the notation We further define and choose the initial condition It can be observed that the following error equations hold where we dropped the dependence of F , G on x, t for simplicity. Taking the test functions ϕ = ξ E and ψ = ξ U , and summing up two equations, we obtain From our definition of projections (3.17)-(3.20), we can conclude that A K ( E , U ; ξ E , ξ U ) = 0. Following the exact same analysis as in (3.14), we can have In addition, the term M 1 : · ξ U dxdW t on the right-hand side is an martingale, and its expectation equals to zero. Therefore, summing the equation (3.21) over all the cells K leads to Following the same derivation of Eqs. (2.42)-(2.43), we have By Itô lemma, we have Follow the same process as the steps in (2.33)-(2.39) in one-dimensional case, we have

(3.24)
Summing over all the cells K, combining these results with (3.22), integrating over time t, and taking expectation, we will have after utilizing the optimal projection error in Lemma 3.1. Note that the chosen initial condition satisfies ξ E (x, 0) = ξ U (x, 0) = 0, we have after applying the Gronwall's inequality and the optimal projection error.

Rectangular meshes
In this subsection, we will investigate DG methods on cartesian meshes, and study the stability and error estimate of the proposed methods. The stability result and its proof is similar to the case of triangular meshes, but different technique is needed for the proof of optimal error estimate, which will be discussed in details.
The two-dimensional rectangular computational domain Ω is set to be I × J, and we consider the rectangular partition with the cells denoted by ] for i = 1, 2, · · · , N x and j = 1, 2 · · · , N y . Let , and y j = 1 2 (y j− 1 2 + y j+ 1 2 ) . Furthermore, we define the mesh size in both directions as h , 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 V k h as follows: , i = 1, 2, · · · , N x ; j = 1, 2, · · · , N y }.
The DG scheme for (3.1) is as follows: where, for simplicity, the following operators are defined: for α ∈ R, p, q ∈ V k h , with the numerical fluxes defined as follows: Note that the numerical fluxes (3.28) can be viewed as a special case of (2.13) with β 1 = β 2 = 0.
They are chosen such that the optimal error estimate can be easily analyzed.
We first provide the following semi-discrete energy law satisfied by numerical solutions of the proposed DG methods on rectangular meshes. The proof is identical to that of Theorem 3.2, and is skipped here.
Before we study the error estimate, some preparations on projections are provided. Let us define the generalized Radau projection as where P is the standard one-dimensional L 2 projection, and P α,0 x , P β,0 y are the one-dimensional generalized Radau projections (2.23)-(2.26) defined as follows: On the cell I i , J j and for any function q(x), r(y), the projections P α,0 x , P β,0 y into the space V k h are given by (P β,0 y r − r(y))ψ(y)dy = 0, ∀ψ(y) ∈ P k−1 (J j ), and ({P β,0 y r} + α[P β,0 r]) j+ 1 2 = r(y j+ 1 2 ).
The following lemmas are studied in [20] and will be useful in our error estimate analysis.
Lemma 3.2 (Superconvergence property). Let P α,β be the projection defined in (3.30) with α, β = 0. For any function w(x, y) ∈ H k+1 , denote = P α,β w − w. For any φ ∈ V k h , there exists some constant C such that Now we turn to the optimal error estimate.
Proof. We start by introducing which leads to We choose the initial conditions as follows: After taking expectation on equation (3.35), summing over all the cells, and integrating over t, we have Following the same derivation of Eqs. (2.42)-(2.43), we have after applying the superconvergence property in Lemma 3.2. We utilize the projection error property and Young's inequality to obtain The optimal error estimate (3.31) follows from applying the Gronwall's inequality and the optimal projection error.

Temporal discretization
After DG spatial discretization, the semi-discrete DG methods (2.11)-(2.12), or (3.9)-(3.10) can be obtained. To solve the resulting stochastic differential equations, we present Taylor 2.0 strong scheme [17] as the temporal discretization in this section.
Let us consider the general matrix-valued stochastic differential equations of the form with i = 0, 1, · · · , m and j = 0, 1, · · · , M + 1. Suppose that Y i,j n is a numerical approximation of X i,j tn at the time t n , and Y i,j 0 = x i,j 0 . The recurrent equation between Y i,j n and Y i,j n+1 is derived below. Define the following operators: and where f : R (m+1)×(M +2) → R is a twice differentiable function. As studied in [15,Section 10.5] and [17, Appendix A], the following Taylor scheme has 2.0 order of convergence: where   at the moment s = 1. The system (4.3) can be solved numerically to obtain the approximation of ∆W , ∆Z and ∆U . We refer to [17,Appendix A.2.2] for the details of implementation.

Numerical Experiment
In this section we present numerical results of the proposed scheme for the one-dimensional and two-dimensional stochastic Maxwell equations with multiplicative noise. We use polynomials of degree k in our proposed DG methods for spatial discretization, and Taylor 2.0 strong scheme for the temporal discretization. The accuracy tests are provided for both 1D case and 2D case to demonstrate the convergence rate. The time history of the discrete energy is also presented for both examples. The Monte-Carlo method is used to compute the stochastic term, and the expectation is computed by averaging over all the samples.
We first consider the following one-dimensional equations with periodic boundary conditions and the exact solutions given by v = (sin(x − t) + cos(x + t))e Wt− 1 2 t , u = (sin(x − t) − cos(x + t))e Wt− 1 2 t .

(5.1)
The computational domain is set to be [0, 2π], and the final time T = 0.5. Nx and Nt are used to denote the number of space steps and time steps, and we use uniform meshes as spatial discretization. The numerical initial condition is taken as the projection of the exact solutions (5.1) at t = 0. We apply Monte Carlo simulation with 3000 samples to approximate the expectation. Table 5.1  and Table 5.2 show the convergence rate when k = 1 and k = 2 respectively. In both examples ∆t is chosen small enough to ensure the spatial error dominates, and we can observe the optimal error estimate (the expected (k + 1)-th order of convergence), which is consistent with the result in Theorem 2.3 for the semi-discrete DG method. Next we consider the two-dimensional stochastic Maxwell equations with periodic boundary conditions. The exact solutions take the form      E = sin(x + t) − cos(y + t) e Wt− 1 2 t , S = sin(x + t)e Wt− 1 2 t , T = cos(y + t)e Wt− 1 2 t .

(5.2)
The space domain is [0, 2π] 2 . The numerical initial condition is taken as the projection of the exact solutions (5.2) at t = 0. Nx, Ny and Nt are used to denote the number of space cells in x, y directions, and the number of time steps, and uniform rectangular meshes are considered as spatial discretization. We run the simulation until final time T = 0.1. Monte Carlo simulation with 500 samples are used to approximate the expectation. Table 5.3 and Table 5.4 present the convergence rate for the case of k = 1 and k = 2 respectively. We can observe that in both cases, the optimal convergence rate is achieved. The discrete energy law satisfied by the numerical solutions was studied in Theorem 2.2 for one-dimensional system, and in Theorems 3.2 and 3.4 for two-dimensional system. In Figure 5.1, the time history of averaged energy is shown for two cases.

Conclusion Remarks
In this paper we applied high order DG scheme for one-and two-dimensional stochastic Maxwell equations with multiplicative noise. We provide the semi-discrete energy law for both cases. Optimal error estimate of the semi-discrete method is obtained for one-dimensional case, and twodimensional case on both rectangular meshes and triangular meshes under certain mesh assumptions. The semi-discrete method is combined with strong Taylor 2.0 temporal discretization, and numerical results are presented to validate the optimal error estimates and the growth of energy.