The interior resonance problem of time domain integral equations (TDIEs) formulated to analyze acoustic field interactions on penetrable objects is investigated. Two types of TDIEs are considered: The first equation, which is termed the time domain potential integral equation (TDPIE), suffers from the interior resonance problem, i.e., its solution is replete with spurious modes that are excited at the resonance frequencies of the acoustic cavity in the shape of the scatterer. Numerical experiments demonstrate that, unlike the frequency-domain integral equations, the amplitude of these modes in the time domain could be suppressed to a level that does not significantly affect the solution. This is achieved by increasing the numerical solution accuracy through the use of a higher-order discretization in space and the band limited approximate prolate spheroidal wave function with high interpolation accuracy as basis function in time. The second equation is obtained by linearly combining TDPIE with its normal derivative. The solution of this equation, which is termed the time domain combined potential integral equation (TDCPIE), does not involve any spurious interior resonance modes but it is not as accurate as the TDPIE solution at non-resonance frequencies. In addition, TDCPIE's discretization calls for treatment of hypersingular integrals.

## I. INTRODUCTION

Many applications in engineering and physical sciences call for simulations of acoustic scattering from penetrable objects, i.e., scatterers that internally support nonzero velocity potential and pressure field.^{1–11} An acoustic scattering problem involving penetrable objects is also known as an acoustic transmission problem. The time-harmonic (frequency-domain) acoustic transmission problem can be analyzed by solving a set of integral equations enforced on the surface of the scatterer.^{12} This set of equations is obtained by using the Kirchhoff-Helmholtz theorem to express the scattered fields of the exterior and interior problems in terms of (unknown) velocity potential on the surface of the scatterer and its normal derivative. The exterior and interior problems involve the unbounded domains with the material properties (density and wave speed) of the background medium and the scatterer, respectively. Numerical schemes developed to solve these integral equations discretize only the surface of the scatterer and implicitly enforce the radiation condition at infinity,^{12} offering advantages over finite element and finite difference methods that directly solve the Helmholtz equation, and require a volumetric discretization of the whole computation domain and use absorbing boundary conditions to approximate the radiation condition.

On the other hand, traditional integral equation formulations suffer from so-called “interior resonance” problem.^{13–21} This problem is observed when the excitation frequency approaches any one of the resonance frequencies of the acoustic cavity that is in the shape of the scatterer and has the density and the wave speed of the background medium. At these frequencies, the surface integral operator has a null space and the corresponding equation does not have a unique solution. Several approaches have been proposed to address the interior resonance problem of the frequency-domain integral equations. Examples of these include the combined Helmholtz integral equation formulation^{13} and the Burton-Miller scheme.^{14}

Even though interior resonance problem is well-studied for the frequency-domain integral equations, there are only a couple of studies that investigate the spurious interior resonance modes in the solution of time domain integral equations (TDIEs) of acoustics.^{22–25} In Ref. 22, a spurious resonance-free Burton-Miller-type time domain combined field integral equation is formulated to analyze acoustic scattering from sound-rigid bodies. In Ref. 26, interior resonance modes observed in the solution of the time domain electric field integral equation (of electromagnetics) that is enforced on perfect electrically conducting scatterers are investigated. Theoretically, TDIEs should not admit any interior resonance modes since their solution is obtained under zero initial condition and the interior resonance modes do not satisfy this initial condition.^{26,27} However, the interior resonance modes are still observed in the time domain solutions. It is discussed in Ref. 26 that this is because of the numerical errors introduced due to discretization and matrix inversions carried out during time marching.

