{'fraction_non_alphanumeric': 0.032720063735571486, 'fraction_numerical': 0.01474811479813981, 'mean_word_length': 4.549455068887518, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A multiscale continuous Galerkin method for stochastic simulation and robust design of photonic crystals

We present a multiscale continuous Galerkin (MSCG) method for the fast and accurate stochastic simulation and optimization of time-harmonic wave propagation through photonic crystals. The MSCG method exploits repeated patterns in the geometry to drastically decrease computational cost and incorporates the following ingredients: (1) a reference domain formulation that allows us to treat geometric variability resulting from manufacturing uncertainties; (2) a reduced basis approximation to solve the parametrized local subproblems; (3) a gradient computation of the objective function; and (4) a model and variance reduction technique that enables the accelerated computation of statistical outputs by exploiting the statistical correlation between the MSCG solution and the reduced basis approximation. The proposed method is thus well suited for both deterministic and stochastic simulations, as well as robust design of photonic crystals. We provide convergence and cost analysis of the MSCG method, as well as a simulation results for a waveguide T-splitter and a Z-bend to illustrate its advantages for stochastic simulation and robust design.


Introduction
A field that has attracted significant interest in recent years is the study of light propagation through photonic crystals, which are heterogeneous materials engineered to exhibit properties that cannot be found in homogeneous materials. Photonic crystals are assembled by combining conventional materials in lattice structures usually at the microscopic level. Electromagnetic wave propagation through photonic crystals is characterized by band gaps, that is, ranges of frequencies for which light waves are not allowed to propagate through the periodic optical nanostructure [49]. Introducing defects in these periodic structures allows us to create waveguides, since waves traveling at frequencies within the band gap are exponentially attenuated within the crystal, and thus can only propagate along the defect. Waveguides exploiting this mechanism can be much more efficient than traditional waveguides based on total internal reflection (TIR). Photonic crystals have applications in fibers [34,46], waveguides [29,32] and superlenses [37,41].
The simulation of wave propagation phenomena in heterogeneous media is an active research field and several techniques have been proposed. One of the most widely used approaches is the finite-difference time-domain (FDTD) method [35,40,52,53]. The FDTD method is simple and efficient, albeit not well suited to treat complex geometries, irregular domains and multiple length scales. Another family of approaches are the finite-volume time-domain (FVTD) methods [3,4], which enable mesh adaptation and refinement and exhibit low dissipation and dispersion due to the use of high-order Godunov schemes. The finite element (FE) method both in the time and frequency domain [28,48,57] is also a popular alternative for solving wave propagation problems thanks to its ability to handle complex geometries and inhomogeneous materials, as well as allow h/p adaptivity. Furthermore, material interface conditions and boundary conditions can be implemented in a natural manner. Although low order methods are often used due to their simplicity, high order methods [2,20] are more accurate and efficient if high accuracy is required. All the above methods have been used to study the propagation of waves in photonic crystals [18,21,27,36,44,50].
The ability to accurately simulate wave propagation in photonic crystals presents unique challenges. Problems of interest typically involve complex geometries and a mismatch in critical length scales, which can be of several orders of magnitude. Resolving the small scales using uniform discretizations requires a prohibitive large number of grid points. An additional difficulty stems from the fact that mathematical models (e.g., based on Helmholtz and Maxwell's equations) may not capture the real physical phenomena accurately enough due to simplifications and uncertainty in the model data, such as geometry, material properties, and boundary conditions. Several methods have been developed to deal with multiple scales in composite materials [33]. Homogenization methods [6] allow the treatment of multiscale features by solving a coarse-scale model which has been modified to account for the small scales. These approaches can be quite successful at predicting the global macroscopic behavior but lack the detailed description of the physics at the smaller scales. The multiscale FE method [23,24] and the mixed multiscale FE method [14] are also alternatives that have been successfully applied to multiscale elliptic problems. The main idea is to construct multiscale FE basis functions on a coarse-scale that capture the fine-scale features of the solution. These multiscale functions are then coupled globally into a linear system, whose unknowns are the solution values on the coarse grid. The finescale solution can be recovered combining the multiscale functions with the coarse-grid nodal values. The main drawback of these methods is their strong dependence on the boundary conditions of the subproblems and, for problems where strong heterogeneities are present at the interfaces, one needs to develop adaptive boundary conditions to avoid small scale resonances [24].
In our approach, we advocate for not modifying the original multiscale problem, but instead obtaining computational efficiency via domain decomposition. A successful domain decomposition approach is the mortar element method [9,10]. This method allows for an independent discretization of each subdomain and enforces continuity of the solution across subdomains weakly. Since the meshes across the subdomain boundaries do not need to match, this approach allows great flexibility in the definition of the subdomains, and becomes particularly attractive for problems requiring mesh adaptivity and complex geometries. Other methods include the multiscale DG method [1], which blends the multiscale FE method and imposes weak continuity at the subdomain interfaces, the hybridized multiscale DG method [39], the geometric multiscale FEM [11], and the method of polarized traces [58] used for high-frequency problems. For an extensive review of multiscale FE methods, we refer the reader to [16].
In this paper, we present the multiscale continuous Galerkin (MSCG) method to simulate wave propagation problems on structured materials. This method is an extension of the hybridized continuous Galerkin method (CG) introduced in [15] and the hybridized multiscale DG method [39]. The multiscale CG method possesses considerable advantages over other simulation methods. First, the multiscale discretization results in a significantly smaller global linear system than that of the original problem, owing to static condensation of all the degrees of freedom corresponding to the domain interiors. Second, it exploits repetitive patterns in the structure to rapidly construct the global linear system by only solving a small number of local subproblems. As a result, the method provides fast and accurate simulations of wave propagation in structured media beyond the capabilities of current numerical methods.
Furthermore, we also consider the simulation and optimization of wave propagation problems in the presence of fabrication uncertainties. For complicated structures the interaction of electromagnetic waves with het-erogeneous materials can be very sensitive to geometry, thus small perturbations in the microscopic structure may significantly degrade the performance of the device. To that end, we augment the MSCG method with the following ingredients: (1) a reference domain formulation [42] that allows us to treat geometric variability; (2) a reduced basis approximation [5,43,45] at the subdomain level [25,26] that drastically accelerates the solution of the parametrized local subproblems; (3) an adjoint technique for computing gradients of an output functional; and (4) a model and variance reduction technique [54,55] that enables the accelerated computation of statistical outputs by exploiting the statistical correlation between the MSCG solution and the reduced basis approximation.
This article is organized as follows. In Section 2, we introduce the wave propagation problem and present the MSCG method. In Section 3, we extend our approach to stochastic simulation and robust design problems. In Section 4, we present numerical results to demonstrate the performance of the proposed method. Finally, we summarize the main conclusions of this article in Section 5.

Problem statement
The propagation of the time-harmonic electric and magnetic fields E, H is governed by the frequency-domain Maxwell's equations. However, for two-dimensional problems the linearity of Maxwell's equations implies the transverse magnetic (TM) and transverse electric (TE) polarizations decouple, hence it is only necessary to solve the scalar Helmholtz equation in a domain Ω ∈ R 2 , with Lipschitz boundary ∂Ω = ∂Ω D ∪ ∂Ω N where Dirichlet and Neumann conditions are both prescribed, namely: The auxiliary parameters ρ and κ enable the compact definition of the above non-dimensional equation for both polarizations. In TM or E-polarization, the magnetic field is confined to the plane of propagation (which we assume the x − y plane), that is H = (H x , H y , 0), and the electric field is perpendicular to this plane E = (0, 0, E z ), hence ρ = 1, κ 2 = (ω/c) 2 ε(x) and u = E z . Conversely, in TE or H-polarization, the electric field is confined to the plane and the magnetic field is perpendicular to it, that is E = (E x , E y , 0) and H = (0, 0, H z ), hence ρ = ε(x) −1 , κ 2 = (ω/c) 2 and u = H z . The dimensionless relative permittivity field ε(x) is a piecewise-constant function that prescribes distinct permittivity values to different spatial regions within the photonic crystal, thus enabling the simulation of heterogeneous materials. The angular frequency of the wave ω is non-dimensionalized with the speed of light c and the periodicity of the photonic crystal a as ω = ωa/(2πc).
The Helmholtz equation above assumes a finite computational domain. Nonetheless, wave propagation problems often occur in unbounded domains. In order to simulate unboundedness while restricting our simulation to a finite domain, we resort to Perfectly Matched Layers (PMLs) [7] surrounding the photonic crystal. PMLs can be encoded in the auxiliary parameters with the transformations where s x is a complex-valued frequency-dependent damping function that prescribes the attenuation of the x-propagating waves inside the PML, and analogously for s y . We refer to [30] for a detailed description of PMLs, their numerical implementation and definition of damping functions.

Structured heterogeneous materials
As mentioned earlier, our goal in this paper is to develop a method to efficiently simulate wave propagation through heterogeneous structured materials. Heterogeneity is encoded in ε(x) by assigning different permittivity values to distinct spatial regions forming the photonic crystal. As a motivating example, we consider the waveguide shown in Fig. 1 (left) consisting of a regular arrangement of dielectric rods (ε(x) > 1) in air (ε(x) = 1), and a row defect without rods. Our work relies on decomposing the domain into smaller subdomains that belong to distinct classes, hence within each class, the geometry of the subdomains is the same. In this case, only two different classes of subdomains are needed: the subdomains with the rod and subdomains without the rod. Unfortunately, in real photonic applications the rods are not identical, mainly due to manufacturing variability that arise as a consequence of the extreme-scale fabrication techniques as sketched in Fig. 1 (right). Thus, in order to achieve truly predictive simulations one must take into account these variations. In this paper we consider only circular dielectric rods, but the extension to other types/shapes of rods is straightforward. Figure 1: Reference domain (left), physical domain (right) and subdomain mapping.

Reference domain formulation
We assume that the geometry of the physical domain Ω is parametrized by z, defined on a compact set, and that we want to solve (1) for many different realizations of z. In this scenario, it is much more convenient to map the physical domain Ω onto a fixed reference domain Ω r . For the waveguide example in Fig. 1, the reference domain Ω r is the structure with perfectly circular rods in Fig. 1 (left), whereas the physical domain Ω is the real waveguide with irregular rods in Fig. 1 (right). Furthermore, we can leverage the structure of the photonic crystal partitioning both the physical and the reference domain into M nonoverlapping subdomains, that is Ω = Solving (1) for a given z first requires either re-discretizing the transformed geometry or deforming the original discretization, which can be complex and cumbersome. Alternatively, we follow Persson et al. [42] and use the diffeomorphism G that maps x r ∈ Ω r to x ∈ Ω to transform the governing equations (1) into a new set of equations defined on the reference domain Ω r that remains unchanged. To derive these equations, we first integrate on a volume and use the divergence theorem. The mapping G has deformation gradient G = ∇ r G and Jacobian g = det G, and for an arbitrary control volume v ⊂ Ω (with Lipschitz boundary) and its associated v r ⊂ Ω r , the transformed equation on the reference domain reads where we have invoked the relationships ∇ = G −T ∇ r , dv = gdv r and nds = gG −1 n r ds r with outwardpointing normals n, n r . Transforming the boundary conditions, we arrive at where G = gG −1 G −T and g s is the restriction of the Jacobian g to the face. Thus, for a given z instead of solving (1) on the deformed domain Ω, we may alternatively solve (3) on the original domain Ω r .

MSCG method
Let us introduce the MSCG method to solve the transformed governing equations (3) on the fixed reference domain Ω r . We first define a discretization T m h on each subdomain Ω m r , consisting of disjoint regular elements T that partition the subdomain. The discretization of the entire reference domain is thus T h = M m=1 T m h . Subdomain discretizations are disjoint between subdomains, thus one can seek for instance a triangulation for one subdomain and a quad discretization for its neighboring subdomain. The division of a sample reference domain into five subdomains with its corresponding discretizations is illustrated in Fig. 2a.  Let f be a subdomain face given by either ∂Ω m r ∩ ∂Ω n r (n = m) or ∂Ω m r ∩ ∂Ω r . The collection of subdomain faces of T h is {f } L =1 . Furthermore, each subdomain face f may be subdivided into N elements f i , 1 ≤ i ≤ 5 N . Thus we define the set of global face elements as F f = {f i , 1 ≤ i ≤ N , 1 ≤ ≤ L}, also referred to as the mesh skeleton. We note that, not only the discretizations of the faces and the boundaries of the adjacent subdomains are not required to coincide, but there is also flexibility in the way we define the subdomain faces. For instance, in Fig. 2a, the subdomain faces coincide with the global face elements except for the interface between Ω 5 r and both Ω 3 r , Ω 4 r , wherein the subdomain face f 1 is subdivided into two global face elements f 1 1 (gold) and f 1 2 (blue). We can now introduce our approximation spaces as where P p m (T ) is the space of polynomials of degree at most p m on T ∈ T m h and P p f (f) is the space of polynomials of degree at most p f on f ∈ F f . Note that we allow for polynomial spaces of different degrees for both different subdomains and face elements on the subdomain's boundaries. To impose Dirichlet boundary conditions, we introduce the set where the operator P V h represents the L 2 projection onto the space V h on the boundary ∂Ω r,D . The volume inner products for these finite element spaces are defined as and the surface inner products as We then introduce the auxiliary variable q h which approximates the normal component of the flux q = ρG∇ r u · n r , and seek an approximation The next step is to eliminate the unknowns of (u h , q h ) to obtain a formulation only in terms of the variables λ h defined on the subdomain interfaces. Invoking linearity and superposition, we consider two subproblems for each subdomain Ω m r : the first subproblem maps the function where the operator P W m h represents the projection onto the restriction of W m h on the subdomain boundary. Finally, we find the Lagrange multiplier λ h ∈ V h (u D ) as a unique solution of the weak formulation One of the main features of the MSCG method is the flexibility of mesh generation and h/p adaptivity at the subdomain level, since the local subproblems are independent and are only coupled through the interface. Hence, the local subproblems may be statically condensed and expressed in terms of the variables defined on the subdomain faces. The idea is to relax the solution continuity at the subdomain interface and impose it back through a set of Lagrange multipliers λ, which correspond to the unique solution of the variational formulation on the mesh skeleton. Once the global solution is obtained, the solution for the local subproblems can be recovered independently for each subdomain. Another important advantage of the MSCG method is that the solution at the interfaces can be approximated with a different polynomial space than that of the subproblems, for instance using high-order Lagrange interpolation polynomials defined on the Chebyshev nodes. This approach results in a drastic reduction of the size of the global system that would otherwise need to be solved had (u h , q h ) not been eliminated. Finally, continuity of the approximate solution across subdomains is enforced by imposing Dirichlet boundary conditions for each subdomain determined by the value of the global solution at the interfaces.

Implementation
Let the total number of global (resp. local) degrees of freedom in the problem be given by N G (resp. N L ), and let N m G (resp. N m L ) denote the number of global (resp. local ) degrees of freedom associated to the m-th subdomain. Let the space V h be spanned by a set of global basis functions For each global face element f ∈ F f , these global basis functions are chosen to be nodal Lagrange polynomials with the nodes placed at the Chebyshev points, as shown schematically in Fig. 2b. where The formation of the global stiffness matrix K ∈ C N G ×N G and load vector F ∈ C N G is detailed below. Firstly, at the subdomain level we obtain u f m h and u is the set of global basis functions that have non-zero support on ∂Ω m r . A total of N m G + 1 linear systems of dimension N m L × N m L are therefore solved for the m-th subdomain. The corresponding local stiffness matrix and load vector are computed as We then assemble K m , F m for m = 1, . . . , M on the global stiffness matrix K and load vector F through the standard finite element assembly procedure. After computing the global unknowns Λ from the linear system (4), the solution on each subdomain Ω m r may be recovered as All in all, in the absence of geometry deformations (that is, G = 1), only the local problem (5)-(6) for one representative of each subdomain class needs to be solved, since all the subdomains within the same class are identical. This feature is one of the key advantages of the MSCG method, which makes it particularly attractive for problems with periodic or repetitive structures. Indeed, a judicious subdomain partition enables the simulation of complex photonic crystal structures just by solving a handful of small local subproblems and a single global linear system (4) defined only on the skeleton of the mesh.
However, for real applications the geometry of each subdomain within a class may not be exactly the same, see Fig. 1, hence it would appear that this advantage is lost. In the next section, we introduce a reduced basis formulation that allows us to consider a more general parametrized problem that is, in fact, able to account for the small variations in the geometry that occur from subdomain to subdomain, and thus to partially retain the computational advantage of the MSCG method.

3-D extension
The formulation and implementation described above for the Helmholtz equation can be naturally extended to three-dimensional structures. However, there are some aspects of going from 2-D to 3-D that should be discussed.
The main obstacle is the increased computational cost that stems from considering linear systems that are larger and denser, that is going from 2-D to 3-D in solving the local subproblems (5) and from 1-D to 2-D in solving (4) on the global face elements. In addition, since the Lagrange polynomials that approximate the global solution are now defined on surfaces, for each subdomain we have more global face elements and Dirichlet conditions per face. More precisely, using a tensorial grid on each global face squares the number of Dirichlet conditions. Even though sparse interpolation grids [51] can help mitigate this issue, the cost of (5b) augments significantly.
Finally, if the simulation of full three-dimensional photonic crystals is sought we can no longer reduce Maxwell's equations to Helmholtz, since the polarizations are coupled. Changing from Helmholtz to Maxwell is not straightforward from both the formulation and computation standpoints, since we now have to solve for six coupled three-dimensional fields instead of one, while ensuring the formulation is curl-conforming. To that end, a multiscale method for Maxwell's equations has been proposed and applied to 3-D optical fibers and 2.5-D photonic crystals in the dissertation [47], where the reader is referred for further details.

Stochastic MSCG
In this section, we introduce an approach for computing the statistics of quantities of interest, such as the transmission power of photonic devices. The main ingredients of our approach include: stochastic modeling to deal with geometry variability, a reduced order model [5,43,45] to economize the simulation at the subdomain level, a variance reduction method [19,54,55] to accelerate the convergence of the statistical outputs, and an adjoint technique to compute derivatives.
In photonic crystals, the most widely used patterns are either dielectric circular rods in air or air holes in a dielectric slab [29]. Ideally, all of these circular structures would be circles and have the same radius R 0 . However, because of technological limitations at micro and nano scales, they are never perfectly round. For this work, we focus on photonic crystals with circular shapes, and model only geometry variations in the radial direction, albeit more geometries and complex deformations can be easily accommodated within our framework. The assumptions and mapping used for the examples in this article are thoroughly described in Appendix A.

Reduced basis approximation
A key component of our proposed approach is the use of a reduced basis (RB) approximation to efficiently compute the subdomain solutions for different rod geometries. For electromagnetism there has been considerable work to incorporate geometry as a parameter of the RB model [13,21,22,56]. Indeed, if the high-fidelity solutions lie in a low-dimensional manifold induced by the parametric dependence, then they can be accurately approximated using RB solutions at only a fraction of the cost, due to its rapid convergence. However, incorporating a reduced basis approximation in our multiscale framework is not straightforward, thus its construction is outlined below.
The bilinear forms involved in the local problems (5)-(6) depend non-affinely on the parameter vector z since the deformation mapping G is non-affine. To circumvent this issue, we use the empirical interpolation method (EIM) [5] at the subdomain level. For the m-th subdomain we approximate G m and g m with the following affine expansions In our context, we compute for multiple values of z m the fields G m and g m evaluated on the Gaussian quadrature points of the subdomain discretization, required to evaluate the volume integrals of the weak formulations (5)- (6). Proper orthogonal decomposition (POD) [8] is then applied to attain the above expansions, keeping the modes necessary to retain at least (1 − 10 −8 )% of the system's energy, thus guaranteeing the approximation error is not dominated by the empirical interpolation. Naturally, the affine expansions may be reused for subdomains in the same class.
We then substitute (7) into (5) and seek u f m h ∈ W m h (0) and u for all w ∈ W m h (0). For subproblem (8b) assume that we are given orthonormalized basis functions ζ i n , 1 ≤ n ≤ N max , such that (ζ i n , ζ i n ) W m h = Ω m r ζ i n ζ i n = δ n,n . The associated hierarchical RB space is defined as W m i,N = span{ζ i n , 1 ≤ n ≤ N }. In our case this spaces are constructed by computing a POD on a set of solutions. The main caveat of this formulation is that we typically have many Dirichlet conditions, which can impact the efficiency of the method since we must solve multiple RB systems. Instead, we propose to develop a single RB space W m N , where the basis functions ζ n are computed with a weighted POD on a set of solutions for all possible Dirichlet boundary conditions. In other words, we treat the Dirichlet condition as an additional parameter, hence include the solutions for the various boundary conditions as independent snapshots. This strategy allows us to assemble and solve only one reduced basis system per subdomain group. For the first subproblem (8a), a similar process is followed if f has spatial dependence. In such case, the RB needs to be constructed for all possible variations of the source on the subdomains. Nonetheless, if the source is constant, the RB for (8a) may be obtained using the standard RB procedure, and we denote the associated RB space as W m N (0).
We then apply a Galerkin projection to find N m The local stiffness matrix K m ij and load vector F m i are then approximated as The affine parametric dependence enables an offline-online computational strategy where the parameter-independent instances are precomputed and stored beforehand, and the cost to obtain the approximate u  [38] can be leveraged to regain well-posedness, at the expense of a higher computational cost. However, in the numerical experiments presented here the natural Galerkin projection described above the RB is well-posed, mainly due to the low frequencies of interest in photonic crystal design, thus the minimum-residuals formulation has not been used for this work.

Model and variance reduction
Our approach exploits the structure of the problem using the MSCG method and the RB for subdomains under a deformation mapping. We now study how these two ingredients can be combined to produce fast yet accurate estimates of the statistics of a quantity of interest. We propose a multilevel variance reduction method that exploits the statistical correlation between the different reduced basis approximations and the high-fidelity MSCG discretization to accelerate the convergence rate of Monte Carlo simulations. The multilevel variance reduction method provides efficient computation of the statistical outputs by shifting most of the computational burden from the high-fidelity MSCG approximation to the reduced basis approximations. The methodology presented below is the direct extension of our previous work [54,55] for the multiscale setting.
We first introduce a probability space (Υ, Σ, P ), where Υ is the set of outcomes, Σ is the σ-algebra of the subsets of Υ, and P is the probability measure. If Z is a real random variable in (Υ, Σ, P ) and υ a probability event, we denote its expectation by E[Z] = Υ Z(υ)dP (υ). For an arbitrary subdomain Ω m r , we shall consider functions w ∈ L 2 (Ω m r × Υ) equipped with the following norm We assume that, for a given subdomain Ω m r , the parameters that define the geometric mapping z m d (υ) for d = 1 . . . , D are mutually independent random variables with zero mean. In addition, we assume that each of the z The solution u of the original problem can be written as a function of z ∈ Γ. Now let s be a bounded functional. We introduce a random output s defined as s(z) = s(u(·, z)).
We are interested in evaluating the expectation and variance of s as Since the exact output cannot be computed, we introduce the MSCG and RB outputs defined as s h (z), s N (z) respectively. The high-fidelity output s h (z), computed by the MSCG method, is a very accurate approximation to s(z) for any z ∈ Γ. Conversely, s N (z) denotes the output computed when we use an N -dimensional RB to evaluate (9)-(10) on the subdomains with parametric dependence, as well as the high-fidelity CG formulation (5)-(6) for the remaining subdomains. Moreover, since the RB spaces provide a rapidly convergent approximation to the CG solution, we expect a high statistical correlation between both outputs.
We now apply the above idea to compute an estimate of E[s h ]. To achieve this goal, we introduce where s N (z) is the RB output for some N ∈ [1, N max ]. Because s N (z) generally approximates s h (z) very well, the two outputs are highly correlated. The expectation may be recast as The underlying premise here is that the two expectation terms on the right hand side can be efficiently computed owing to variance reduction and model reduction: the first term requires a small number of samples because its variance is generally very small, while the second term is less expensive to evaluate because it only involves the RB output. The model and variance reduction (MVR) estimate of the expectation is computed through Monte Carlo simulations of the terms in (11) with (I 0 , I 1 ) i.i.d. samples, that is which is an unbiased estimator of the expectation.
The application of the CLT enables us to derive a posteriori estimate for the error in the expectation, that is where β > 0 is the confidence level. For instance, in order to guarantee that E[s h ] − E I0,I1 [s h ] is bounded by a specified error tolerance tol with a high probability (say, greater than 0.95), we need to take β ≥ 1.96 according to the CLT. Similarly, the estimate of the variance is defined as where the auxiliary variables are The error bound for the variance is defined analogously. The above derivations, the extension to the multilevel context and the optimal choice of (I 0 , I 1 ) and N are thoroughly discussed and analyzed in [55].

Computation of gradients
In the optimization context, the usage of first order optimization algorithms usually leads to accelerated convergence to the optimum value, as it guarantees a more efficient exploration of the design space. We now review how to obtain the gradients for the MSCG method.
Following our assumptions, the long vector of parameters z has at most dimension M D, since the subdomains not containing rods do not require any parameters. We now derive an adjoint approach to compute gradients in the multiscale context. The solution to the MSCG problem in matrix form is u = U ϕ Λ + U f , where U ϕ ∈ C N L ×N G is a matrix that contains the solutions to the Dirichlet subproblems and U f is a vector of dimension N L that contains the solution to the source problem. Note that both U ϕ , U f are complex-valued, since the governing equations on PML subdomains are characterized by imaginary values of ρ and κ 2 per (2).
The derivatives of s(u, z) can be recovered as We invoke the adjoint technique to solve the last part, namely where † is the conjugate transpose operator, and the adjoint variable ψ is the solution of thus allowing us to compute the (at most) M D derivatives with only an additional adjoint solution.
The derivatives of U ϕ , U f , K, F with respect to z are sparse since the parameters only have local influence, and can be evaluated with (5)-(6) for the MSCG case and with (9)-(10) for the RB case. An additional advantage of RB is the computation of these derivatives, since the affine dependence on the parameters that arises from the EIM greatly simplifies the computation.

Numerical results
The MSCG algorithm introduced above is implemented in Matlab. To maximize computational efficiency, the routines to calculate the local stiffness matrix and load vector for CG (6) and RB (10) are fully vectorized. In order to present fair time comparisons, the same script is used to solve the entire multiscale problem, and the only difference is the routine used for local subproblems, either CG or RB. Thus, the global problem is assembled and solved identically for both CG and RB. Time estimates are obtained averaging the wall time for 500 runs for each task using a single processor of a 512GB Linux 18.04 machine with 32 AMD Opteron(tm) Processors 6320x15. For the results presented here, the multiple different local subproblems are solved in series, and parallelization is applied to the MVR computation of statistical moments given by (12) and (15).

Convergence test and cost analysis
We start by studying the convergence of the MSCG method and compare it to that of the standard CG method for a simple problem. We consider problem (1) in a unit square Ω = (0, 1) 2 for ρ = κ 2 = 1, where the Dirichlet condition and source term are chosen such that the exact solution is u(x) = x 2 + y 2 + sin(k(x cos θ + y sin θ)).
The exact solution is a plane wave of wavenumber k propagating in the θ direction. The results reported here correspond to θ = π/4 and k = 6.
For the CG discretization, we consider a homogeneous discretization of 2n 2 triangular elements with uniform element size of h = 1/n. For the MSCG, we subdivide the domain into q 2 uniform squares, each of them with 2n 2 /q 2 triangles of uniform element size of h = q/n. To represent the solution, we consider both polynomial order p = 1, 2 for all subdomains; for the interfaces, besides considering a maximum polynomial order sufficient to capture the frequency of the problem, for instance 10, we also require for any boundary face that the local discretization is finer than the global discretization in order to capture the boundary conditions prescribed by the Lagrange polynomial. Hence we take p f = min{10, np/q}, and discretize each global face with a single high-order element.
We compute the L 2 (Ω) errors for both methods for n = 8, 16, 32, 64, 128. Since the global polynomial order suffices for the frequency of interest, the error is dominated by the subproblems, where both methods have the same discretization, hence both the CG and MSCG errors for different q values coincide, and yield the expected convergence rate of O(h p+1 ) for smooth solutions, see Table 1. Another relevant study is the comparison of degrees of freedom for both methods. A uniform triangular mesh of 2n 2 elements of order p renders N CG = (np + 1) 2 high-order nodes. Conversely, for each of the q 2 subdomains of the MSCG we have 2(n/q) 2 elements and N L = (np/q +1) 2 degrees of freedom. Furthermore, if we assume a homogeneous polynomial order p f for the Lagrange polynomial at the interfaces, it can be shown that N G = (q + 1)(2qp f − q + 1). If we adopt p f = min{10, np/q}, which suffices to represent the solution accurately for this case, we obtain the degrees of freedom presented in Table 2 for p = 2. 8  81  93  25  145  9  225  289  16  289  117  81  305  25  513  1089  32  1089  117  289  385  81  1089  4225  64  4225  117  1089  385  289  1377 16641  128 16641  117  4225  385 1089 1377 66049 This analysis shows that, for this problem where there is only one subdomain type, the best strategy is to use more subdivisions per direction as we refine the mesh. Indeed, a multiscale configuration is deemed optimal whenever the number degrees of freedom for both the local and global problems are similar, thus the computational burden is divided evenly. Furthermore, notice that the benefit of the MSCG as opposed to regular CG becomes apparent as more elements are used in the discretization. The results presented here correspond to the very simple case where only one local computation needs to be performed, therefore if there is more than one subdomain type one needs to account for the multiple unique local subproblems.

Waveguide T-splitter simulation
We now focus on the simulation of waveguides in photonic crystals. Photonic crystals are assembled by combining different materials, giving rise to periodic nanostructures that exhibit the bandgap phenomenon. That is, there exists broad bands of frequencies for which wave propagation through the crystal is disallowed.
Common examples of such structures are obtained by uniformly placing dielectric rods in air or drilling regular patterns of holes in a dielectric slab. The bandgap response of these structures has been studied extensively [29]. Waveguiding arises as an application of photonic crystals whenever symmetry of the lattice is broken, for instance when a line of rods is removed from the crystal. In this scenario, if the crystal is illuminated with a wave whose frequency is in the gap, the wave will only be allowed to travel along the defect and will decay exponentially away from the defect. The MSCG presented in this paper is therefore an attractive candidate to simulate photonic crystal applications, since by solving a small number of subproblems it enables the simulation of large lattice structures by exploiting the repeated patterns. The first example is a waveguide T-splitter consisting of Gallenium Arsenide rods (ε = 11.4) in air of radii R 0 = 0.4a, where a is the periodicity of the crystal. For the TM polarization, this structure presents a first bandgap for ω ∈ (0.36, 0.40) and a second bandgap for ω ∈ (0.52, 0.55). In order to numerically simulate the splitter with the MSCG method, we first have to identify the subdomains in which to split the computational domain. The subdomains should be invariant to rotation and translation to ensure that only a small number of subproblems are solved. The size of the subdomain along with the frequencies of interest determine the polynomial order chosen to approximate the solution at the global interfaces. Conversely, the discretization needed at the subdomain level is governed by the details of the geometry. Furthermore, the MSCG approach offers the possibility of using different polynomial approximation order for different subdomains, which can result in greater efficiency.
For the waveguide splitter shown in Fig. 3 (left), we choose three different classes of subdomains: (1) the rod subdomain, where we use p = 2, curved elements and adaptive mesh size to approximate the fine detail of the rod curvature; (2) the air subdomain, where we use an homogeneous mesh of straight-sided elements of order p = 2; and (3) the PML subdomain, with the same discretization as type 2. Although not represented in Fig. 3, three PML subdomains surround the computational domain in each direction to ensure outgoing waves are properly attenuated. The subdomain discretization is chosen through a grid convergence study on progressively refined meshes, until the relative error of the optical intensity at both output ports is below 0.1%, compared to an extremely fine subdomain discretization. For each global interface, we choose a single element of polynomial order p f = 10, which gives a sufficient resolution of more than 20 points per 14 wavelength. The degree of freedom count is reported in Table 3. Note that solving the exact same problem with regular CG would suppose solving a linear system of more than 6M dof as opposed to a system of 26K dof, thus for structured geometries the MSCG is a very competitive approach.  In Fig. 3 (right) we show the amplitude field E z for a frequency in the first bandgap, a frequency between the two bandgaps and a frequency in the second bandgap. An attractive feature of waveguiding with photonic crystals are the low losses that occur even for sharp bends, thus enabling the efficient manipulation of electromagnetic waves.

Waveguide Z-bend optimization
Here we consider the simulation and design of a photonic crystal consisting of a silicon (ε = 12.1) slab with drilled air holes of radius R 0 = 275a/800 forming a lattice with triangular symmetry. This structure presents a broad bandgap for ω ∈ (0.26, 0.34) for TE waves, see [17]. A waveguide is generated by opening a Z-shaped defect and illuminating the crystal with a wave impinging at the left input port. The quantity of interest is the intensity of the optical power at the right output port in the x-direction, namely where the subdomains at the output port Ω out comprise the line defect and one hole subdomain above and below. The schematics of the slab are shown in Fig. 4a (left), where we have ensured that the bends are sufficiently separated such that they do not interact. The objective here is to apply the methodology introduced above to design the bending region such that transmission is enhanced for certain frequencies of interest within the bandgap of the regular lattice [17], see Fig. 4b.
The different subdomain classes identified are shown in Fig. 4a (right). The subdomain classes 1 and 2 correspond both to the hole subdomain, but are treated independently since class 1 will be subject to shape optimization, whereas the subdomains from class 2 remain unaltered. The other two subdomain types correspond to the line defect and the PML respectively. The dof count for the different types and the global problem is collected in Table 4.  In the next sections, we present optimization results for three different scenarios: single frequency, range of frequencies and range of frequencies combined with the presence of geometric imperfections in the holes. In  The optimizations are performed using the nlopt [31] optimization package. For each objective function, we first apply the multi-level single-linkage global optimization algorithm, and then choose the computed optimum as the starting point for a local gradient-based optimization using the preconditioned truncated Newton, in order to perform a better local exploration, see references in [31]. Relative tolerances for both the objective function value and the design parameters are employed as termination criteria, and are set to 10 −2 for the global optimization and 10 −4 for the local optimization. Finally, for each objective function we compute multiple global-local optimization cycles with different initial guesses, in order to avoid getting stuck in local minima. The results shown here correspond to the optimization cycle that achieved the lowest objective function value.

Single frequency radii optimization
The first step towards efficient optimization is to develop a RB for subdomain class 1, which reduces to finding a RB for the Dirichlet subproblem (5b) as we have no source. The mapping for subdomain class 1 is described in Appendix A, whereas the mapping for the other subdomain classes is the identity. The application of the EIM on the three components of G (represents a second-order symmetric tensor) and the Jacobian g enables the affine parametrization of the mapping, for a total of Q G = 6 and Q g = 4 interpolation elements. We then compute snapshots to (5b) for the different Dirichlet conditions, where the parameters are θ and the frequency ω ∈ [0.3045, 0.3055] -thus the RB can be reused for the next example. Finally, the RB is constructed by compressing the snapshots for all the boundary conditions to form a single POD basis.
We now proceed to optimizing the radii of the holes for a single frequency ω = 0.305, for which intensity at the outlet is poor. For subdomains 2, 3 and 4 we use CG to compute (5b)-(6) in 0.42, 0.15 and 0.30 seconds respectively. Thus for each new θ we only need to solve for the 58 subdomains from class 1. The deterministic optimization problem reads where the expensive CG local solution (5b)-(6) is substituted for its inexpensive RB counterpart (9)- (10) to accelerate the computations. We use a reduced basis of size N = 100, rendering a wall time for all 58 suproblems of 0.12 s, as opposed to 25.60 s required by the high-fidelity CG method. The global system is finally solved in 0.47 s. Overall, the reduced basis approach delivers the value of s N in just 1.45 s, that is 18 times faster than computing s h with CG on all subdomains, and without compromising accuracy as shown below.
Convergence is achieved after 1200 iterations. A frequency sweep for s N (θ * ) is shown in Fig. 5a (dashed red), showing a much higher transmission than the un-optimized version, see Fig. 4b. Lastly, accuracy of RB is verified by re-optimizing (16) with s h instead of s N , setting the maximizer θ * as initial guess. The new optimum θ * h is attained only after 3 iterations, since s N is a faithful representation of s h . The relative error of the output for the same maximizer, s h (θ * ) − s N (θ * ) /s h (θ * ) = 5.8e − 4, as well as the relative error of the maximizer, θ * h − θ * 2 / θ * h 2 = 8.9e − 4, are below the acceptable engineering standards, therefore concluding the RB effectively approximates the high-fidelity CG solutions of the local subproblems.

Frequency-range radii optimization
The next step is to seek a design that is robust with respect to a range of frequencies, which is desirable when the optimum presents a sharp peak with a rapid decay away from the single frequency optimum. This optimization problem can be recast in a stochastic optimization framework as where γ controls the weight assigned to the variance minimization. Since the stochastic space is unidimensional, the output statistics can be computed with simple Gauss-Legendre quadrature, reducing (17) to the weighted evaluation of the output s N at selected frequencies given by Legendre points. We choose a 20-point quadrature that integrates exactly polynomials up to degree 39, selected as the coarsest quadrature that renders a mean relative error on (17) of < 0.1% for multiple θ values, compared to the results of a 100-point quadrature. We reuse the RB constructed in the previous section for subdomain class 1, again with 100 modes, and precompute (5b)-(6) for subdomain classes 2, 3 and 4 at the selected Legendre frequencies in the range ω = [0.3045, 0.3055] to expedite the optimization process. A frequency sweep for s h ( θ) (dashed green) with γ = 1 is shown in Fig. 5a. The radii configuration θ is indeed less sensitive to frequency variations in the prescribed range, albeit attaining a lower output intensity at ω = 0.305 compared to the single-frequency optimum θ * (dashed red).
Similarly as before, we can assert the RB is accurate by inspecting both the relative error on the optimum if the objective function is evaluated with CG ( s h ) instead of RB ( s N ) for the subproblems belonging to class 1, namely s h ( θ) − s N ( θ) / s h ( θ) = 1.8e − 4, and the relative error of the optimum θ h when maximizing s h initialized with θ, that is θ h − θ 2 / θ h 2 = 7e − 4.

Frequency-range robust radii optimization
Finally, we analyze the robustness of the solutions with respect to geometry errors on a range of frequencies.
The geometry errors, which will be considered for subdomains in classes 1 and 2, are given by (A.1), where we append the constant deformation θ m for the optimization, namely r m 0 → R m 0 [1 + θ m ] + δR m 0 . For the results below, we select D = 11, σ = 0.02, L c = 1/16 and z m ∈ [− √ 3, √ 3] D in (A.1), for which the 95% confidence interval gives δR m 0 /R m 0 < 3% if θ = 0. Note we encapsulate both the optimization and the stochastic parameters in the same expression, and thus the same RB for both classes of subdomains can be used by simply setting θ m = 0 for the subdomains in class 2.
In order to achieve robust designs accounting for variation in both the frequency and the geometry parameters, we formulate the following stochastic optimization problem The stochastic dimension of the problem under consideration is large, since we have 11 geometry parameters per subdomain, and there are 665 subdomains of classes 1 and 2 combined. Hence, we shall resort to the MVR method described in Section 3.2 to evaluate the objective function (18). We first develop a RB for the Dirichlet subproblem combining the EIM with POD for the non-homogeneous radius variation. We set θ m ∈ [−0.127, 0.047]a as before for the optimization, which in this case leads to Q G = 61 and Q g = 23 elements in the interpolation basis. We again use N = 100 for the RB model, but since have considered 11 additional parameters to represent the geometry we can expect this RB to be less accurate than the one developed for the homogeneous radii variation. The main difference here is that the RB is seen as a surrogate that correlates with the high-fidelity model, not as a substitute, thus a coarser basis still casts excellent results.
Subdomains 3 and 4 are solved with CG in 0.15 and 0.30 seconds respectively, and the global system in 0.47 s. Unfortunately, the CG computation for subdomain class 2 can no longer be precomputed and reused due to the geometry errors, thus requiring 665 high-fidelity CG linear system solutions for subdomains 1 and 2 combined. This calculation takes 301.3 s in our implementation, for a total MSCG wall time of 302.2 s for a single tuple (θ, z). Additionally, computing the output statistics at each optimization step typically demands multiple output evaluations for different z, thereby greatly hindering the computational feasibility of the optimization. To circumvent this issue we employ the RB on subdomains 1 and 2, reducing the wall time to 3.22 s for a total MSCG time of 4.16 s. The MVR method combines the fact that s N is more efficiently evaluated than s h (a factor of 70), with high variance reduction between outputs V (s h − s N )/V (s h ) > 50, therefore conveniently allocating many more samples to the RB output than to the variance reduction term (I 1 = 15K and I 0 = 355 for ∆ E I0,I1 = 10 −3 at the optimum) in order balance the computational cost of (12). In Fig. 5a we show in dashed line, for a set of discrete frequencies in the interval ω = [0.3045, 0.3055], the transmission s h for the several configurations considered: single frequency optimization θ * , frequency-range θ and frequency-range robust design θ. To assess the robustness with respect to geometry, for each design we introduce geometry variations on the subdomain classes 1 and 2, and evaluate the expected value of the output and its 95% confidence interval as E G [s h (θ)]±∆ E in Fig. 5a, as well as the variance V G [s h (θ)]±∆ V in Fig. 5b, both represented in solid line. The single frequency optimum θ * produces the highest transmission for ω = 0.305, but it degrades significantly in the presence of geometric errors. Conversely, the other optima, despite achieving lower peak outlet intensities, maintain a consistent performance for all the frequencies in the range. Moreover, the robust optimum θ outperforms the other two optima both in expected value and in variance, see Fig. 5b. These results show the importance of accounting for geometry variations in the objective function if robust designs with respect to geometry are sought. Finally, the optima configurations for the three scenarios considered are depicted in Fig. 5c, expressed as a variation on the nominal radius.  (c) Single frequency θ * (left), range of frequencies θ (middle) and robust design θ (right) optima.

Conclusions
We have presented a multiscale continuous Galerkin method for computing wave propagation phenomena in photonic crystals. The method relies on partitioning the computational domain into subdomains, and computing the local solution in the subdomains as a function of the Lagrange multipliers at the interfaces. Thus, only a linear system for the global Lagrange multipliers is solved. The solution in the interior of the subdomains is recovered after the global problem has been solved. Furthermore, the MSCG method is especially advantageous for problems that exhibit repeated patterns, since the local subproblems only need to be solved once for each type of subdomain. In addition, we propose the construction of a reduced order model for the local subproblems that allows us to consider geometric variation within the subdomains, leveraging a reference domain formulation of the governing equation. We then use the model and variance reduction method to compute statistical outputs of stochastic wave propagation problems. Finally, we verify the convergence rate of MSCG, carry out a deterministic simulation of a waveguide splitter to illustrate the advantages of MSCG in photonic crystal, and a robust optimization of a photonic slab to demonstrate the performance of the MSCG method to attain robust designs.
We now describe the mapping G m used for subdomains containing a dielectric rod. We consider square subdomains as shown in Fig. A.6b, and require that the subdomain boundaries remain fixed, as well as a small box surrounding the origin (to prevent the mapping becoming singular). Since we are only considering deformation in the radial direction, the mapping can be expressed as (x, y) = G m (x r ) with x = r(x r ) cos α and y = r(x r ) sin α.
In order to define r(x r ) for the interior points in the domain, we assume that the perimeter of the rod is deformed according to (A.1), and that the remaining points are linearly deformed according to The derivatives ∂r/∂x i may be computed with the chain rule. Note there are two distinct types of boundaries, see Fig. A.6a, and its distance to the origin can described as B = 1 + (Y /X) 2 (red) and B = 1 + (X/Y ) 2 (black). The distance B to the boundary of the origin box is similarly computed.