  This article presents novel numerical algorithms based on pseudodifferential
operators ($\Psi$DO) for fast, direct solution of the Helmholtz equation in
one-, two- and three-dimensional inhomogeneous unbounded media. The proposed
approach relies on an Operator Fourier Transform (OFT) representation of
$\Psi$DO which frame the problem of computing the inverse Helmholtz operator,
with a spatially-dependent wave speed, in terms of two sequential applications
of an inverse square root $\Psi$DO. The OFT representation of the action of the
square root $\Psi$DO, in turn, can be effected as a superposition of solutions
of a pseudo-temporal initial-boundary-value problem for a paraxial equation.
The OFT framework offers several advantages over traditional direct and
iterative approaches for the solution of the Helmholtz equation. The operator
integral transform is amenable to standard quadrature methods and the required
pseudo-temporal paraxial equation solutions can be obtained using any suitable
numerical method. A specialized quadrature is derived to evaluate the OFT
efficiently and an alternating direction implicit method, used in conjunction
with standard finite differences, is used to solve the requisite component
paraxial equation problems. Numerical studies, in 1, 2, and 3 spatial
dimensions, are presented to confirm the expected OFT-based Helmholtz solver
convergence rate. In addition, the efficiency and versatility of our proposed
approach is demonstrated by tackling nontrivial wave propagation problems,
including 2D plane wave scattering from a geometrically complex inhomogeneity,
3D scattering from turbulent channel flow and plane wave transmission through a
spherically-symmetric gradient-index acoustic lens. All computations, even the
latter lens problem which involves solving the Helmholtz equation with more
than one billion complex unknowns, are performed in a single workstation.