In this work, the interior resonance problem of two different TDIEs formulated to analyze the transient acoustic transmission problem is investigated: The first equation, which is termed time domain potential integral equation (TDPIE) (in unknowns velocity potential and its normal derivative) here, is the time-domain equivalent of the frequency-domain integral equation that is traditionally used in the literature to solve the acoustic transmission problem.^{6,9,10,28} This equation suffers from the interior resonance problem. Its solution is replete with spurious modes that oscillate (without any decay) with the resonance frequencies of the acoustic cavity that is in the shape of the scatterer and has the density and the wave speed of the background medium.^{15} These modes are excited when the band of excitation includes their resonance frequency. In this work, it is demonstrated that unlike the frequency-domain integral equations, the amplitude of these modes in the time domain could be suppressed to a level, which does not significantly affect the solution, by increasing the accuracy of the numerical solution. This is achieved by using a higher-order Nyström discretization^{29–33} in space and the band limited approximate prolate spheroidal wave (APSW) function^{34–36} with high interpolation accuracy as basis function in time. The APSW function is a “two-sided” interpolator, which means that the resulting TDIE solver is non-causal. To address this and enable time-marching, a highly-accurate extrapolation scheme^{37,38} is used to “restore” the causality.

On the other hand, the second equation investigated in this work, which is termed time domain combined potential integral equation (TDCPIE), completely eliminates the interior resonance problem. The frequency-domain counterpart of TDCPIE has been introduced in Ref. 8 and it has been theoretically shown that its frequency-domain solution is unique at the interior resonance frequencies. Numerical results in this work verify that TDCPIE does not admit any interior resonance modes. It should be noted here that TDCPIE is obtained by linearly combining TDPIE with its normal derivative. Coupling parameters of this combination are carefully selected to enable the computation of the strongly-singular and hypersingular integrals that appear in the expressions of the matrix elements resulting from the Nyström discretization in space.^{29–33} In addition to this more-complicated singularity treatment method, the TDCPIE solution is, in general, less accurate than the TDPIE solution (except at the interior resonance frequencies of TDPIE).

Note that a preliminary version of this work has been described in a conference contribution.^{39} This paper includes a significantly more detailed formulation of TDPIE and TDCPIE and their discretization as well as a systematic investigation of the effects of solution accuracy on the spurious resonance modes associated with TDPIE. The remainder of the paper is organized as follows. In Sec. II, TDPIE and TDCPIE are derived. Section III describes the spatial and temporal discretization schemes and the marching-on-in-time method that is used to solve the resulting matrix system. In Sec. IV, numerical results are presented to validate the accuracy of TDPIE and TDCPIE solutions and demonstrate the relationship between the numerical errors and the interior resonance modes observed in the solution of the TDPIE. Section V concludes the paper with a short summary.

## II. FORMULATION

Let $\Omega 2$ and $\Omega 1$ denote the supports of an acoustically penetrable scatterer and the unbounded background medium in which this scatterer resides, respectively (Fig. 1). Let *S* represent the surface that separates these two domains, i.e., the surface of the scatterer. The wave speed and the density in $\Omega k,\u2009k\u2208{1,2}$, are $ck$ and $\rho k$, respectively. An acoustic field with velocity potential $\phi i(r,t)$ is incident on *S*. It is assumed that $\phi i(r,t)$ is band limited to maximum frequency $fmax$ and vanishingly small for $t\u22640$ on $r\u2208S$. In response to this excitation, scattered fields with velocity potentials $\phi ks(r,t)$ are generated in $\Omega k$. Total velocity potentials in $\Omega 1$ and $\Omega 2$ are expressed as $\phi 1(r,t)=\phi i(r,t)+\phi 1s(r,t)$ and $\phi 2(r,t)=\phi 2s(r,t)$, respectively.

Using the Kirchhoff-Helmholtz theorem,^{40} $\u2202t\phi 1(r,t)$ and $\u2202t\phi 2(r,t)$ are expressed as^{28}

Here, $\u2202t$ denotes the temporal derivative, $\u2202n=n\u0302(r)\xb7\u2207,n\u0302(r)$ is the outward pointing unit normal at point $r\u2208S$, and the spatiotemporal integral operators $Sk[x](r,t)$ and $Dk[x](r,t)$ are given by

where “$*$” denotes temporal convolution and

is the time domain Green function. Note that since the Green function is in the form of a Dirac delta function $\delta (t\u2212t0)$, the temporal convolutions in operators $Sk[x](r,t)$ and $Dk[x](r,t)$ reduce to retarded time integrals.

