  We develop a novel method of virtual sources to formulate boundary integral
equations for exterior wave propagation problems. However, by contrast to
classical boundary integral formulations, we displace the singularity of the
Green's function by a small distance $h>0$. As a result, the discretization can
be performed on the actual physical boundary with continuous kernels so that
any naive quadrature scheme can be used to approximate integral operators.
Using on-surface radiation conditions, we combine single- and double-layer
potential representations of the solution to arrive at a well-conditioned
system upon discretization. The virtual displacement parameter $h$ controls the
conditioning of the discrete system. We provide mathematical guidance to choose
$h$, in terms of the wavelength and mesh refinements, in order to strike a
balance between accuracy and stability. Proof-of-concept implementations are
presented, including piecewise linear and isogeometric element formulations in
two- and three-dimensional settings. We observe exceptionally well-behaved
spectra, and solve the corresponding systems using matrix-free GMRES
iterations. The results are compared to analytical solutions for canonical
problems. We conclude that the proposed method leads to accurate solutions and
good stability for a wide range of wavelengths and mesh refinements.
