FDTD: Solving 1+1D delay PDE in parallel

Abstract

We present a proof of concept for solving a 1+1D complex-valued, delay partial differential equation (PDE) that emerges in the study of waveguide quantum electrodynamics (QED) by adapting the finite-difference time-domain (FDTD) method. The delay term is spatially non-local, rendering conventional approaches such as the method of lines inapplicable. We show that by properly designing the grid and by supplying the (partial) exact solution as the boundary condition, the delay PDE can be numerically solved. In addition, we demonstrate that while the delay imposes strong data dependency, multi-thread parallelization can nevertheless be applied to such a problem. Our code provides a numerically exact solution to the time-dependent multi-photon scattering problem in waveguide QED. Program summaryProgram Title: FDTD: solving 1+1D delay PDEProgram Files doi:http://dx.doi.org/10.17632/mmyw3fgjxh.1Licensing provisions: MITProgramming language: C (C99)Nature of problem: This program solves an unconventional 1+1D delay PDE that emerges in the study of waveguide quantum electrodynamics. The delay PDE is complex-valued and has a non-local delay term, and the solution to it provides the full dynamics of the system consisting of a few 1D photons and a two-level system in front of a mirror.Solution method: The finite-difference time-domain (FDTD) method is adapted. Given the initial condition of the system, the corresponding boundary condition is generated, and then the FDTD solver marches through the entire spacetime grid. Multiple solvers are supported using either OpenMP (wavefront) or pthreads (swarm).Additional comments including restrictions and unusual features: 1. Depending on the input parameters the memory and disk usages of the program can be excessive, so the users should choose the parameters wisely (see main text). 2. The multi-thread support using OpenMP is turned on by default. See README for how to turn it off and switch to pthreads instead. 3. As a by-product, a numerical routine is provided for evaluating the incomplete Gamma function γ(n,z) with nonzero positive integers n≥1 and complex-valued z.