## Numerical Method
### Evaluation of [A][x] Pre-Corrected Fast Fourier Transform Method  
An alternative approach is to use the pre-corrected fast Fourier transform (PFFT) accelerated algorithm, which can evaluate the product \([A][x]\) with \(O(NlogN)\) operations without explicit construction of the matrix \([A]\). The PFFT-BEM, which reduces the computational effort to \(O(NlogN)\), has been applied in various fields, such as the design of interaction circuits, the modeling of microelectromechanical systems, wave hydrodynamics problems, and acoustic scattering problems (e.g., [9; 10; 11]). Here, we briefly introduce the PFFT technique. For illustration, we take a radiation sphere as an example in each iteration step of GMRES. As shown in Figure 2a, using the boundary element method, we discretize the surface of the sphere into \(N\) small triangular elements with surface sizes of \(S_{j},j=1\),..., \(N\). We assume \(p_{n}\) is unknown on the sphere surface and we try to evaluate  
\[[A][x]=\sum_{j=1}^{N_{i}}\left(p_{n}\right)_{j}\iint\limits_{E_{j}}G(\mathbf{x} \mathbf{\xi})dS_{j}(\mathbf{\xi}). \tag{19}\]  
In the step of grid definition, we define \(N_{g}=N_{x}\times N_{y}\times N_{z}\) grids in the 3D domain covering the entire surface of the sphere, as shown in Figure 2b. We refer to the cube made up by nearby grids as a cell. Based on a previous study using PFFT [10], the size of the cells \(h\) and the elements \(\Delta l\) are assumed to be the same. In the step of projection, the normal gradient of the pressure \(p_{n}\) at element \(j\) is projected to the grids of the cell that contain element \(j\). The formula for projection is  
\[q_{jm}=\left(p_{n}\right)_{j}\iint\limits_{E_{j}}H_{m}(\mathbf{\xi})d\mathbf{\xi}(\bm {\xi}), \tag{20}\]  
where \(H_{m}\) is the spatial interpolation function (see, e.g., [10]) for the \(m^{th}\) grid of the cell covering element \(j\), and \(q_{jm}\) is the strength of the projected \(\left(p_{n}\right)_{j}\) at the \(m^{th}\) grid of the cell covering element. Figure 2: Demonstration of the grid definition, projection, convolution, and interpolation of the pre-corrected fast Fourier transform (PFFT)-boundary element method (BEM) numerical scheme. (**a**) boundary element mesh; (**b**) fast Fourier transform (FFT) grid definition over the element mesh; (**c**) demonstration of (1) projection, (2) convolution, and (3) interpolation. In the step of convolution, we evaluate the total influence at grid point \(l\) due to all the projected \(\left(p_{n}\right)_{j}\) at other grids as  
\[\left(IF_{s}\right)_{lm}=\sum_{m=0}^{N_{S}-1}q_{m}G_{l-m},\ \ \ l=0,1,\ldots,N_{S}-1, \tag{21}\]  
where \(G_{l-m}=\frac{1}{\left[Y_{l}-\xi_{m}\right]}e^{ik\left(Y_{l}-\xi_{m}\right]}\). We note that the summation in Equation (21) is in convolution form and can be computed using FFT with a computation operation of O(\(N_{S}\)\(log\)\(N_{S}\)) and a memory requirement of O(\(N_{S}\)). \[\left(IF_{s}\right)_{lm}=FFT^{-1}(FFT(q_{m})\cdot FFT(G_{l})),\ l=0,1,\ldots,N _{S}-1. \tag{22}\]  
In the step of interpolation, we evaluate the total influence \(\left(IF_{s}\right)_{ij}\) at element \(i\) (\(i=1,\ldots,N\)) due to all \(\left(p_{n}\right)_{j}\) using the interpolation function \(H_{m}\) and the grid values \(\left(IF_{s}\right)_{lm}\). The steps of projection convolution and interpolation are demonstrated in Figure 2c as steps (1)-(3). It is clear that the accuracy of the interpolation of the influence between elements \(j\) and \(i\) deteriorates when \(\frac{\left[N-\xi\right]}{h}=O(1)\) or smaller, where \(h\) is the size of the cells.