Wave propagation in random waveguides

We study uncertainty bounds and statistics of wave solutions 
through a random waveguide which possesses certain random 
inhomogeneities. The waveguide is composed of several homogeneous 
media with random interfaces. The main focus is on two homogeneous 
media which are layered randomly and periodically in space. 
Solutions of stochastic and deterministic problems are compared. The 
waveguide media parameters pertaining to the latter are the averaged 
values of the random parameters of the former. We investigate the 
eigenmodes coupling due to random inhomogeneities in media, i.e. 
random changes of the media parameters. We present an efficient 
numerical method via Legendre Polynomial Chaos expansion for 
obtaining output statistics including mean, variance and probability 
distribution of the wave solutions. Based on the statistical 
studies, we present uncertainty bounds and quantify the robustness 
of the solutions with respect to random changes of interfaces.


1.
Introduction. Many classes of problems in computational science and engineering are characterized by uncertainties. Even if the mathematical equations are a perfect representation of physical reality, its solutions will not describe reality unless we have perfect knowledge about the initial and boundary conditions which are often difficult to prescribe with high accuracy. Some uncertainties are inherited from the data, or from a lack of knowledge of the underlying physics or from the impossibility of fully describing the solutions of a system possessing a large number of degrees of freedom (see e.g. [1], [15], [18]). The goal of uncertainty quantification is to investigate the impact of errors, or uncertainty, on the solutions. These uncertainties come from incomplete knowledge of physical parameters, initial and boundary conditions, location of media interfaces in electromagnetic and acoustical waveguides, errors in measurements and other sources. To predict or control solutions from governing equations or to quantify uncertainty in the solutions, we 148 CHANG-YEOL JUNG AND ALEX MAHALOV should take uncertainty quantification into account in order to fully understand simulation results and subsequently the true physics. Problems with uncertainties include Stochastic Partial Differential Equations (SPDE) (see e.g. [11], [6], [3], [4]) and PDE supplemented with uncertain data (initial/boundary conditions, parameter values) etc. (see e.g. [29], [26], [30], [2], [14], [16]). In this paper we investigate statistical properties of the solutions of hyperbolic equations in randomly layered media characterized by the presence of random interfaces. In different contexts, i.e. different type of equations, you may refer to [7], [23], [24], [25], [27], [22] and [32].
We aim to construct a method for obtaining output statistics of solutions of onedimensional wave equation in a periodic waveguide pertaining to the random media parameter c: where the initial conditions g, h are deterministic and 2π-periodic smooth functions, and the parameter c, which depends on the media properties, is defined with uncertainties. More precisely, c = c ξ1,ξ2 (x) = c 0 for x ∈ (0, ξ 1 ) ∪ (ξ 2 , 2π) and c 1 for Fig. 1), and ξ 1 , ξ 2 are random parameters which are defined in (3) below. The media, in particular, waveguide, can be polluted or damaged for some reasons and then the properties of some parts of media might change. However, in general, such changes or location of imperfections and impurities are uncertain. We thus would like to compute the output solutions taking into account such uncertainties. The equation (1) is not an energy conservative form. However, the work presented here can be extended to a conservative form, i.e. u tt = (c 2 u x ) x .
A waveguide, e.g. optical fiber, is a material or structure which guides electromagnetic or sound waves and they vary from slab waveguides to fiber or channel ones. Although our model problem is posed on one-dimensional periodic space, the numerical techniques which take into account randomness of waveguide media properties can be extended to 2 or 3 dimensional spaces in which real waveguides can be constructed.
Let us first consider two kinds of homogeneous media characterized by the parameters c 0 , c 1 , resp., which interface randomly and periodically in space. We assume that the random parameters ξ 1 , ξ 2 are independent and satisfy and we also assume that ξ 1 , ξ 2 are uniformly distributed over (l 1 , l 2 ), (l 3 , l 4 ), respectively. We can then write probability density functions f 1 (ξ 1 ), f 2 (ξ 2 ) as where χ A (ξ) is the characteristic function of A. Random parameters with nonuniform distributions will be discussed in future studies. We approximate u as follows where P l1,l2 k , P l3,l4 l are shifted Legendre polynomials such that and P k are the standard Legendre polynomials with span (−1, 1). We call (4) Legendre polynomial chaos (PC) expansions and the coefficients u ij (x, t) PC mode (see [9], [10], [31]). Thanks to the orthogonality of Legendre polynomials, we can explicitly obtain the expectation and variance. That is, (u ij (x, t)) 2 (2i + 1)(2j + 1) .

