Consider the scattering of a time-harmonic electromagnetic plane wave by a periodic structure, all known as a grating in diffractive optics. The scattering theory in periodic structures has many applications, particularly in the design and fabrication of optical elements such as corrective lenses, antireflective interfaces, beam splitters, and sensors [15, 16]. Throughout, it is assumed that the medium is invariant in the -direction. Thus the model problem of the three-dimensional time harmonic Maxwell equations can be reduced to a simpler model problem of the two-dimensional Helmholtz equation. The more complicated biperiodic problem will be considered in a separate work.
The scattering problems in periodic structures have been studied extensively on both mathematical and numerical aspects. We refer to  and references therein for the mathematical studies of existence and uniqueness of the diffraction grating problems. Numerical methods can be found in [13, 14, 4, 25] for either an integral equation approach or a variational approach. A comprehensive review can be found in [27, 5] on diffractive optics technology and its mathematical modeling as well as computational methods.
This paper is concerned with the numerical solution of the inverse problem, which may be described as follows: given the incident field, determine the periodic structure from a measured reflected field a constant distance away from the structure. The inverse problem arises naturally in the study of optimal design problems in diffractive optics, which is to design a grating structure that gives rise to some specified far-field patterns [15, 16, 3]. A number of numerical methods have been developed to solve these inverse problems. Ito and Reitich  proposed a high-order perturbation approach based on the methods of variation of boundaries. In , Arens and Kirsch applied the factorization method to scattering by a periodic surface. Iterative regularization methods were developed by Hettlich in  based on the shape derivatives with respect to the variations of the boundary. Bruckner and Elschner  gave a two-step optimization algorithm to reconstruct the grating profile. See also Elschner  for a reconstruction algorithm of a two-dimensional periodic structure based on finite elements and optimization techniques. More recently, Bao  presented an efficient continuation method to capture both the macro and micro structures of the grating profiles with multiple frequency data. The method was further extended to the case of phaseless data in 
. Related uniqueness results and stability estimates for the inverse diffraction problem were obtained in[23, 2, 7, 11].
Existing studies of the inverse problem mostly assume that the periodic structure is deterministic and only the noise level of the measured data is considered for the inverse problem. In practice, however, there is a level of uncertainty of the scattering surface, e.g. the grating structure may have manufacturing defects or it may suffer other possible damages from regular usage. Therefore, in addition to the noise level of measurements, the random surface itself also influences the measured scattered fields. Considering the inverse diffraction problem by random surfaces is closer to reality. In fact, surface roughness measurements are of great significance for the functional performance evaluation of machined parts and design of micro-optical elements. In this scenario, the analytical methods are insufficient any more to deliver desired error tolerances. Little is known in mathematics or computation about solving inverse problems of determining random surfaces. One challenge lies in the fact that the scattered fields depend nonlinearly on the surface, which makes the random surface reconstruction problem extremely difficult. Another challenge is to understand to what extend the reconstruction could be made. In other words, what statistical quantities of the profile could be recovered from the measured data? Until now, only some initial progress has been made. In , an asymptotic technique has been applied to obtain the analytical solutions with reasonable accuracy in the case of a small perturbation of a plane surface. B-slines have been employed in  to recover the unknown rough surface, which is represented by the control points for the spline function. The Monte Carlo method in  has been adopted to calculate the statistics of the scattered waves. Shi  have used an ultrasonic methodology to evaluate the correlation length function of the randomly rough surfaces in solids.
In this work, we propose an efficient numerical method to reconstruct the random periodic structure from mutli-frequency scattered fields away from the structure, in the sense that three critical statistical properties, namely the expectation, root mean square, and correlation length of the random structure may be reconstructed. Our method (MCCUQ method) is based on a novel combination of the Monte Carlo technique for sampling the probability space, a continuation method with respect to the wavenumber, and the Karhunen-Love expansion of the random structure. Throughout, the periodic structure is assumed to be a stationary Gaussian process, which is reasonable for many practical situations. For the representation of the random structure, we employ the Karhunen-Love expansion in this work , due to the fact that the expansion is optimal in the sense that the mean-square error of the truncation of the expansion after finite terms is minimal . Also, our work may be viewed as an extension of Bao , where an effective continuation method was first introduced to reconstruct the (deterministic) periodic structure from the multiple frequency data.
The rest of the paper is organized as follows. The model problem is introduced in Section 2. The Monte Carlo - Continuation - Uncertainty Quantification reconstruction method is presented and discussed in Section 3. Several numerical results are presented in Section 4 to demonstrate the accuracy and reliability of the method. The paper is concluded with some general remarks.
2 Model problem
Consider the scattering of an incident field by a periodic random surface
where the rough surface is a stationary Gaussian process as shown in Figure 1. Here, for a probability space , denotes the random sample, are the spatial variables, and the random surface is the sum of a deterministic funciton and a stationary Gaussian process with a continuous and bounded covariance function . The deterministic function is further assumed to be -periodic, that is, , and the covariance operator of the random surface is defined by
For the representation of the random structure, we choose the Karhunen-Love (KL) expansion[18, 30]. The KL expansion of a square-integrable stochastic process is a particularly suitable spectral decomposition for uncertainty quantifications since the mean square error of the truncation of the expansion after finite terms is minimal. In addition, the expansion with respect to the orthonormal basis enjoys many useful properties. Specifically, we have
where is a
-periodic deterministic function. The orthonormal eigenfunctions of the covariance operatorare defined by
where the corresponding eigenvaluesare arranged in a descending order, and
is a random variable with zero mean and unit covariance. It follows that the KL expansion of the random processmay be written as
where , and are independent and identically distributed random variables with zero mean and unit covariance.
where is the root mean square of the surface, and is the correlation length.
We next state the model problem. Consider a plane wave incident on a periodic structure which is characterized by the wavenumbers ruled on a perfect conductor. It is assumed that the medium and material are invariant in the direction, and periodic in the direction with a period . Above the structure, the medium is assumed to be homogeneous. In this two-dimensional setting, for simplicity, we consider the transverse electric (TE) polarization case.
More precisely, let
be an incident plane wave , where , is the angle of incidence with respect to the positive -axis. Let the total field
where is the scattered field and . Then, the time-harmonic Maxwell equations can be reduced to the two dimensional Helmholtz equation 
In addition, as usual, due to the uniqueness consideration, we seek for the quasi-periodic solutions which satisfy
Moreover, since the medium below is perfectly electric conducting, from physics, the following surface condition holds
For each given sample , it follows from the Rayleigh expansion that can be written as
with the boundary given by
This work is to study the inverse problem of reconstructing a perfectly reflecting random periodic structure.
Inverse Problem: For each sample , given the incident wave , determine the KL expansion random structure as well as the two important statistical properties of the structure, the root mean square and the correlation length , from the multiple frequency measurements of the scattered fields .
3 Reconstruction method
For each given sample , we begin with the single-layer potential representation for the scattered field
where is an unknown periodic density function and is a quasi-periodic Green function given explicitly as
Since is quasi-periodic, we can expand
On the other hand, the periodic function with periodicity can be expanded as
Therefore, we have
Define the operator :
In practice, in order to restrict the exponential growth of noise, we take the regularization  here, that is,
where is some positive regularization parameter.
3.1 The nonlinear problem
Recalling that in Eq. 12, we may study the nonlinear problem
Substituting the expansion for the operator , we have
Since the summation decreases exponentially with respect to , we choose a sufficiently large , truncate it into finite summation, and consider the numerical solution of the following nonlinear equation
3.2 The Monte Carlo technique
In order to approximate the solution of the nonlinear equation Eq. 23, the Monte Carlo technique is employed to sample the probability space, which allows us to derive useful statistical properties for the reconstruction.
Let be a positive integer which denotes the number of realizations. For each sample , we denote the reconstructed random surface by . Thus the algorithm yields an approximation of given by
Similarly, the standard deviation of the sampled solution can be computed by using the formula
the expectation of which is an approximation of the root mean square of the random structure.
when the sample size is large enough, and the profile is a real -period function. Thus the mean value admits a Fourier series expansion. Without loss of generality, we take the period to be from now on. Hence has the following expansion
Therefore, the first step of our method is to determine each coefficient to reconstruct the mean structure profile of the random surface.
3.3 The Monte Carlo - continuation (MCC) method
Because of the nonlinearity of Equation Eq. 23, we propose a Monte Carlo - continuation (MCC) method to recursively reconstruct these Fourier coefficients for all samples.
It is evident that the finite expansion of Equation Eq. 26 can reasonably approximate the random periodic grating profile. Here, we choose a prescribed wavenumber not exceeding the wavenumber , and the number of Fourier expansion modes is taken not larger than the prescribed wavenumber.
For each sample , we set the initial approximation and , where
is taken to be the largest integer that is smaller than or equal to the prescribed wavenumber. Denote the vector
We first choose an initial value for the wavenumber , and seek for an approximation to the profile by the Fouries series
where is taken to be the largest integer that is smaller than or equal to the wavenumber . Denote the first terms of , i.e. with and .
For the incident angle , , , define
Denote , then the nonlinear equation Eq. 23 can be reformulated as
where . In order to reduce the computational cost and instability, for each sample, we consider the nonlinear Landweber iteration 
where is a relaxation parameter dependent on the wavenumber, the initial vector , and the Jacobi matrix
can be computed explicitly. In order to update the reconstructed grating profile function by recursive linearization, we set the resulting coefficients as . Namely, update the first terms of by using the above nonlinear Landweber method at the wavenumber .
Next, by increasing the wavenumber to (), we seek for a new approximation to the profile by the Fourier series
and determine the coefficients , where is again taken to be the largest integer that is smaller than or equal to the wavenumber . Denote where
and we repeat the above Landweber iteration with the approximation to the profile at the previous smaller wavenumber as our starting point, i.e.
We may repeat the above procedure until the prescribed wavenumber (smaller than ) is reached. The resulting coefficients are deduced by taking the expectation of , i.e.
Therefore, the mean grating profile of the random periodic structure is reconstructed as the following finite Fourier series expansion
Next, we discuss some key practical implementation issues of our MCC method. On the one hand, we use multi-frequency scattered field data for a range of incident angles to realize the reconstruction of random periodic structure. At each recursion, the initial value of the high frequency reconstruction depends on the previous lower frequency information. Our numerical experiments have shown that the present MCC method converges for a larger class of rough surfaces than the usual Newton’s method using the same initial value. On the other hand, we take the relaxation parameter dependent on the wavenumber. In fact our algorithm converges for a larger class of surfaces, provided that the relaxation parameter is chosen to be smaller as the frequency increases.
Our algorithm seems to be convergent and stable. As the result, a good approximation to the random periodic structure at a higher frequency can be obtained at each recursion. The precise description of our method to reconstruct the mean profile is given in Algorithm 1. It should be pointed that since all of the samples in the probability space are independent and identically distributed, the Monte Carlo sampling process can be performed well in parallel.
3.4 Monte Carlo - continuation - uncertainty quantification (MCCUQ) method
Based on the mean grating profile reconstructed in Section 3.3, we can realize the reconstruction of the random process including the two key statistics characterizing the random structure. Recalling that
we have is a stationary Gaussian process.
By simple calculations, we can see that if the eigenvalues of the covariance function are in a descending order, then we only need the first several terms to represent the random profile with an accuracy of . For example, when the period is , the root mean square is , and the correlation length is , Fig. 2 shows the decay of the eigenvalues.
Hence, we can truncate the KL expansion and further reconstruct the profile of the random surface only considering the first several terms. Since the eigenfunctions of the covariance operator are orthonormal and are independent and identically distributed random variables with zero mean and unit covariance, the eigenvalues of the covariance function are just that of the corresponding covariance matrix where
Therefore, the second step of our reconstruction method is given as follows (Algorithm 2). Note that the grating profile of a random periodic structure is reconstructed in the form of the sum of a finite Fourier series expansion and a truncated KL expansion represented by the reconstructed eigenvalues by this MCCUQ method.
Since the correlation length is greatly smaller than , it can be deduced after simple calculation that the eigenvalues are only related to the root mean square and the correlation length satisfying
Thus we can further approximate these two statistics which characterize the random surface only by the eigenvalues reconstructed above. Specifically, we just need the first two eigenvalues and , which leads to the following recovery formula
In fact, the correlation length and root mean square may be determined by any two eigenvalues.
4 Numerical results
To test the stability of our method, some random noise is added to all the measured scattering data, i.e., for each given sample, the scattering data takes the form
where rand represents uniformly distributed random numbers inand is the noise level of the measured data. The near-field measurements in this work are simulated by using an adaptive finite element (AFE) method with a perfectly matched absorbing layers[14, 4]. Since each sample is independent and identically distributed, our numerical method can be performed in parallel.
For all of the numerical examples presented below, the number of eigenvalues determined by the covariance operator is chosen to satisfy the random surface with an accuracy of ; the noise level and the regularization parameter are taken as and respectively; the relaxation parameter dependent on the wavenumber at each stage is chosen as ; and the truncation of the infinite summation in Equation Eq. 23 is set to . Finally, we set the number of realizations in sampling the probability space.
4.1 Accuracy of the algorithm
Recalling that the random surfaces are modeled by the KL expansion with the covariance function
to test the accuracy of our algorithm, we first show three realizations of randomly rough surfaces in one period when , the correlation length , and the root mean square , respectively, as seen in Fig. 3. It is clear that the deviation of the sample profile from the center profile becomes larger as the root mean square increases.
To reconstruct a random periodic surface that is a sum of a deterministic function
and a stationary Gaussian process with the covariance function
with the correlation length . Graphs of the original profile and the reconstructed profiles with different wavenumbers are shown in Fig. 4 with and , respectively. Since the original deterministic profile consists only of a few Fourier modes, only a small number of iterations are needed to get a good reconstruction both for and .
To reconstruct a random periodic surface that is a sum of a deterministic function
and a stationary Gaussian process with the covariance function
again with the correlation length . Graphs of the original profile and the reconstructed profiles with different wavenumbers are shown in Fig. 5. As the original deterministic profile contains more Fourier modes, incident waves with higher frequencies are needed to get a better resolution of the reconstruction.
Recalling the formula of the expectation and standard deviation given in Section 3.2, we denote as the deviation of the reconstruction profile and the original deterministic profile for each sample. The mean error and standard deviation are denoted by and respectively. In Example 1, the mean error and standard deviation are when the root mean square , and when the root mean square . In Example 2, when the root mean square , and when the root mean square . Thus as shown in Fig. 4-Fig. 5, all the reconstructions at least have an accuracy of by our method (Algorithm 1 and Algorithm 2) proposed in this paper. Also, we can further evaluate the root mean square of the random structure directly by the standard deviation formula.
Therefore, we can reconstruct the unknown random periodic structure with an accuracy of and give an evaluation of the root mean square value to certain degree, which shows a better reconstruction than that of the deterministic grating profile in .
4.2 Performance of the correlation length
We now discuss the influence of changes of the correlation length on our proposed method. For this, we consider random periodic surfaces with different correlation lengths. Fig. 6 shows the realizations of random periodic surfaces when and respectively. It is noticeable that the surface becomes rougher as the correlation length decreases.
We still consider the random periodic surfaces in Example 1 and Example 2 mentioned above. By simulation, it is deduced that all the reconstructions of the random periodic surfaces for the correlation length and the root mean square in Example 1 at the wavenumber and in Example 2 at the wavenumber have an accuracy of . That is, our proposed algorithm has the convergence and stability for both the rougher profile with and the larger root mean square with , which shows that the method proposed in this work is reliable and efficient.
4.3 Evaluation of the statistics of the random periodic surface
We still consider the random periodic surfaces in Example 1 and Example 2 mentioned above. Fig. 7 shows the reconstruction of eigenvalues in Example 1 when and respectively. Fig. 8 shows the reconstruction of eigenvalues in Example 2 when and respectively. As shown in Fig. 7-Fig. 8, all the eigenvalues reconstructed by our MCCUQ method achieve at least an accuracy of . By doing more experiments, we observe that with the smaller root mean square or the larger correlation length, the reconstructed eigenvalues can reach a higher accuracy.
From the reconstructed eigenvalues, we can calculate the statistics of the random periodic surface by the recovery formula discussed in Remark 1. It can be calculated directly that , when , and