On *S*, the acoustic pressure field and the normal component of the velocity field are continuous, i.e., the velocity potential satisfies the following boundary conditions:^{28}

### A. TDPIE

Taking the limit of Eqs. (1) and (2) as $r$ approaches *S* from $\Omega k,\u2009k\u2208{1,2}$ and inserting Eqs. (3) and (4) into the resulting equations yield TDPIE as^{28}

Here, “∼” on top of $D\u0303k$ means that the space integral in $Dk$ is evaluated in the principal value sense.^{41}

Note that the time-derivative form of TDIEs is used in this paper because this form of TDPIE is the one originally formulated in Ref. 28 to analyze the transient acoustic transmission problem. Ideally, one could remove the time derivative and still solve the resulting equation. However, as explained in Refs. 42–45, for the solution of a TDIE to be stable, basis and testing functions should be selected from appropriate Sobolev spaces. More specifically, the basis function should reside in the domain space of the integral operator and the testing function should reside in the space dual to the range of the integral operator. Following the analysis carried out in Refs. 42 and 43, the testing function for TDPIE without the time derivative should be selected as the time-derivative of the Dirac delta function. This means that a marching method, which relies on point-testing in time (i.e., testing with Dirac delta function) as described in Sec. III of this paper, does not yield a stable solution when it is used to solve TDPIE without the time derivative.

### B. TDCPIE

where the spatiotemporal integral operators $D\u2032k[x](r,t)$ and $Nk[x](r,t)$ are given by

and $\u2202nn\u20322=\u2202n\u2202n\u2032$ denotes the double normal derivative. Taking the limit of Eqs. (7) and (8) as $r$ approaches *S* from $\Omega k,\u2009k\u2208{1,2}$ and inserting Eqs. (3) and (4) into the resulting equations yield the normal derivative of TDPIE in Eqs. (5) and (6) as

Linearly combining Eqs. (5) and (6) and (9) and (10) as *α*_{1} [Eq. (5)] +$\alpha 2\rho 2/\rho 1$ [Eq. (6)] and $\beta 1$[Eq. (9)] +$\beta 2$[Eq. (10)] yields TDCPIE as

## III. NUMERICAL SOLUTION

To solve Eqs. (5) and (6) and (11) and (12) numerically, first, *S* is discretized into a mesh of curvilinear triangles and surface unknowns $\phi 1(r,t)$ and $\u2202n\phi 1(r,t)$ are expanded in space and time as

In Eqs. (13) and (14), $Nt$ is the number of time steps, $Np$ is the number of curvilinear triangles, $Nn$ is the number of interpolation nodes on each triangle, $\u2113qn(r)$ is the Lagrange interpolation function defined at $rqn$ (node *n* on triangle *q*),^{29}^{,}$\u03d1(r)$ is the inverse of the Jacobian of the coordinate transformation between the unit right triangle and the Cartesian coordinate system, *T*(*t*) is the temporal basis function, which is constructed using the APSW function,^{34–36} $\Delta t$ is the time step size, and $Ii1|qn$ and $Ii2|qn$ are the unknown expansion coefficients to be solved for.

Substituting Eqs. (13) and (14) into Eqs. (5) and (6) and (11) and (12) and point-testing the resulting equations in space at $rpm,\u2009p=1,\u2026,Np,\u2009m=1,\u2026,Nn$ (i.e., Nyström discretization in space), and in time at $t=j\Delta t$ yield the following system of matrix equations:

Here, $Nhw$ is the half-width of the APSW function *T*(*t*), $Ij=[Ij1\u2009\u2009Ij2]T$, where $Ij1$ and $Ij2$ store the unknown expansion coefficients in Eqs. (13) and (14), $Vj=[Vj1\u2009\u2009Vj2]T$ store the tested excitation vectors at time $t=j\Delta t$, and

