In general I work at the intersection of numerical analysis and mathematical physics. In much of my work physical insight has allowed radical improvements in numerical methods; in turn, the development of efficient numerical methods has allowed a deeper investigation and understanding of physical phenomena, particularly in quantum systems. I am also working on a medical imaging inverse problem, where due to the complexity of the problem, many levels of approximation are necessary and physical insight is essential. For this noisy real-world application, a statistical understanding of the imaging process is important---this contrasts much existing work in the field of inverse problems.

Also here are some recent (up to Nov 2006) pages I have put up on certain special topics:

- Global basis sets for Helmholtz boundary value problem
- Quantum mushroom billiards
- Random plane waves

Finite Element or finite-difference type methods are impractical when
there is a large number of wavelengths across the system, therefore
boundary methods are necessary.
I have been developing variants of
the `Method of Particular Solutions' (MPS), which can be viewed as
approximating trial eigenfunctions by a linear sum of basis functions
which satisfy the PDE (the Helmholtz equation at some trial
energy parameter *E*) in the domain.
Approximations to the eigenvalues are then found by searching for *E*
values where the minimum achievable
*boundary norm* becomes very small.
The MPS was largely abandoned by the numerical analysis community in the
1980's, however many exciting independent developments have occured in
the physics (quantum chaos) community since then, using similar ideas,
upon which I base my work.

In my thesis work I placed the domain norm and boundary norm on an
equal footing, thus tackling the problem that there can exist nontrivial
basis combinations which are exponentially small everywhere in the domain.
I also derived formulae allowing the domain norm to be evaluated using
boundary integrals alone, improving efficiency.
I also gave the first account of the working of the powerful
`scaling method' of
Vergini, which can be seen as an accelerated MPS, whose relative
efficiency grows like *O(k)* where *k=E ^{1/2}*.
The scaling method
produces

I am now analysing the intrinsic errors of eigenfunction and eigenvalues
computed by the scaling method, using the same tool.
Despite its use within the quantum chaos community, such errors
have not been yet understood.
In both the scaling method, and the improved inclusion bounds, the
role of the matrix *Q _{ij}* of
inner products of the eigenfunction
normal derivative functions, under the boundary weight function

Basis sets for the linear space of Helmholtz solutions, necessary for the MPS and scaling methods, are not well understood. In the original MPS work regular Bessel functions were used, and the physics community has based much work on sets of plane waves. I have developed basis sets of Neumann functions (fundamental solutions) placed on a curve exterior to the domain, which appear to represent all eigenfunctions accurately, even for nonconvex domains with corners. I am analysing approximation properties of such bases, in the semiclassical limit. Can bases for re-entrant corners be developed? Are there shapes for which convential BIE will be the only option?

The Quantum Unique Ergodicity (QUE) conjecture of Rudnick-Sarnak, made
in the context of negatively-curved manifolds, is that
there is *no exceptional subset* in the QET.
There is much recent interest in QUE within pure mathematics (number theory and
ergodic theory).
QUE has recently been proven by Lindenstrauss for automorphic forms,
where the infinite number of symmetries allows an analytic approach.
In billiards no such analytic machinery is known, however
I have found strong numerical evidence for QUE in a uniformly-hyperbolic
billiard system.

To what extent do eigenfunctions behave like random sums of plane waves
(the prediction of a Random Matrix Theory [RMT] model)?
Many open issues remain in this area.
Correlations between eigenfunctions cannot be predicted by RMT,
since they involve the dynamics of the system.
For instance, the
variance of off-diagonal matrix elements *A _{ij}* is
proven to converge to a classical autocorrelation function.
I have tested this convergence and again found it suprisingly slow.
The prediction for diagonal variance is related by a factor

There are several applications of this work.
In quantum physics, such matrix elements form the backbone of study of
driven systems and transition rates (Fermi Golden Rule).
Off-diagonal matrix element variance corresponds to *dissipation*
(heating) rate under periodic driving.
In `quantum dots' (cold confined 2D electron gases) chaotic effects
have been instrumental in understanding fluctuations in electrical
conduction.
Quantum dots are now model systems for study of spins ('spintronics')
and are promising candidates for quantum computers.

** Forward problem**:
DOT requires numerical simulation of the forward problem,
that is, the passage of photons through a medium of complex
geometry, in our case the human head.
We approximate the time-dependent transport equation by the
diffusion equation, since it can be solved much more rapidly.
In tests against Monte Carlo simulations
we have found diffusion to be adequate.
In iterative or Bayesian approaches to inverse problems,
the forward model must be called a large number of times during the
optimization (fitting process)
or exploration of the posterior probability density function.
The scalp/skull/brain system is approximately a layered medium,
for which we were surprised to find a lack of efficient diffusion solvers.
I therefore have developed an efficient layered solver
which makes use of the symmetry, a 1D finite-difference method in the *z*
direction, and a logarithmic timestepping scheme
motivated by the physical diffusion mode decay, giving ~0.1 sec CPU run time
in relevant regimes.
There are many pressing goals such as applying this timestepping scheme
to more general 3D finite-difference discretizations.
Since MRI data is cartesian, and triangulation in 3D is a tricky problem,
we expect this to be superior to currently-used Finite Element schemes.
Discretization errors due to poor representation of the smooth head
surface by a cartesian grid should be overcome by the inclusion of
immersed interface methods into the finite-difference scheme.
The *gradient* (Jacobean) of model outputs with respect to inputs
is vital to the optimization procedure, therefore I intend to apply adjoint
differentiation techniques to these forward models.

**Inverse problem**:
To make extraction of useful brain activity tractable,
our key approach, proposed by David Boas and his group in Boston,
is to use structural MRI to define the subject's head geometry, then use
some kind of segmented optical model with a reduced number of
parameters (of order 10 or 100 unknowns) in order to fit the optical signals.
Use can be then made of information about each tissue type:
optical parameters of the skull are fixed, and the scalp
oxygenation correlates strongly with `global' physiological parameters.
Thus changes in the brain can be isolated with more certainty.
I have performed a simulated feasibility study of this segmented approach,
fitting *baseline* brain parameters, a nontrivial problem in itself.
This used a Bayesian approach, enabling the full probability density
function (PDF) in the parameter space to be estimated. This gives
marginal PDFs on each parameter, and correlations. Ultimately in an
imaging problem with many unknowns, low resolution,
and noisy (single photon counting) data collection,
I argue that such statistical approaches will be the best way forward,
despite the higher computational cost.

Currently in my work PDFs are explored using Markov chain Monte Carlo (MCMC) methods, such as Metropolis algorithms, however there is much to be gained by testing hybrid MCMC methods and combining with gradient-based optimization. In the limit of small noise, the posterior PDF peak approaches Gaussian, and we have established what range of imaging parameters this approximation is adequate. With a Gaussian PDF approximation, expensive MCMC can be replaced with much cheaper optimation followed by estimation of the Hessian (inverse of covariance) matrix. I have implemented this method on layered medium fits, with various numbers of layers, where the layer thicknesses themselves are unknowns. Questions such as whether estimation of skull/scalp layer thickness by DOT or via MRI give more accurate DOT optical parameter fits, and the effect of introduction of (unknown) calibration parameters on the fit are vital to answer in order to advance the reliability of the technology.