Remark 2.
Although the data g(x), h(x) are deterministic, their stochastic versions can be also taken into account by writing (4) for u = g, h and then the vector of their PC modes is replaced for the right-hand sides of (8) 3,4 .
So far we have considered two random interfaces. In the next section we demonsrate how to reduce the computational size by decomposing the parameter c = c ξ1,ξ2 (x) and we will also consider n-random interfaces below.
Using the Kronecker product we can write where Noting that where U = U (x, t) = (u ij ) and U S is the stacking of U and we can then rewrite (8): Using the fact (B T ⊗ A)X S = (AXB) S and noting that L 1 (x), L 2 (x) are symmetric we can write and thus The computational size can then be substantially reduced by writing that, for 0 ≤ Comparing with (8) the matrix size involved is reduced from (N 1 N 2 ) 2 (for Λ(x)) to either N 2 1 or N 2 2 (for L 1 (x), or L 2 (x)). The numerical plots in Figs. 2-3 show respectively mean, variance, PC modes and the decay of PC modes. We observe that the PC modes u ij , i, j = 0, 1, · · · , 10, decay exponentially and they decay faster with small variance (or fluctuation) than with large variance. The PC modes in the Figures are computed from the matrix equation (17) or (18)   the initial conditions u = v, u t = 0 at t = 0 where v is an eigenvector given by (24) below with the eigenvalue ω 2 = 1.7444405 for (a), 2.0806212 for (b).

Remark 3.
It is noteworthy to investigate the entries of Λ(x) as either l 1 or l 2 (similarly, either l 3 or l 4 ) is close to x. We first observe that for x ∈ (l 1 , l 2 ), setting Hence we can see that . This implies that Λ(x) tends to be diagonalized as the lengths l 1 , l 2 , l 3 , l 4 tend to be small, i.e. the variances of ξ 1 , ξ 2 approach zero. Here we infer that as the variances of ξ 1 , ξ 2 are getting larger, the matrix Λ(x), x ∈ (l 1 , l 2 ) ∪ (l 3 , l 4 ), tends to be dense and the PC modes broadly interact each other in the x-space.

2.
Robustness of wave solutions with respect to random changes of waveguides interfaces. In this section we investigate some statistical properties of the wave solutions of (1) and the robustness of the solutions with respect to random changes of interfaces. For this, we compute the stochastic solutions of (1), and as in Figure 4 we compare their averages u 00 with the deterministic solutions of (1) with c 2 = µ 2 (x) = E(c ξ1,ξ2 (x) 2 ) (see Fig. 1), i.e. (35) below. In the deterministic problem we consider a special solution, i.e. an eigenmode u * = e iωt v(x), where v is an eigenvector (or carrier) of the (deterministic) Helmholtz equation (see [17], [5]):  Note that u * is the solution of (35) below with initial conditions u * = v(x) and u * t = 0 at t = 0. Comparing with the stochastic solutions of (1) with the same initial conditions we investigate the robustness of the eigenmode u * through the waveguides. The eigenvalue ω 2 can be computed from a discrete problem described in (39) below.
We first obtain the expectation and variance of c ξ1,ξ2 (x) 2 in the lemma below.

Homogeneous media.
It is interesting to analyze the behaviors of solutions in homogeneous media, i.e. the media parameters are constant. In particular, we consider the averageμ as in (32) below which does not depend on x and the stochastic parameters ξ 1 , ξ 2 . By elementary calculations, we can evaluate the averagē µ 2 over a space period: We first compare the solutions of two deterministic problems (1) with c = µ(x) and c =μ respectively and obtain the error between the two solutions. More explanations follow in Remark 5 below. Lemma 2.2. Let u * be a solution of the deterministic problem (1) with c 2 = µ(x) 2 = E(c ξ1,ξ2 (x) 2 ) andū of the deterministic problem (1) with c 2 =μ 2 = E(c ξ1,ξ2 (x) 2 ). Then we have the following estimates for u * −ū: Proof. Setting w = u * −ū we can write the equation for w: w tt = µ(x) 2 w xx + (µ(x) 2 −μ 2 )ū xx . Multiplying by w t and integrating over Ω, thanks to Poincaré inequalty, see e.g. [5], we find that Noting µ 2 (x) ≥ d 2 > 0 and from (25) we find that and applying the Gronwall inequality the estimate (33) follows.
Remark 5. The deterministic problem (1) with c = µ(x) (inhomogeneous media) (see the plot of µ 2 (x) in Fig. 1) corresponds to the problem (8) with N 1 = N 2 = 0, "mean dynamics". From (8) we can write The robustness of the stochastic solutions with respect to randomness can be measured by comparing the two mean solutions u 00 of (8) respectively with N 1 = N 2 = 0 (i.e. (35)) and N 1 , N 2 sufficiently large. The numerical results are shown in Figure 4 for various eigenvectors and variance (fluctuation) of the random interfaces. The two solutions are averaged outputs; the first one is a deterministic solution with averaged stochastic parameters ξ 1 , ξ 2 and the second one is an average of the stochastic output solution u of (1). For homogeneous media, i.e. c =μ, thanks to d'Alembert 's formula, the explicit solutions are given bȳ The errors between u * andū as in (33) imply that the error is only up the inhomogeneities which introduce coupling of eigenmodes. The eigenmode coupling due to random inhomogeneities is described in the next section.