store the discretized retarded time integrals between nodes of mesh elements. Expressions of elements of $Vj$ and $Zj\u2212i$ for TDPIE and TDCPIE are provided in Secs. III A and III B, respectively. Note that the system of matrix equations in Eq. (15) is not causal, i.e., $Ij$ can not be solved for without knowing “future” unknowns $Ij+1,\u2009Ij+2$,…, $Ij+Nhw$ [see the second summation on the right-hand side of Eq. (15)]. An extrapolation scheme^{34,37,38} can be used to represent these future unknowns $Ij+1,\u2009Ij+2$,…, $Ij+Nhw$ in terms of “past/current” unknowns $Ij\u2212N+1$,…, $Ij\u22121,\u2009Ij$:

Here, *N* is the number of samples used in the extrapolation and $pn$ stores the extrapolation coefficients. Inserting this expression into Eq. (15) converts it into a causal form as

where the modified matrices $Z\xafj\u2212i$ are expressed as

The extrapolation coefficients $pn$ are computed using the method described in Refs. 37 and 38. One can also use the scheme described in Ref. 34, but this comes with a slight decrease in accuracy.

The system of matrix equations in Eq. (16) is now in a form that can be recursively solved for the unknown coefficient vectors $Ij,\u2009j=1,\u2026,Nt$ *via* time marching as briefly described next. For *j =* 1, $I1$ is found by solving Eq. (16) with the right-hand side $V1$ (contribution from the summation is zero at the first time step). For *j =* 2, the right-hand side of Eq. (16) is computed by subtracting $Z\xaf1I1$ (only term coming from the summation) from $V2$. Then, $I2$ is found by solving Eq. (16) with this right-hand side. For $j=3,$ the right-hand side of Eq. (16) is computed by subtracting $Z\xaf2I1\u2009+\u2009Z\xaf1I2$ from $V3$. Then, Eq. (16) is solved for $I3$. This recursive time marching algorithm is continued until all $Ij,\u2009j=1,\u2026,Nt$ are obtained.

The marching-on-in-time method described above, thanks to the higher-order discretization in space and the use of the APSW function as basis function in time, yields a highly-accurate numerical solution. Indeed, results presented in Sec. IV show that, compared to the traditional marching-on-in-time methods, which uses polynomial temporal basis functions, it significantly suppresses the amplitude of the spurious non-decaying oscillations due to the interior resonance modes.

### A. Elements of V_{j} and Z_{j−i} for TDPIE

The elements of $Vj1$ and $Vj2$ are given as $Vj1|pm=\u2202t\phi i(rpm,t)|t=j\Delta t$ and $Vj2|pm=0$, respectively. The elements of $Zj\u2212i$ are expressed as

Here, $R=|rpm\u2212r\u2032|,\u2009Sq$ is the surface of the curvilinear triangle *q*, and $\delta pq=1$ for *p *=* q*, and $\delta pq=0$ for $p\u2260q$.

### B. Elements of V_{j} and Z_{j−i} for TDCPIE

The elements of $Vj1$ and $Vj2$ are given as $Vj1|pm=\alpha 1\u2202t\phi i(rpm,t)|t=j\Delta t$ and $Vj2|pm=\beta 1\u2202t\u2202n\phi i(rpm,t)|t=j\Delta t$, respectively. The elements of $Zj\u2212i$ are expressed as

### C. Computation of singular integrals

There are four kinds of integrands in Eqs. (17) and (18): $Gk(R,t)*T(t\u2212i\Delta t),\u2009\u2202nGk(R,t)*T(t\u2212i\Delta t),\u2009\u2202n\u2032Gk(R,t)*T(t\u2212i\Delta t)$, and $\u2202nn\u20322Gk(R,t)*T(t\u2212i\Delta t)$. When *p *=* q*, and as $r\u2032$ approaches $rpm$, the integrals of $Gk(R,t)*T(t\u2212i\Delta t)$ and $\u2202nGk(R,t)*T(t\u2212i\Delta t)$ become weakly-singular, the integral of $\u2202n\u2032Gk(R,t)*T(t\u2212i\Delta t)$ becomes strongly-singular, and the integral of $\u2202nn\u20322Gk(R,t)*T(t\u2212i\Delta t)$ becomes hypersingular.^{46} The weakly-singular integrals are computed using the Duffy transformation.^{47} The strongly-singular integrals in Eq. (17) are computed using the approach described in Ref. 30. The strongly-singular and hypersingular integrals in Eq. (18) are treated by extending the singularity subtraction and cancellation methods^{48–51} to account for the temporal discretization, more specifically the temporal basis function. This extension is briefly described next.

