Computational high frequency waves in heterogeneous media

1. Computational high frequency waves in heterogeneous media Shi Jin University of Wisconsin-Madison

2. High frequency waves • Fig. 1. The electromagnetic spectrum, which encompasses the visible region of light, extends from gamma rays with wave lengths of one hundredth of a nanometer to radio waves with wave lengths of one meter or greater. • High frequency waves: wave length/domain of computation <<1

3. Difficulty of high frequecy wave computation • Consider the example of visible lights in this lecture room: wave length: » 10-6 m computation domain » m 1d computation: 106» 107 2d computation: 1012» 1014 3d computation: 1018» 1021 do not forget time! Time steps: 106» 107

4. An Example: Linear Schrodinger Equation

5. If one can find the asymptotic (semiclassical) limit as  -> 0 then one can just solve the limiting equation numerically. This allows coarse (underresolved) computation Semiclassical limit of the linear schrodinger equation

6. Semiclassical limit in the phase space Wigner Transform A convenient tool to study the semiclassical limit (Lions-Paul; Gerard, Markowich, Mauser, Poupaud)

7. Moments of the Wigner function The connection between W and  is established through the moments

8. The semiclassical limit

9. Other high frequency waves • The Wigner method is generic: linear wave—geometrical optics elastic waves electromagnetic waves, etc • We always end up at Liouville equation or coupled Liouville systems Ryzhik-Papanicoulou-Keller • Most recent numerical methods for high frequency waves are based on Liouville equations

10. Liouville equation ft + r H¢rx f - rx H ¢r f = 0 f(t, x, ) is the density distribution of a classical particle at position x, time t, with momentum  H=H(x, ) is the Hamiltonian The bicharacterisitcs of this equation is a Hamiltonian system: dx/dt = r H d/dt = -rx H Classical mechanics: H=1/2 ||2+V(x) (=> Newton’s second law) Geometrical optics: H = c(x) || The Liouville equation can be solved by method of characteristicsif H is smooth

11. Discontinuous Hamiltonians • H=1/2||2+V(x): V(x) is discontinuous-potential barrier, • H=c(x)||: c(x) is discontinuous-different index of refraction • quantum tunneling effect, semiconductor devise modeling, plasmas, geometric optics, interfaces between different materials, etc. • The semiclassical limit breaks down at barrier/interface

12. Analytic issues ft + r H¢rx f - rx H ¢r f = 0 • The PDE does not make sense for discontinuous H. What is a weak solution? dx/dt = r H d/dt = -rx H • How to define a solution of systems of ODEs when the RHS is discontinuous or/and measure-valued?

13. Numerical issues • for H=1/2||2+V(x) • since V’(x)= 1 at a discontinuity of V, one can smooth out V then Dv_i=O(1/x), thus poor resoultion (for complete transmission) wrong solution (for partial transmission) t=O( x )

14. Mathematical and Numerical Approaches Q: what happens before we take the high frequency limit?

15. Snell-Decartes Law of refraction • When a plane wave hits the interface, the angles of incident and transmitted waves satisfy (n=c0/c)

16. Our Idea: An interface condition • We introduce an interface condition for f that connects (the good) Liouville equations on both sides of the interface. • T, R defined from the original “microscopic” problems • This gives a mathematically well-posed problem that is physically relavant • We can show the interface condition is equivalent to Snell’s law in geometrical optics • A new method of characteristics (bifurcate at interfaces) f(x+, +)=Tf(x-,-)+R f(x+, -+) for +>0 H(x+, +)=H(x-,-) R: reflection rate T: transmission rate R+T=1

17. Solution to Hamiltonian System with discontinuous Hamiltonians T R • Particles cross over or be reflected by the corresponding transmission or reflection coefficients (probability) • Based on this definition we have also developed particle methods (both deterministic and Monte Carlo) methods

18. Key idea in Hamiltonian-preserving schemes • consider a standard finite difference approximation V: piecewise linear approximation—allow good CFL fI,j+1/2, f-i+1/2,j ---- upwind discretization f+i+1/2, j ---- incorporating the interface condition (motivated byBerthame & Simeoni: kinetic scheme for shallow water equations with bottom topography)

19. Scheme I (finite difference formulation) • If at xi+1/2 V is continuous, then f+i+1/2,j= f-i+1/2,j; • Otherwise, For j>0, f+i+1/2,j = f(x+i+1/2, +) = T f-(x-i+1/2, -) +R f(x+i+1/2, -+)

20. The transmitted and reflected fluxes • if the particle is transmitted 2) If the particle is reflected:

21. The New CFL condition • Note the discrete derivative of V is defined only on continuous points of V, thus t=O( x, )

22. Positivity, stability, l1-convergence • for first order scheme (forward Euler in time + upwind in space), under the “good” CFL condition if fn >0, then fn+1 > 0; k fn+1kl1 (x, )·k fnkl1 (x, ) k fnk1 ·Ck f0k1 (except for measure-valued initial data) l1-convergence

23. Curved interface

24. Geometrical optics • The same idea has also been extended to geometric optics H = c(x) || with partial transmission and reflection We build in Snell’s Law into the flux References: J-Wen, semiclassical limit of Schrodiger, Comm Math Sci ’05 J-Wen, geometrical optics, JCP 06, SINUM

25. Quantum barrier

26. A semiclassical approach for thin barriers(with Kyle Novak, SIAM MMS, JCP) • Barrier width in the order of De Broglie length, separated by order one distance • Solve a time-independent Schrodinger equation for the local barrier/well to determine the scattering data • Solve the classical liouville equation elsewhere, using the scattering data at the interface • Ben Abdallah-Degond-Gamba (1d quantum-classical coupling) • Ben Abdallah-Tang (1d, stationary computation)

27. A step potential ( V(x)=1/2 H(x) )

28. Resonant tunnelling

29. Circular barrier (Schrodinger with =1/400)

30. Circular barrier (semiclassical model)

31. Circular barrier (classical model)

32. Diffraction grating:

33. Semiclassical

34. Semiclasical vs Schrodinger (=1/800)

35. Computation of diffraction (with Dongsheng Yin)

36. Transmissions, reflections and diffractions (Type A interface)

37. Type B interface

38. Hamiltonian preserving+Geometric Theory of Diffraction • We uncorporate Keller’s GTD theory into the interface condition:

40. A type B interface

41. Another type B interface

42. A type A interface

43. Half plane

44. Other applications/developments • Elastic waves (with X. Liao, JHDE) • High frequency waves in random media (with X. Liao and X. Yang, SISC) • Fast phase-flow particle method for the Liouville equations (with H. Wu and Z. Huang)

45. Computational cost (=10-6) • Full simulation of original problem for  x » t » O()=O(10-6) Dimension total cost 2d, O(1018) 3d O(1024) • Liouville based solver for diffraction  x » t » O(1/3) = O(10-2) Dimension total cost 2d, O(1010) 3d O(1014) Can be much less with local mesh refinement

46. Summary • Unique, physically relevant solution to linear transport equation with discontinuous and measure-valued solution • Probability solution to Hamiltonian system with discontinuous Hamiltonians • Finite-difference, finite-volume methods with interface condition built into the numerical flux • Particle (both deterministic and Monte-Carlo) methods • Able to compute (partial) transmission, reflection, and diffraction (by building geometrical theory of diffraction into the interface condition) for many high frequency waves (geometrical optics, semiclassical limit of Schrodinger, elastic wave, thin quantum barrier, etc.): Liouville equation + interface conditions