Nonlocal operators such as layer and volume potential operators are prevalent in physical sciences. We are looking towards making such evaluators accessible from UFL. As an example application, we consider a wave-structure interaction model for photoacoustic trace gas sensors. The sensors of interest are based on quartz-enhanced photoacoustic spectroscopy or resonant optothermoacoustic detection. The model consists of thermoacoustic waves in the exterior fluid domain and thermoelastic waves in the interior solid domain. We plan to use integral equation method (IEM) for the exterior and FEM for the interior. In this contribution, we propose a novel second-kind integral equation formulation for the exterior solution that solves the Morse–Ingard equation subject to Neumann boundary conditions on the solid surface. The resulting boundary integral equation is solved with GMRES where quadrature-by-expansion with fast-multipole acceleration is used to evaluate the nonlocal integral operators. Since the solution representation naturally satisfies the far-field conditions, domain truncation and perfectly matched layer are not needed. We demonstrate with examples that our solver has \(O(n)\) complexity, is high-order, and can handle complex geometries. For the next steps, we are working towards integrating our Morse–Ingard solver with FEM, and enabling more general IEM-FEM coupling in the process.