$T(t\u2212R/ck)$ is expanded using the Taylor series around *R=0* as

Using this expansion in $\u2202n\u2032{T(t\u2212R/ck)/(4\pi R)}$ and $\u2202nn\u20322{T(t\u2212R/ck)/(4\pi R)}$ yields

Here, $R\u0302=(rpm\u2212r\u2032)/R$ and $O(Ru\u22650)$ represents the higher-order terms that are not singular as $R\u21920$. The singular terms on the right hands of Eqs. (20) and (21) are subtracted from $\u2202n\u2032Gk(R,t)*T(t\u2212i\Delta t)$ and $\u2202nn\u20322Gk(R,t)*T(t\u2212i\Delta t)$, respectively. Then, the integrals of these terms are added back to yield the final expressions for $Zj\u2212i11|pm,qn$ and $Zj\u2212i21|pm,qn$ in Eq. (18) as

The first and the third integrals on the right-hand side of Eq. (22) and the first and the fourth integrals on the right-hand side of Eq. (23) are the integrals of the integrands remaining after the singularity subtraction. Since these integrands are “smooth” (i.e., not singular as $R\u21920$), their integrals can be computed accurately using a Gaussian quadrature rule.^{52} The third and the sixth integrals on the right-hand side of Eq. (23) are subtracted weakly-singular integrals and they are computed using the Duffy transformation.^{47} The second and the fourth integrals on the right-hand side of Eq. (22) are subtracted strongly-singular integrals and they cancel out each other for $\alpha 1=\alpha 2$. Similarly, the second and the fifth integrals on the right-hand side of Eq. (23) are subtracted hypersingular integrals and they cancel out each other for $\beta 1\rho 2=\beta 2\rho 1$.

## IV. NUMERICAL RESULTS

In this section, numerical results, which demonstrate the relationship between numerical errors and interior resonance modes, are presented. TDPIE and TDCPIE are used to analyze acoustic scattering from a penetrable unit sphere that resides in an unbounded backgroud medium. It is assumed that the sphere is centered at the origin. The wave speed in the background medium and inside the sphere is $c1=300$ and $c2=200\u2009m/s$, respectively. The ratio of the densities in these two media is $\rho 1/\rho 2=1.5$. In all simulations, the excitation is a plane wave with velocity potential

where $\phi 0$ is the amplitude, $k\u0302i$ is the unit vector along the direction of propagation, and $Gmod(t)=cos\u2009[2\pi f0(t\u2212tp)]\u2009exp\u2009[\u2212(t\u2212tp)2/(2\sigma 2)]$ is a modulated Gaussian pulse with center frequency $f0$, time delay $tp$, and duration *σ*. The excitation parameters are selected as $\phi 0=1\u2009m2/s,\u2009k\u0302i=\u2009z\u0302,\u2009f0=120\u2009Hz,\u2009tp=10\sigma $, and $\sigma =3/(2\pi fbw)$, where the effective bandwidth $fbw=80\u2009Hz$. Note that this definition of *σ* ensures that $99.997%$ of the wave energy is within the frequency band $[fmin,fmax]$ with $fmin=f0\u2212fbw$ and $fmax=f0+fbw$.^{53} Also, this specific selection of $f0=120\u2009Hz$ and $fbw=80\u2009Hz$ ensures that the frequency of the lowest interior resonance mode (the first cavity mode of the Dirichlet problem), $150\u2009Hz$, is within the frequency band $[fmin,fmax]$, i.e., this resonance mode could possibly be excited using the Gaussian pulse described above.^{15} The time step size is chosen as $\Delta t=1/(2\gamma fmax)$ with oversampling factor *γ*. Unless stated otherwise, the APSW function with an half-width of $Nhw=7$ is used to construct *T*(*t*). The surface of the sphere is discretized using $Np=396$ curvilinear triangles with $Nn=6$ interpolation nodes on each triangle. For TDCPIE, the linear combination coefficients in Eqs. (11) and (12) are chosen as $\alpha 1=\alpha 2=\beta 1=1$ and $\beta 2=2/3$. Lower-upper (LU) decomposition is used to solve the matrix system in Eq. (16) (at every time step) to ensure that the error in the matrix solution is at the machine precision level.^{54}

