An FEM-MLMC algorithm for a moving shutter diffraction in time stochastic model

We consider a moving shutter and non-deterministic generalization of the diffraction in time model introduced by Moshinsky several decades ago to study a class of quantum transients. We first develop a moving-mesh finite element method (FEM) to simulate the determisitic version of the model. We then apply the FEM and multilevel Monte Carlo (MLMC) algorithm to the stochastic moving-domain model for simulation of approximate statistical moments of the density profile of the stochastic transients.


(Communicated by Eric Kaufmann)
Abstract. We consider a moving shutter and non-deterministic generalization of the diffraction in time model introduced by Moshinsky several decades ago to study a class of quantum transients. We first develop a moving-mesh finite element method (FEM) to simulate the determisitic version of the model. We then apply the FEM and multilevel Monte Carlo (MLMC) algorithm to the stochastic moving-domain model for simulation of approximate statistical moments of the density profile of the stochastic transients.
1. Introduction. In this article we describe and demonstrate a computational model, using a moving-mesh finite element method (FEM) and the multilevel Monte Carlo (MLMC) algorithm, to simulate statistical moments of a moving-domain induced stochastic quantum transients. For a review on the theory and applications of transients, we refer to [2] and 288 references therein. The earliest of these references are due to Moshinsky. A closely related deterministic stationary-domain version of the model considered in this article was proposed in the seminal 1952 paper by Moshinsky [10], entitled Diffraction in Time (DiT). Even several decades after introduction of the DiT model, there continues to be substantial interest in understanding more practical variants and solvability of the model. For a recent effort on an exactly solvable DiT model, see [6] and references therein. We refer to the above references for large literature on the DiT models, including reasons for the term DiT attributed to a quantum matter transient waves phenomenon related to that in optics.
Transients are processes induced by sudden interaction of particles impinging on a shutter or changes due to initial states or boundary conditions. As described in the review article [2,Page 2], the transient models investigated mainly in the (Physics) literature are those for which analytical solutions can be constructed using classical functions (such as Fresnel integrals). Such analytic constructions impose restrictions Figure 1. An illustration of the DiT model governed by a moving shutter with speed γ (t) and a continuous initial state. The shutter on the left represents the closed shutter's position at t = 0 and the right shutter is a closed shutter at the position determined by γ(t).
The deterministic case, with a constant γ(t) (and a discontinuous initial state) of the stochastic model considered in this article reduces to the stationary shutter DiT model introduced by Moshinsky [10], which models the deterministic quantum dynamics of a suddenly released beam of independent particles of mass m. In this special stationary case, the time evolution of a quantum particle with wave function ψ(x, t) satisfying the Schrödinger equation can be written as the Moshinsky function, M (x, k, t) [2,10], and the concept of DiT can be tied to the quantum density profile |ψ(x, t)| 2 [2,10]. Here, k is the wavenumber and λ = 2π/k is the wavelength, so that /λ and /(mλ) (with being the Planck's constant) are respectively the nominal momentum and velocity, of a quasi-monochromatic atomic beam impinging on a stationary shutter. These constants occur together as a coefficient of the diffusion part of the Schrödinger equation [10]. In this article we use the notation α(k) to denote the coefficient.

1.1.
A moving-domain stochastic DiT model. The main focus of this article is developing and demonstrating an algorithm to compute statistical moments of a quantity of interest (QoI), determined by the stochastic quantum density profile |ψ(x, t; ω)| 2 , which are crucial for quantifying uncertainties in the QoI. The stochastic quantum density profile integrand is determined by the evolution of the unknown wave function ψ satisfying the following stochastic Schrödinger model on a moving domain D(t) = (β, γ(t)), where γ(t) represents the position of the moving shutter at time t: In the above model, for convenience, we assume homogeneous Dirichlet boundary conditions at x = 0 and at x = γ(t; ω), for all t ≥ 0. In this article we take the non-constant speed of the shutter to be non-deterministic, of the form γ (t) = ωf (t), where ω is from a one-dimensional probability space (Ω, F, P) and f : [0, T ] → R is a real-valued function. We denote the resulting moving non-deterministic domain as D(t; ω) := (β, γ(t; ω)), for t ∈ [0, T ], with T being the final simulation time of interest. We note that for the domain to be well defined, γ(t; ω) ≥ β for all t ∈ [0, T ] and ω ∈ Ω. In practice at time t = 0, depending on the choice of f (t), the initial domain D(0; ω) need not depend on ω.
In addition, we augment the complexity of the variant stochastic version of the Moshinsky's problem by taking the initial state to be a random field, depending on a random vector ω. Here, ω is an element of Ω ⊆ R d , and (Ω, F , P) is the associated probability space with Euclidian sample space Ω and Borel σ-algebra F . The probabilities of the events in F are given by the probability measure P. For notational convenience, in (1)-(2) and throughout the article, we use ω = (ω, ω) to denote a random vector in Ω = Ω × Ω ⊂ R d+1 . The associated product probability measure in Ω is P = P × P.
Because of the complexity of the stochastic model on the moving domain, in this article, we take the initial random state of the model in (2) to be continuous in the spatial variable. Developing an efficient moving-mesh and stochastic realization algorithm of this continuous initial state model is a challenging problem. In a future article, we plan to investigate the extension of the model to allow for random discontinuous initial states on moving domains. Such discontinuous data based stochastic models are challenging even for fixed spatial domain cases.
In this setting, with a large d-stochastic dimension, the released initial state of the beam of particles is a square-integrable random process and the covariance operator of the random process induces the associated Karhunen-Loéve (KL) series expansion of the field [9]. In practice, assuming decay in the Fourier coefficients of the field, the initial state may be approximated by a truncated KL expansion with d-terms. This leads to our d + 1-dimensional stochastic model induced by the perturbation of the wavenumber dependent deterministic initial state µ(x)g k (x), x ∈ [β, γ(0)], by a random field R(x; ω). In the stationary shutter deterministic Moshinsky model on (−∞, ∞) from [10], g k (x) = exp(ikx) (a plane wave) and µ(x) = 0 for x ∈ [0, ∞). In our simulations, we also choose g k to be the plane wave and choose µ such that the initial state satisfies the boundary conditions at x = β and at x = γ(0).
To illustrate the decay assumption in the random field, say with mean zero, consider for example with x ∈ [0, 1], the eigenfunctions of the second differential operator −d 2 /dx 2 (with domain restricted to functions satisfying homogeneous Dirichlet boundary conditions at x = 0 and at x = 1) given by sin(jπx), and associated eigenvalues j 2 π 2 , for j = 1, 2, . . . . The integral covariance operator inducing the random field, C = [−d 2 /dx 2 ] −2 , has the same eigenfunctions, but the eigenvalues of C decay with O(j −4 ) [3]. Hence, depending on the appropriate accuracy required in the representation of the random field, we can choose a parameter d to truncate the KL expansion and represent the random field using d random variables ω j , for j = 1, . . . , d (say uniformly distributed on [−0.5, 0.5]) as For example with d = 64 (which we use in our numerical experiments), the truncated terms in the series KL expansion are less than O(10 −4 ). In our model, since the domain at time t = 0 is [β, γ(0)], with β < 0, the arguments of the sin(j·) functions need to be appropriately scaled using change of variables, to satisfy the homogeneous Dirichlet conditions at x = β and at x = γ(0). That is, we use the representation (4) In this article, we are interested in efficiently computing statistical moments of the QoI, Q( ω) at the final simulation time T , where For computation of an approximation to the first statistical moment (expected value) of the QoI in (5), we need to approximate the d + 1 dimensional integral An approximation of the high-dimensional integral in (6) requires a large number, say N M C , of sampling/realization points in the (d + 1)-dimensions and use of non-standard quadratures, such as the low-order accurate Monte Carlo (MC) [5] approximation. The advantage of the MC approximation is in its simplicity of choosing the realization points in the probability space ( Ω, F , P) and that the rate of convergence of the MC approximation of the integral in (6) is independent of the high stochastic dimension. However, because the accuracy of the MC approxima- in general a large number of samples, N M C , are required. For each such fixed sample, say ω * , we need to simulate the associated deterministic moving-domain DiT model to compute Q( ω * ) in (5) In the next section, we describe a moving-mesh FEM for efficiently solving the deterministic moving-domain DiT model to compute ψ * (x, t) and demonstrate the second-order accuracy of the algorithm. The Crank-Nicolson algorithm is unconditionally stable, requiring no CFL condition on the time discretion step, and in our simulation we take the spatial and time discretization steps to be of the same size. Our simulation of the deterministic FEM model suggests that a fine time-dependent spatial-mesh, say h f ine (t), of the moving domain is required, in conjunction with the matching fine time-mesh, to obtain better than O(10 −3 ) accuracy. In particular, we aim for this accuracy to be better than that which is used for the approximate representation of the random field R, using the truncated KL expansion. Hence, computing approximate expected values of the QoI using the standard MC algorithm requires a large number of realizations of the deterministic model, with a fine space-time mesh for each realization.
In the next section, we tabulate deterministic moving-mesh FEM simulation results for various mesh sizes, with (approximate second-order) increase in accuracy as the mesh becomes finer. These results also provide a framework to apply a multilevel version of MC that facilitates an appropriate trade-off between a large number of MC samples to a fewer number of MC samples in conjunction, respectively, with coarse and fine mesh, depending on the level. The initial level of MLMC [4] may be considered as standard MC, but with coarser space-time mesh and the later levels may be regarded as MC with fewer samples and a finer mesh.
In Section 3, we describe and apply the standard MC and MLMC algorithms to simulate the expected value of the QoI and demonstrate that the computational cost of the MLMC algorithm is substantially better than the MC algorithm, for the stochastic moving shutter DiT model. Throughout this article, even for a fixed time t, we use the operator notation ∂ ∂t acting on a time-variable dependent function to denote the first partial derivative of the function w.r.t. the second variable. Thus, for example, ∂ ∂t ψ * (x, t) means that we first take the partial derivative of the function ψ * w.r.t. the second variable and then evaluate the resulting function at the fixed time t.

2.
A deterministic moving-mesh computational DiT model. High-order finite difference and element methods for the Schrödinger equation have been applied in many areas over several decades, see for example [12,11,7]. However, these developed methods were mainly restricted to the Schrödinger equation on fixed domains. For recent work on deterministic moving-mesh models, see [8] and references therein, including the 2009 survey article [1] (with a large number of references) on the development of moving-mesh methods over the last four decades. To our knowledge, this article is the first work on applying moving-mesh techniques for simulation of the deterministic and stochastic DiT models.
In this section, for a fixed realization ω * = (ω * , ω * ), we describe an algorithm to efficiently approximate ψ * (x, t), satisfying the deterministic moving-domain Schrödinger model, induced by a monochromatic beam of particles moving parallel to the x-axis impinging on a shutter that moves perpendicular to the beam in its frame, with a deterministic speed γ (t) = ω * f (t): For each fixed t ∈ [0, T ] and v 1 , v 2 ∈ L 2 (D * (t)), we use the notation to denote the standard inner product of square integrable functions defined on D * (t). Using the homogeneous Dirichlet boundary conditions, it is easy to show that a weak form of the deterministic model (7)- (9) is: For For , spanned by continuous splines of degree p ≥ 1 defined on a spatial mesh of D * (t), say with uniform size h(t), satisfying the approximation theory property that for With this abstract setting, a theoretical way of defining a semi-discrete Galerkin spline FEM is to find ψ * h(t) (·, t) ∈ V h(t) satisfying (10) and (11), with ψ * , H 1 0 (D * (t)), and µg k + R(·, ω * ) replaced, respectively, with ψ * h(t) , V h(t) and P h(t) (µg k + R(·, ω * )). Here, P h(t) (µg k + R(·, ω * )) is a projection of the initial data µg k + R(·, ω * ) onto Given that our model problem requires a moving domain, there is a need to efficiently setup the time-dependent non-standard stiffness, mass, and related hybrid matrices. In particular, unlike the fixed-domain case, the linear spline functions depend nonlinearly on the time-mesh and hence efficient setting up of the movingmesh FEM algebraic system, including an appropriate Jacobian of the spline basis functions, is required.
To this end, below we give concrete details of a fully-discrete moving-mesh FEM, which includes setting up the system matrices and the resulting algebraic system that needs to be solved at each discrete time step. We restrict our attention below to the continuous linear splines (p = 1) and the approach can be generalized for high-order splines and non-uniform mesh parameters.
and for j = 2, . . . , N h , the basis functions are defined as In addition, we define φ 1 and φ N h +1 appropriately, depending on the boundary conditions. For m = 0, . . . , M ∆t , with h m = h(t m ), to compute ψ * hm (·, t m ; x m ) := ψ * hm (·, t m ) ∈ V h(tm) , we use the ansatz: In (17), using the initial condition (14), we obtain its projection onto V h(t0) by taking For m = 1, . . . , M ∆t , we obtain the N h − 1 unknown coefficients Ψ i (t m ) by forcing the representation in (17) to satisfy a fully-discrete version of (13) by a second-order discretization of ∂ ∂t using the mid-time-nodal points t m−1/2 = (t m + t m−1 )/2 with mesh-width ∆t/2. Because of the nonlinear dependence of the basis functions φ i in the time variable, it is convenient to first consider a simplified representation of ∂ ∂t ψ * hm (·, t m ). Below we use the notation x m j to denote the j-th component of the vector x m , for j = 1, . . . , N h + 1. Using (17) we then obtain where for the bracketed terms, in the penultimate step we used (17) and the final step can be obtained by calculating derivatives of the basis functions in (16), w.r.t. the components in x m and relating to the derivative w.r.t. the spatial variable x. Thus, using (19), The above simplified representation yields, in addition to the standard fixed-mesh mass matrix, terms involving discretization of moving mesh terms and an associated matrix with entries composed of both the basis functions and their derivatives.

A fully-discrete moving mesh FEM and complex algebraic systems.
In this section, we use the following simplified notation to approximate the time derivatives in (20) at t m−1/2 = (m − 1/2)∆t, for m = 1, . . . , M ∆t with second order accuracy obtained by letting Following a Crank-Nicolson type approach, for i = 1, . . . , N h + 1, using (16), let φ Using the anstaz (17), we consider the approximation, for x ∈ D * (t m−1/2 ), Using (20), (21) and (22), our fully-discrete moving-mesh FEM approximation to (10) can be described as follows. Using the data Ψ 0 i ∈ C in (18), for m = 1, . . . , M ∆t , i = 1, . . . , N h + 1, find the coefficients Ψ m i ∈ C of the ansatz (17) of the Using (22) and expanding further, we obtain N h ] T , and consider the following (N h − 1) × (N h − 1) mass, stiffness, and hybrid matrices, with (i, j)-th entries being Using these matrices in (23), the moving-mesh FEM algebraic system, for each m = 1, . . . , M ∆t , can be written as Using the initial dataΨ 0 , we can use direct complex algebraic system solvers for (25), to compute the unknown coefficients in (17) and hence simulate the approximate wave function ψ * hm (·, t m ; x m ) = ψ * hm (·, t m ) ∈ V h(tm) , for m = 1, . . . , M ∆t . 2.2. Numerical experiments: Deterministic simulation. In this section, using the KL expansion in (4) with d = 64, we demonstrate the accuracy of the moving-mesh FEM model to compute ψ * hm (x, t m ; x m ) in (17), for 0 < t ≤ 2 and for one realization sample ω * = (ω * , ω * ) ∈ Ω ⊂ R 65 with each of the 64 components of ω * being a uniformly distributed random variable in [−0.5, 0.5] and ω * is the non-deterministic speed of the shutter, which is a uniformly distributed random variable in [1,2]. Using the resulting deterministic setting, we approximate and simulate the FEM discrete version of the deterministic continuous model (7)-(9) on [β, γ(t)] induced by a moving shutter with speed γ (t) = 2ω * t and the wavenumber dependent initial state induced at β = −2. Thus, the initial computational domain at t = 0 is [−2, 0] and grows until the computational domain is [−2, 4ω * ], at the final time T = 2.
For the randomly perturbed initial state in (8), we choose wavenumber k = π so that g k (x) = exp(ikx) and α(k) = k 2 , and to impose the homogeneous Dirichlet boundary conditions at the boundaries of D * (0), we choose µ(x) = (x−β)(x−γ(0)). For and hence compute the maximum nodal error with respect to both the time and spatial grid: (27) We start with the coarse spatial mesh and double the mesh size until N h = 10 * 2 10 and using the associated Err max (N h ) in (27), we compute the estimated order of convergence (EOC) of the moving-mesh FEM algorithm. We note that the theoretical FEM convergence rates are usually w.r.t. the maximum mesh size, and not w.r.t. the degrees of freedom used in the computation. Given that we fix N h for all increasing time-dependent sizes of the domain (and hence change the maximum mesh size) we will compute the EOC with respect to the spatial degrees of freedom N h .
Results in Table 1 demonstrate the accuracy of our moving-mesh FEM computational DiT model for a fixed realization of the initial state and shutter speed in the 65-dimensional stochastic space. In Figure 3, we visualize the computed density profile of the simulated solution ψ of the deterministic Schrödinger model on the moving domain, for three discrete time steps.    (5) on ω are its expected value given in (6) and variance These give, respectively, measures of the mean and spread of the QoI. A standard tool for uncertainty quantification is the Monte Carlo method, which yields an approximation to the expected QoI where ω j for j = 1, . . . , N M C are independent samples of the random vector ω ∈ Ω.
For each sample ω j , the QoI Q( ω j ) can be computed by a deterministic movingmesh DiT FEM model, with a fixed configuration described by ω j . The Monte Carlo method also yields an approximation to the variance In particular, for the E MC [Q; N M C ] approximation of E [Q] in (6), convergence of the MC approximation is w.r.t. the variance error The disadvantage of the MC approximation error is in its slow O(N −1/2 M C ) convergence rate. Hence, depending on the variance of the QoI, a large N M C number of samples are required to obtain a few matching digits of accuracy as the approximate E MC [Q; N M C ] values are computed, using an increasing number of samples. In addition, the terms in Q( ω j ) in (29) need to be further approximated by, say Q ∆t,h∆t ( ω j ), using the moving-mesh FEM algorithm described in the last section, with mesh parameters M ∆t and N h . The approximation Q ∆t,h∆t ( ω j ) is then obtained by replacing ψ(x, T ; ω j ) in (5) with ψ ∆t,h∆t (x, T ; ω j ), for j = 1, . . . , N M C .
Based on the results in Table 1, to ensure that the moving-mesh FEM error is less than the MC error, we took the fine space-time grid with M f ine ∆t = N f ine h = 20480 and performed the MC simulation for various values of N M C as described in Table 2.
In particular, we computed approximations to E MC [Q; N M C ] as to obtain the results in Table 2  We note that the MC algorithm does not provide an approach to choose an accuracy , for which the algorithm adaptively selects various numbers of samples to achieve a similar accuracy. The values in Table 2 are based on the estimated theoretical rate of convergence, assuming the variance of the QoI (occurring in the order constant) is of moderate size so that MC approximations achieve accuracy.
We could achieve such a large number of MC based realizations results in Table 2 because of our message passing interface (MPI) based parallel code written in modern versions of Fortran (90+). We ran the MC-FEM parallel code using multiple cluster nodes with each node composed of dual octa-core Intel X5670 2.93GHz processors. In order to compare the actual computational cost of our naturally parallel MC simulation with a relatively cheaper MLMC algorithm (described in Section 3.1), we provide the total core computational cost (in seconds) in Figure 6, which reflects the cost of the results in Table 2 for various values of . Figure 6 also demonstrates substantial reduction in the cost by instead using the MLMC algorithm.
3.1. Multilevel Monte Carlo simulation. The MLMC algorithm was introduced less than two decades ago and we follow the framework in the survey article [4] to further develop our parallel MPI FEM-MC code, which will substantially reduce the computational cost of the MC simulation to compute approximate expected values of the QoI.
Under certain assumptions [4, Theorem 1], it is possible to choose desired accuracy values (such as = [N M C ] −1/2 in our MC simulations). This can be achieved by performing a few more, say, L-levels of MC simulations, with integrands having smaller variances, as we increase from level 0 to level L. The reduction in the computational cost can be achieved using certain trade-offs between finer and coarser realizations of the space-time grids.
In particular, at level = 0 the MLMC algorithm is based on simulations using the largest number of adaptively chosen N Hence corresponding to (6), we obtain the telescopic sum relation The key outcome of this identity is that while the variance of Q ( ) for = 0, . . . , L may not be small, the variance of the difference Q ( ) − Q ( −1) terms in (34) are expected to reduce in magnitude as the level increases. The expected reduction in variances follows from the control variate analysis [5], as explained in [4,Page 3]. In Figure 4 we demonstrate this (using a semi-logy plot and L = 4) for our moving-mesh FEM based QoI and for demonstration purposes, the variance values were computed using a large number of test samples. Because of the substantially MLMC samples used at the initial level. Results in Figure 4, with very low variances of the difference integrands, also suggests that L = 3 is sufficient for our MLMC simulation.
Using the tabulated values in Table 1 Unlike the MC algorithm, the MLMC approach allows the user to provide an expected accuracy value and based on the coarse to fine grid space-time mesh, such as that in (35) and certain initial number of MC samples at each level (with prescribed minimum and maximum number of levels), the MLMC framework [4] adaptively chooses the final number of samples N ( ) MLMC to approximately satisfy the required accuracy. In particular, for our simulation we chose four values that correspond to the MC simulation values in Table 2 and the associated dependent values of N ( ) MLMC are given in Figure 5. We conclude this article with results in Table 3 that demonstrate the power of our FEM-MLMC algorithm to compute approximate expected values of the QoI with several matching digits at lower computational cost as demonstrated in Figure 6. Based on our further MLMC simulations for other values of = 0.01 * 2 −j for j = 0, 1, 2, 3, 4, we observed an approximate value of E[Q] is 4.6276, as demonstrated in Table 3