Research
My areas of research can be broadly categorized as multi-scale numerical methods for hyperbolic PDEs. More specifically, I work on efficient numerical algorithms for electromagnetic wave propagation. My doctoral dissertation addressed fractional relaxation models (FRMs) for electromagnetic (EM) pulse propagation through dispersive dielectric media, which required numerical implementation of fractional derivatives. My postdoctoral work has focused on the development of an implicit Maxwell solver, to address the massive scale separation that arises in plasma simulations. Below is a summary of these two projects.
More recently, I have been dabbling into the numerical computation of special functions. I recently proved that an original algorithm due to Lanczos for the Gamma function is interpolatory. This led to a deeper insight to a paper by Nakatsukasa and Trefethen, where the AAA algorithm was used. I am also looking into other related functions, such as the multiple gamma and hyperfactorial functions.
Causley, M. The Gamma Function via Interpolation. Numerical Algorithms, vol. xx, no. x, pp. x-xx, 2021.
Nakatsukasa, Y., Sète, O., and Trefethen, L. N. The AAA Algorithm for Rational Approximation. SIAM J. Sci. Comput., vol 40, no. 3, pp. A1494-A1522, 2018.
Implicit Maxwell Solver
Plasma, one of the four fundamental states of matter, arises when kinetic energy is increased inside of a gas, causing it to ionize. Understanding plasmas better is relevant to many areas of our daily lives as well as active fields of research, ranging from the behavior of the sun and other astral bodies, to energy applications, such as plasma assisted combustion.
From the standpoint of simulation, plasma simulations scale very badly in space, due to the large number of particles required to resolve the plasma dynamics. Usually explicit solvers, such as the finite-difference time-domain (FDTD) method, are used to solve Maxwell's equations, which obey the Courant-Friedrichs-Lewy (CFL) time step restriction. However, due to the inclusion of the plasma particles, this time step restriction forces explicit solvers to run for months or even years to produce an accurate simulation!
However, implicit solvers for Maxwell's equations remove this time step restriction, and can simulate plasmas in a reasonable amount of clock time. First Maxwell's equations, are reduced to solving the wave equation, using the Lorenz gauge. Next, using the method of lines transpose (MOLT), the system is reduced to a coupled set of (elliptic) boundary value problems. Using boundary integral methods, the solution is fully discretized. Our method has three advantages over explicit Maxwell solvers:
1) Large time steps can be taken, to "step over" the smaller time scales of less relevant physics.
2) Point sources can be located anywhere (even off the mesh), and be tracked accurately.
3) Complex geometries can be addressed by embedding them within Cartesian meshes.
However, there are still challenges. For instance, taking large time steps can result in a loss of accuracy. To mitigate this difficulty, we have constructed higher order, A-stable implicit solvers. Using a novel technique, successive convolution in space is used to achieve higher temporal accuracy, without sacrificing stability properties, so that large time steps can be taken, without sacrificing accuracy! Furthermore, this method of approach is quite general, and holds great promise for efficiently solving parabolic problems, such as nonlinear reaction-diffusion systems.
Below are some movies of our solver in action, testing the implementation of point sources, periodic and outflow boundary conditions.
An implicit solution demonstrating point sources.
The spatial mesh is taken to be dx = 0.005 in both directions, and the time step is dt = 0.01. Outflow boundary conditions are prescribed at the right, Dirichlet boundary conditions on the left, and the top and bottom edges are periodic.
A cosine bump is specified in each side of a circular cavity,
which contains cusps at the intersection of the two circles. The spatial mesh is taken to be dx = 0.0006... in both directions, and the time step is dt = 0.012... Dirichlet boundary conditions are prescribed along the boundary of the cusp.
An initial Gaussian is pulse is released in the center of a domain, and periodically traverses up and to the right. For illustrative purposes, 9 identical copies of the domain are plotted, to enhance visual effect.
References
Causley, M.F., and Christlieb, A. Higher order A-stable schemes for the wave equation using a recursive convolution approach. SIAM J. Numer. Anal., vol. 52, no. 1, pp. 220–235, 2014.
Causley, M.F., Christlieb, A., and Wolf, E. Method of Lines Transpose: A Fast Implicit Wave Propagator. submitted, Mathematics of Computation, 2013.
Causley, M.F., Christlieb, A., Van Groningen, L. and Ong, B. Method of Lines Transpose: An Implicit Solution to the One Dimensional Wave Equation. to appear, Mathematics of Computation, 2013.
Dielectric Fractional Relaxation Models (FRMs)
My dissertation focused on the numerical aspects of modeling pulse propagation through dispersive dielectric media that exhibit fractional relaxation. Basically, I designed algorithms for incorporating fractional differential equations (FDE) into the finite-difference time-domain (FDTD) framework.
Fractional differential operators are defined via convolution with singular, slowly decaying kernels. This slow decay is used to describe "memory-dependent" features in physical systems. In the case of FRMs, the dielectric response (the frequency-dependent permittivity) persists for longer times, changing the phase and amplitude of the polarization field.
FDEs are both interesting and challenging from a numerical standpoint. A naive approach to convolution leads to an O(N2) algorithm since at each time step an integral of length N must be computed. This makes the problem intractable. Instead, fast convolution is used, which decomposes the kernel into a finite-rank sum of separable kernels. The error introduced in approximating the true kernel is known a priori, and the computational complexity is now reduced to N log(N). The other issue with FDEs accurate treatment of the singularity, which remains an area of active research.
Below is a movie depicting a fast rise-time sinusoidal electric pulse, propagating through an FRM described by the Havriliak-Negami model, showing how the fractional relaxation rapidly dampens high frequency content in the wave fronts, but then supports characteristic bumps, which decay algebraically, and contain heavy tails.
References
Causley, M.F. and Petropoulos, P.G. On Numerically Solving the Cole-Cole Dielectric Model. in preparation.
Causley, M.F. and Petropoulos, P.G. On the Time-Domain Response of Havriliak-Negami Dielectrics. Transactions on Antennas and Propagation, vol. 61, no. 6, pp. 3182-3189, 2013.
Causley, M.F., Petropoulos, P.G. and Jiang, S. Incorporating the Havriliak-Negami Dielectric Model in the FD-TD Method. Journal of Computational Physics, vol. 230, no. 10, pp. 3884-3899, 2011.