### A. Accuracy of TDPIE and TDCPIE

For the first set of simulations, the oversampling factor is selected as $\gamma =6$ resulting in $\Delta t=0.42\u2009ms$. Sixteen-point Gaussian and 9-point Gauss-Legendre quadrature rules^{52} are used to compute the two-dimensional (2D) surface integral with “smooth” integrand and the one-dimensional (1D) line integral needed for the Duffy transformation,^{47} respectively. First, three different TDPIE solutions are compared: The first solution is obtained using the proposed method where the temporal basis function *T*(*t*) is constructed using the APSW function. The second and the third solutions are obtained using the “traditional” marching-on-in-time method, where *T*(*t*) is constructed using third- and fourth-order Lagrange polynomial (LP) functions, respectively.^{55,56} It should be noted here that all solutions are obtained using the same discretization in space. Also, note that the traditional scheme with LP functions does not require an extrapolation scheme (as described in Sec. III) to ensure the causality of the time marching. Figure 2 plots the magnitude of the expansion coefficient $Ij1|qn$, *q =* 3, *n =* 3 (corresponding to $rqn=(0.78,0.58,0.23)\u2009m$), $j=1,\u2026,Nt$ computed by solving TDPIE using the APSW, the third-order LP (termed “LP3”), and the fourth-order LP (termed “LP4”) functions. The figure clearly shows that all three solutions are corrupted by non-decaying oscillations that are especially visible in the late time. These oscillations correspond to the spurious interior resonance mode with frequency 150 Hz. The figure also shows that the amplitude of this resonance mode is significantly smaller in the TDPIE solution obtained using the APSW function. This is observed because the APSW function leads to a more accurate solution than the LP functions and this observation supports the claim that spurious resonance modes are induced in the TDPIE solution because of the numerical errors. It should be noted here that the solution obtained using the APSW function is more accurate because it is band limited (with a short duration time) and has excellent interpolation properties and continuous derivatives at $t=j\Delta t$ while the LP functions have discontinuous derivatives and a wideband spectrum.^{34,35,37}

Next, the solutions of TDPIE and TDCPIE, both of which are obtained using the APSW function, are compared. Figure 2 plots the magnitude of the expansion coefficient $Ij1|qn$, *q =* 3, *n =* 3 (corresponding to $rqn=(0.78,0.58,0.23)\u2009m$), $j=1,\u2026,Nt$ computed by solving these two equations. Clearly, the solution of TDPIE is corrupted by non-decaying oscillations while the solution of TDCPIE is free from any resonances. Figures 3(a), 3(b), and 3(c) compare the normalized Fourier transform of $\phi 1(r,t)$ [i.e., Fourier transform of $\phi 1(r,t)$ divided by the Fourier transform of $Gmod(t)$], which is computed after solving TDPIE and TDCPIE, in the time range $t\u2208[0,1.67s],\u2009t\u2208[0,0.27s]$, and $t\u2208[0,0.17s]$ to the frequency-domain (time-harmonic) total velocity potential computed at $r=(0.78,0.58,0.23)\u2009m$ using the Mie series solution,^{57} respectively. Note that the normalization is required to ensure that the time-harmonic response (with equal excitation amplitude at each frequency) is obtained from the solutions of TDPIE and TDCPIE. Figure 3(a) clearly shows that both simulations are accurate within the effective band of the excitation, except at $150\u2009Hz$ where the solution of TDPIE is corrupted by the interior resonance mode. Figure 3(a) also shows that TDPIE is more accurate than TDCPIE at other frequencies. This might be explained by the fact that a second-kind surface integral equation (e.g., TDCPIE) is usually less accurate than its first-kind counterpart (e.g., TDPIE).^{58} The results in these three figures show that the visibility of the error due to the spurious resonance mode (in the solution of TDPIE) at 150 Hz decreases with the decreasing length of the Fourier transform window. More specifically, in Fig. 3(c), which presents the results with the shortest Fourier transform window, the error at 150 Hz cannot be distinguished from the numerical errors at other frequencies. This might be interpreted as that the amplitude of the spurious mode is already small enough not to be visible in the main scattering response. But when the Fourier transform window is increased to include the late time response, the presence of the spurious mode becomes more visible. Note that all three figures, as expected show that the solution of TDCPIE is free from any interior resonance modes.