2.2.
Inhomogeneous media due to random inhomogeneities. Due to inhomogeneities in media (note c = c ξ1,ξ2 (x)) we observe eigenmodes coupling of the ones from the wave equations in homogeneous media (see e.g. [8], [6]). Moreover, the randomness of the media interfaces complicates the coupling; due to the randomness, we have already observed the PC modes coupling as indicated in (8).
Here the randomness makes the matrix Λ(x) dense, i.e. the PC modes all interact each other for x in the ranges of random interfaces. The eigenmodes of the wave solutions of (1) in a homogeneous media, e.g. c =μ, are easily found. Setting ) . Taking into account the 2π-periodicity in x we find that ω/c 0 = j, j = 0, 1, · · · . Hence the wave solution u can be written u = ∞ j=0 v j e ij(c0t+x) . However, in an inhomogeneous media, the eigenmodes are coupled. Writing u = ∞ j=0 v j (t)e ijx and substituting for u in (1) we obtain that, for k = 0, 1, which is supplemented with the initial conditions derived from (1) 3,4 . This equation is a system of ordinary differential equations and well-posed. The random inhomogeneities make all eigenmodes coupled with each other.
2.3. Stochastic eigenvalue problems. From the wave equations (1) setting u = e iωt v(x) we can obtain a non-deterministic Helmholtz equation: If the variance of noise is small, we expect to find a good frequency ω which might be close to the one from the deterministic problem (24). Although the eigenvalues and eigenvectors for the Helmholtz equation (24) may not be found analytically, discretizing the problem by Fourier spectral method we obtain a linear system (see [21]) where A is the spectral differentiation matrix for −∂ 2 /∂x 2 , B = diag(1/µ 2 (x j )) and v = (v(x j )). The eigenvalues ω 2 and corresponding eigenvectors v can be found using e.g. eig function in MATLAB. For the stochastic problem (38) we obtain a random matrix B = B ξ1,ξ2 = diag(1/c 2 (x j )) instead. We can now numerically compare the eigenvalues and eigenvectors for (39) and for (39) with B being replaced by B ξ1,ξ2 . We expect that the eigenvector corresponding to the smallest eigenvalue has a low frequency which makes robust with the change of the lengths l 1 , l 2 , l 3 , l 4 (i.e. the variances of ξ 1 , ξ 2 ), whereas the eigenvector with a large eigenvalue tends to be sensitive to the change of the lengths. For more discussion on stochastic eigenvalue problems (38) and (39), the readers are referred to [19], [28] respectively.

2.4.
Numerical examples. Figure 4 shows the effect of random fluctuation of interfaces to the wave solutions of (1). The mean solution u 00 of (8) is very close to the solution u * of the mean dynamics (35) when the variance or random fluctuation is small (first row) and the eigenvector v has a low frequency or small eigenvalue (a). More precisely, the mean u 00 of solution with high frequency (c) is much sensitive to even small random fluctuation whereas the one with low frequency is robust to it. From this we can infer that the uncertainty caused by the random interfaces cannot be ignored only except in the case when the random fluctuation is small and the solutions have low frequency. We have to take into account the uncertainty in the output solutions where we use the PC method (4). It satisfactory describes the output stochastic solutions by solving the modal system (8). Furthermore, the exponential decay of the PC modes as in (3) supports using the PC method. In Figure 5 via Monte Carlo simulations applied to the PC expansions (4) we plot the cumulative distribution function (CDF) for the solutions u(x, t) of Eq. (1) as in [13], [20]. If the variance of the random interfaces pertaining to the parameter c is small as in (a), we also observe small variance in the output u(x, t) whereas the large variance of them results in large variance of the output as in (b). If the input variance of random interfaces is small, we may ignore the output variance and compute the mean dynamics. However, if it is not small, solutions from the mean dynamics are very different from the mean output as indicated in Figures 4, 5. 3. Conclusions. We have constructed a numerical tool to investigate statistical quantities, e.g. CDF, mean and variance, of the wave solutions of Eq. (1) with randomly layered media, more precisely, media with random interfaces. The basic idea of the numerical tool is projecting the stochastic equations (1) on a random field via the so-called Legendre polynomial chaos and computing the PC modes by standard numerical tools. We compute and compare the mean solutions u = u 00 (x, t) repectively with N 1 = N 2 = 0, i.e. the mean dynamics (35), and with N 1 , N 2 sufficiently large. The difference between those two solutions indicates the robustness with respect to random interfaces. It turns out that the eigenmodes with a small[large] eigenvalue is robust[sensitive] with the random changes of media interfaces. This implies that a slow[fast] part of solutions is less[more] likely to be subjected to randomness. We have implemented an efficient numerical tool taking into account two random interfaces due to different media properties. More interfaces will be discussed in a future study. We can also consider inverse problems like finding an interface position based on the desired probability distributions of wave solutions (see e.g. [12]) and these problems will be studied in subsequent articles.