The effect of interior resonances on the scattered velocity potential is investigated by comparing the scattering cross section (SCS) of the sphere computed using the normalized Fourier transformed solutions of TDPIE and TDCPIE to SCS computed using the Mie series solution.^{57} Figures 4(a) and 4(b) plot SCS versus *θ* for $\varphi =0\xb0$ at 120 and 150 Hz, respectively. As shown in Fig. 4(a), SCS computed using TDPIE and TDCPIE solutions at 120 Hz shows good agreement with the Mie results. On the other hand, as shown in Fig. 4(b), the resonance mode at $150\u2009Hz$ dramatically changes SCS computed using the TDPIE solution.

### B. Effect of numerical integration accuracy

In this set of simulations, the effect of the computation accuracy of the integrals in Eqs. (17) and (18) on the amplitude of the interior resonance modes is investigated. Three sets of computation accuracy are considered by using different number of quadrature points to compute the 2D surface integral with “smooth” integrand and the 1D line integral needed for the Duffy transformation:^{47}

16-point Gaussian and 9-point Gauss-Legendre quadrature rules

7-point Gaussian and 5-point Gauss-Legendre quadrature rules

4-point Gaussian and 3-point Gauss-Legendre quadrature rules

In all simulations, the oversampling factor is selected as $\gamma =6$ resulting in $\Delta t=0.42\u2009ms$. To clearly identify the interior resonance mode, Fourier transforms of the late-time data of $Ij1|qn$, *q =* 3, *n =* 3 [corresponding to $rqn=(0.78,0.58,0.23)\u2009m$] computed after solving TDCPIE and TDPIE with integration accuracy sets (i), (ii), and (iii) in the time range $t\u2208[0.8\u2009s,1\u2009s]$ are plotted in Fig. 5. The figure shows that the solutions of TDPIE with all three sets exhibit spurious interior resonance mode at $150\u2009Hz$. However, as expected, no spurious resonance mode is present in the solutions of TDCPIE. Furthermore, the amplitude of the interior resonance mode observed in the solution of TDPIE with set (iii) is stronger than those with sets (i) and (ii). This is because set (iii) yields larger numerical errors than sets (i) and (ii).

To compare the accuracy of the solutions of TDCPIE and TDPIE with sets (i), (ii), and (iii), the $L2$-norm error in the normalized Fourier transform of $\phi 1(r,t)$ (over all 2376 interpolation nodes on the sphere surface) and the $L2$-norm error in SCS (for $\theta =[0\xb0,180\xb0]$ with 361 sample points and $\varphi =0\xb0$) are plotted versus frequency in Figs. 6(a) and 6(b), respectively. Note that the reference data used in the computation of the $L2$-norm errors is obtained using the Mie series solution. Figure 6 shows that the accuracy of the TDPIE solution in the vicinity of $150\u2009Hz$ is significantly affected by the interior resonance mode. As expected, larger numerical errors increase the amplitude of the resonance mode. Figure 6 also shows that the spurious interior resonance mode is not observed in the TDCPIE solution regardless of the integral computation accuracy.

To visualize the interior resonance mode at $150\u2009Hz$, the normalized Fourier transform of $\phi 1(r,t)$ (at all interpolation nodes on the sphere surface) computed after solving TDCPIE and TDPIE with set (iii) are presented in Figs. 7(a) and 7(b), respectively. Figure 7(c) presents the difference in the normalized Fourier transform of $\phi 1(r,t)$ obtained from the TDCPIE and TDPIE solutions. As expected, the pattern of the difference follows the amplitude of the interior resonance mode.^{59}

### C. Effect of time step size

In this set of simulations, the effect of time step size $\Delta t$ (i.e., temporal discretization density) on the amplitude of the interior resonance modes is demonstrated. Three different oversampling factors are used: $\gamma \u2208{6,7.5,10}$ resulting in $\Delta t={0.42,0.33,0.25}\u2009ms$, respectively. In all simulations, the integration accuracy set (i) described in Sec. IV B is used to ensure that the numerical error resulting from the computation of the space integrals in matrix elements is suppressed to be sufficiently small. Figure 8 compares the magnitude of the expansion coefficient $Ij1|qn$, *q =* 3, *n =* 3 [corresponding to $rqn=(0.78,0.58,0.23)\u2009m$], $j=1,\u2026,Nt$ computed by solving TDPIE with $\gamma =6,\u2009\gamma =7.5$, and $\gamma =10$. The figure shows that the amplitude of interior resonance mode as identified by the non-decaying oscillations in the late time decreases by reducing the time step size. Figure 9 plots the Fourier transforms of $Ij1|qn$, *q =* 3, *n =* 3 [corresponding to $rqn=(0.78,0.58,0.23)\u2009m$] computed after solving TDCPIE and TDPIE with $\gamma =6,\u2009\gamma =7.5$, and $\gamma =10$ in the time range $t\u2208[0.8\u2009s,1\u2009s]$. As expected, interior resonance mode is observed in the TDPIE solution at $150\u2009Hz$. Furthermore, the figure shows that using a larger *γ* (or smaller $\Delta t$) reduces the amplitude of the resonance mode. As expected, no interior resonance mode is observed in the TDCPIE solution regardless of the value of *γ* used. Figures 10(a) and 10(b) plot the $L2$-norm error in the normalized Fourier transform of $\phi 1(r,t)$ (over all 2376 interpolation nodes on the sphere surface) and the $L2$-norm error in SCS (for $\theta =[0\xb0,180\xb0]$ at 361 sample points and $\varphi =0\xb0$) computed after solving TDCPIE and TDPIE with $\gamma =6,\u2009\gamma =7.5$, and $\gamma =10$, respectively. Note that the reference data used in the computation of the $L2$-norm errors is obtained using the Mie series solution. Figure 10 shows that interior resonance mode is observed in all of the TDPIE solutions with different *γ* but its amplitude could be significantly suppressed by increasing *γ*. Another point to note here is that, even though the TDCPIE solution does not admit any interior resonance modes, it is usually less accurate than the TDPIE solution within the whole effective band of the excitation except in the vicinity of $150\u2009Hz$.

## V. CONCLUSION

The interior resonance problem of TDPIE and TDCPIE that are formulated to analyze the time domain acoustic field interactions on penetrable scatterers is investigated. Numerical results demonstrate that the solution of TDPIE is corrupted by the spurious interior resonance modes that oscillate (without any decay) with the resonance frequencies of the acoustic cavity in the shape of the scatterer and has the density and the wave speed of the background medium. However, unlike the frequency-domain integral equations, the amplitude of these modes in the time domain can be suppressed by reducing the numerical error. This is achieved by using a higher-order Nyström discretization in space and the band limited APSW function with high interpolation accuracy as basis function in time.

The solution of TDCPIE, which is obtained by linearly combining TDPIE with its normal derivative, is free from spurious interior resonance modes but in general is less accurate than the solution of TDPIE (except at the interior resonance frequencies of TDPIE). In addition, the discretization of TDCPIE calls for treatment of strongly-singular and hypersingular integrals.

## ACKNOWLEDGMENTS

This publication is based upon work supported by the King Abdullah University of Science and Technology (KAUST) Office of Sponsored Research (OSR) under Award No. 2019-CRG8-4056. The authors would like to thank the King Abdullah University of Science and Technology Supercomputing Laboratory (KSL) for providing the required computational resources.