An explicit marching-on-in-time scheme for solving the time domain Kirchhoff integral equation

A fully explicit marching-on-in-time (MOT) scheme for solving the time domain Kirchhoff (surface) integral equation to analyze transient acoustic scattering from rigid objects is presented. A higher-order Nyström method and a PE(CE)m-type ordinary differential equation integrator are used for spatial discretization and time marching, respectively. The resulting MOT scheme uses the same time step size as its implicit counterpart (which also uses Nyström method in space) without sacrificing from the accuracy and stability of the solution. Numerical results demonstrate the accuracy, efficiency, and applicability of the proposed explicit MOT solver.


I. INTRODUCTION
Transient acoustic scattering from rigid objects can be numerically analyzed using either differential equation solvers [1][2][3][4] or time domain surface integral equation (TDSIE) solvers. [5][6][7][8][9][10][11][12][13][14] TDSIE solvers have several advantages over differential equation solvers. 4 (i) They require only a twodimensional discretization of the scatterer surface as opposed to a three-dimensional discretization of the whole computation domain. (ii) They are free from numerical dispersion since they do not discretize spatial/temporal derivatives using finite differences or finite elements. (iii) They do not need approximate absorbing boundary conditions to truncate the unbounded physical domain into the computation domain.
To construct a TDSIE, the scattered pressure field is expressed as a spatio-temporal convolution of the velocity potential generated on the surface of the scatterer (upon excitation) and the Green function of the unbounded background medium. Then, a TDSIE in (unknown) velocity potential is obtained by enforcing the boundary condition of the total pressure field (excitation plus the scattered pressure field) on the rigid surface of the scatterer. One of the classical schemes for solving TDSIEs is the marching-on-in-time (MOT) method. 5,[8][9][10][11]13,14 The MOT-TDSIE solvers expand the unknown velocity potential in terms of (local) spatial and temporal basis functions. Inserting this expansion into the TDSIE and testing the resulting equation in space and time yields system of equations. This system of equations is solved/updated for the unknown expansion coefficients through marching in time. The right-hand side of this system consists of the tested excitation and the scattered pressure field expressed as a discretized spatio-temporal convolution of the Green function and history of the velocity potential.
Despite the advantages of the TDSIE solvers listed in the first paragraph, use of the MOT-TDSIE solvers in analyzing transient acoustic scattering has been rather limited. This can be attributed to the difficulty of obtaining a stable solution as well as the high cost of computing the discretized spatio-temporal convolution. The former bottleneck has been alleviated with the development of new temporal basis functions 15 and testing schemes, 9,16,17 as well as accurate integration methods. [17][18][19][20] To reduce the high computational cost, plane wave time domain (PWTD) algorithm 11,12,21,22 and the fast Fourier transform (FFT)-based acceleration schemes [23][24][25][26] have been developed.
An MOT solver can be either explicit 13,14 or implicit 5-7 depending on the types of (spatial and temporal) basis functions and the testing scheme, as well as the time step size. The implicit MOT schemes are more stable and permit larger time steps without sacrificing from stability. But they call for solution of a matrix equation at every time step. On the other hand, the explicit MOT schemes do not require a matrix equation solution, however, they have to use a smaller time step size to ensure the stability. To overcome this drawback, a quasi-explicit MOT scheme has been developed for solving the time domain magnetic field integral equation (TD-MFIE) of electromagnetics. 27 This scheme employs Rao-Wilton-Glisson (RWG) functions 28 for spatial discretization and can use time step sizes as large as its implicit counterpart while maintaining the stability and the accuracy of the solution. However, it still needs solution of a sparse but unstructured (Gram) matrix equation at every time step because the RWG basis functions have overlapping supports/domains.
In this work, relying on the idea behind this quasiexplicit MOT scheme, a fully explicit MOT scheme is formulated and implemented to solve the time domain Kirchhoff (surface) integral equation (TDKIE) 10,29 for analysis of transient acoustic scattering from rigid objects. The TDKIE is cast in the form of a first-order ordinary differential equation (ODE) in time. The unknown velocity potential generated on the surface of the scatterer is expanded in terms of higher-order polynomials in space with time dependent (unknown) coefficients. Inserting this expansion into the TDKIE and point testing the resulting equation in space yield a system of ODEs (i.e., Nystr€ om method in space). This ODE system is then integrated in time using a PE(CE) m scheme for the (samples of) the unknown expansion coefficients. The resulting MOT scheme does not call for a matrix equation solution since the Gram matrix (resulting from the point testing of the higher-order polynomials) is diagonal. This explicit MOT scheme uses the same time step size as its implicit counterpart (which also uses Nystr€ om method in space) without sacrificing from the stability and the accuracy of the solution. Furthermore, it is faster than the implicit MOT scheme under low frequency excitation (i.e., for large time step sizes) since the matrix equation solved by the implicit scheme becomes fuller with the increasing time step size. Indeed, numerical results demonstrate the accuracy, efficiency, and applicability of the proposed explicit MOT solver.
The remainder of the paper is organized as follows. Section II A briefly formulates the TDKIE. Section II B provides the higher-oder spatial and temporal expansions used for representing the unknown velocity potential. Sections II C and II D detail the implicit and explicit MOT schemes, respectively, and Sec. II E compares the computational complexities of these two schemes under high and lowfrequency excitations (for small and large time steps). Sections III A and III B provide numerical results that demonstrate the accuracy and efficiency, and the applicability of the proposed explicit solver, respectively. Finally, Sec. IV summarizes the work detailed in the paper and outlines future research directions.

A. TDKIE
Let S represent the surface of an acoustic scatterer object residing in an unbounded homogeneous background medium with fluid density q 0 and wave speed c 0 . An acoustic wave with velocity potential u i ðr; tÞ and pressure field p i ðr; tÞ ¼ Àq 0 @ t u i ðr; tÞ is incident on S; it is assumed that u i ðr; tÞ and p i ðr; tÞ are band-limited and vanishingly small for t 0; 8r 2 S. Let u s ðr; tÞ and p s ðr; tÞ ¼ Àq 0 @ t u s ðr; tÞ represent the velocity potential and pressure field of the wave scattered from S. The total velocity potential is the sum of u i ðr; tÞ and u s ðr; tÞ : uðr; tÞ ¼ u i ðr; tÞ þ u s ðr; tÞ. Similarly, total pressure field is the sum of p i ðr; tÞ and p s ðr; tÞ: pðr; tÞ ¼ p i ðr; tÞ þ p s ðr; tÞ. One can obtain the TDKIE in (unknown) uðr; tÞ for the exterior problem (in the background medium, outside the volume enclosed by S) using the Kirchhoff-Helmholtz theorem, 30 @ t u i ðr;tÞ ¼ 1 2 @ t uðr; tÞ À @ t ð S uðr 0 ; tÞ Ãn 0 ðr 0 Þ Á r 0 GðR;tÞds 0 þ@ t ð Sn 0 ðr 0 Þ Á r 0 uðr 0 ; tÞ Ã GðR; tÞds 0 ; r 2 S: (1) Here, GðR; tÞ ¼ dðt À R=c 0 Þ=ð4pRÞ is the Green function in the background medium, R ¼ jr À r 0 j is the distance between the observation point r and the source point r 0 ; dðÁÞ is the Dirac impulse function, andn 0 ðr 0 Þ is the outward pointing unit normal to S at r 0 , and @ t and * denote the temporal derivative and convolution, respectively. Note that the first integral on the left hand side of Eq. (1) should be evaluated in the principal value sense. For an acoustically rigid object, the normal of the (total) velocity vanishes, 30 i.e., uðr; tÞ satisfies the boundary condition nðrÞ Á ruðr; tÞ ¼ 0; r 2 S; wherenðrÞ is the outward pointing unit normal to S at r. Inserting Eq. (2) into Eq. (1) yields the TDKIE for a rigid object, 10,30 1 2 @ t uðr; tÞ À @ t ð S uðr 0 ; tÞ Ãn 0 ðr 0 Þ Á r 0 GðR; tÞds 0 ¼ @ t u i ðr; tÞ: B. Higher-order spatial representation To solve Eq. (3) numerically, S is discretized into N p curvilinear triangular patches, and u(r, t) is spatially expanded/approximated using higher-order interpolation functions as [31][32][33] uðr; tÞ ¼ IðtÞ È É ði;pÞ # À1 ðrÞL ði;pÞ ðrÞ: Here, N u is the number of interpolation points on each patch, L ði;pÞ ðrÞ is the Lagrange interpolation function defined at r ði;pÞ (ith interpolation point on pth patch), 33 and IðtÞ È É ði;pÞ is the time-dependent unknown expansion coefficient at r ði;pÞ . In Eq. (4), #ðrÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð@ n r Á @ n rÞð@ g r Á @ g rÞ À ð@ n r Á @ g rÞ 2 q is the Jacobian of the coordinate transformation between the Cartesian coordinate system and a two-dimensional space with variables n and g. Here, derivatives @ n and @ g should be applied separately to the components of the threedimensional vector r. This quadratic transformation enables the description of the curvilinear triangular patch (defined by three nodes at its vertices and another set of three nodes at the mid-points of its edges) [ Fig. 1(a)] using a (flat) rightangle triangular patch generated in the (n, g) space [ Fig.  1(b)]. 3 Figures 1(a) and 1(b) also show the interpolation points r ði;pÞ in the Cartesian coordinate system and their mapped locations, which are denoted by r ðn;gÞ , in the (n, g) space, respectively.

C. Implicit MOT scheme
Inserting Eqs. (4) and (5) into Eq. (3), point testing the resulting equation in space and time, i.e., evaluating it at spatial points r ðj;qÞ ; j ¼ 1; …; N u ; q ¼ 1; …; N p , and at time samples t ¼ hDt; h ¼ 1; …; N t , yield the implicit MOT system, Here, the entries of the tested excitation vector V i h and the implicit MOT matrices Z imp hÀl are given by where S p is the surface of the pth patch, R ¼ r ðj;qÞ Àr 0 ; R ¼ jRj, The method that is used for computing the entries of Z imp hÀl is described in the Appendix.
After the matrices Z imp hÀl (including Z imp 0 ) and the vector V i h are computed and stored, unknown vectors I h ; h ¼ 1; …; N t , are obtained recursively by time marching, as briefly described next. First, I 1 at t ¼ Dt is found by solving the MOT system in Eq. (6) with right-hand side V i 1 . Then at t ¼ 2Dt; I 1 is used to compute the matrix-vector product Z imp 1 I 1 [there is only one term in the summation on the right hand side of Eq. (6)], which is subtracted from V i 2 to yield the complete right-hand side. I 2 is found by solving the MOT system in Eq. (6) with this right-hand side. Then, at t ¼ 3Dt, I 1 and I 2 are used to compute the summation Z imp 2 I 1 þ Z imp 1 I 2 , which is subtracted from V i 3 to yield the complete right-hand side. This permits the computation of I 3 and so on. At the end of time marching, all I h ; h ¼ 1; …; N t , are obtained known. Solution of the MOT system in Eq. (6), which is required by this time marching algorithm at every time step, is done using an iterative solver. In this work, a transpose-free quasi-minimal residual (TFQMR) scheme 36  are the solution vectors at two consecutive iterations (n and n -1) and v TFQMR is a user-defined tolerance/accuracy parameter.

D. The explicit MOT scheme
Inserting Eq. (4) into Eq. (3) and point testing the resulting equation in space, i.e., evaluating it at spatial points r ðj;qÞ ; j ¼ 1; …; N u ; q ¼ 1; …; N p , yield a linear system of ODEs, Here, the entries of time dependent vectors V i ðtÞ and V s ðtÞ, and the entries of the Gram matrix G are given by G f g ðj;qÞði;pÞ ¼ 1 2 # À1 ðr ðj;qÞ Þd ij d pq ; V i ðtÞ È É ðj;qÞ ¼ @ t u i ðr ðj;qÞ ; tÞ; Equation (9) is a system of ODEs and the samples of the unknown coefficient vector, i.e., I h ¼ IðhDtÞ, is obtained by integrating this ODE system in time using a PE(CE) m -type linear multistep scheme. This approach calls for sampling Eq. (9) in time, where h , one has to account for the retarded-time integral in Eq. (12); this is done by using temporal interpolation on where the entries of the explicit MOT matrices Z exp hÀl are given by Equation (14) is a system of ODEs, which relates samples of the unknown's temporal derivative, _ I h , to the samples of the unknown, I h . A PE(CE) m method enhanced with an SOR scheme 27,37 is used to for integrating Eq. (14) in time for the unknown samples I h ; h ¼ 1; …; N t . Steps of this method are described next.
Step 0: Before the time marching starts, compute and store G À1 .
Step 1: Compute the fixed part (the part that does not change within the steps of one time step) of the right-hand side of Eq. (14), Step 2: Predict I h using K past values of I l and _ I l , Step 3: Evaluate _ I h using V h exp and the predicted value of I h , Step 4: Set _ I and K past values of I l and _ I l , Step 4 Step 4.4: Check convergence, End for loop over n.
Step 5: Assuming convergence happens at the mth iteration, set Several comments about the explicit MOT scheme described above are in order.
(i) At Step 4.4, v PECE is a user-defined tolerance/accuracy parameter. (ii) At Steps 2 and 4.1, p and c are vectors of length 2K and 2K þ 1, which store the predictor and corrector coefficients of the PE(CE) m scheme. These coefficients can be obtained by either applying polynomial interpolation between time samples (e.g., Adam-Moulton, Adam-Bashforth, and backward difference methods 38 ) or using a numerical scheme that approximates the solution in terms of decaying and oscillating exponentials. 39 (iii) Initial values of I l and _ I l are assumed to be zero, i.e., I l ¼ 0 and _ I l ¼ 0, for l ¼ 0,…,K -1. This does not introduce any significant error in the solution since u i ðr; tÞ and p i ðr; tÞ are vanishingly small for t 0; 8r 2 S. If another type of excitation, which increases faster at the beginning of time marching, is used, the initialization of I l and _ I l can be done using the Euler method or spectral-deferred correction type methods. 38,40 (iv) Comparing the expressions of the entries of Z exp hÀl in Eq. (15) to those of Z imp hÀl in Eq. (8), it can be seen that Z exp hÀl and Z imp hÀl are the same except the "self-term" expression # À1 ðr ðj;qÞ Þ@ t TðtÞj t¼ðhÀlÞDt d ij d pq =2 presented in Z imp hÀl . This term contributes only to the diagonal entries of Z imp hÀl ; h À l ¼ 0; 1; …; P, where P is the order of T(t), since @ t TðtÞj t¼ðhÀlÞDt ¼ 0 for h À l 6 2 f0; 1; …; Pg and d ij ¼ 0 for i 6 ¼ j and d pq ¼ 0 for p 6 ¼ q. Since the same diagonal entries of Z exp hÀl are already non-zero, the sparseness levels of Z exp hÀl and Z imp hÀl are the same for all the values of l, l ¼ 1,…,h. The effect of this on the computational complexity of the implicit and explicit MOT schemes is discussed in Sec. II E. The method, which is used for computing the entries of Z imp hÀl and described in the Appendix, is also used for computing entries of Z exp hÀl .

E. Computational complexity analysis
Numerical results presented in Sec. III show that the explicit MOT scheme maintains its stability using a Dt as large as that would be used by the implicit scheme. The results also show that for the same Dt, both schemes achieve the same level of accuracy. Therefore, while comparing the computational complexities of these two schemes, it is assumed that they use the same Dt and they are executed for the same N t . To this end, in this section, the computational costs of one time step for both schemes are compared.
At a given time step, the implicit MOT scheme com- (16) and the fixed part of the right hand side of Eq. (14)]. The sparseness levels of Z exp hÀl and Z imp hÀl are the same for l ¼ 1; …; h as discussed in Sec. II D. This means that the costs of computing the summations in V imp h and V exp h are the same. It should also be noted here that this computation could be accelerated using the PWTD algorithm 11,12,21,22 or the FFT-based methods; 23-26 both the explicit and implicit MOT scheme benefit from the use of these methods. Therefore, with or without an acceleration method, the cost of computing V imp h and V exp h is the same and it is represented by C rhs .
At every time step, the implicit MOT scheme iteratively solves the matrix system in Eq. (6). The computational complexity of this operation scales as Here, N s ¼ N p N u represents the total number of unknowns, N imp iter is the number of iterations, F imp iter is the number of matrix-vector multiplications at each iteration, and c is the "denseness" factor of Z imp 0 , i.e., (average) number of nonzero entries in one row/column of Z imp 0 . The computational cost of the explicit MOT scheme is derived by counting the number of operations in steps 1-4 of the PE(CE) m scheme. The cost at step 1, i.e., the cost of computing V exp h , is already discussed above. The computational cost of the remaining steps 2-4 is given by Here, the first term is the cost of computing the right hand side of Eq. (17), the second term is the cost of computing the right hand side of Eq. (19) for a total of m iterations. The third term is the total cost of multiplications Z exp 0 I h in Eq. (18) and Z exp 0 I ðnÞ h (m times) in Eq. (21). Note that c is the denseness factor of Z exp 0 (same as that of Z imp 0 ). The last term is total cost of multiplication with G À1 in Eqs. (18) and (21).
To compare C imp and C exp , the values of m, K, c, N s ; N imp iter , and F imp iter are needed to be studied. In the case of high frequency excitations (small Dt; Dt ( D max =c 0 , where D max is the maximum distance between any two points in S), both Z imp 0 and Z exp 0 are very sparse, i.e., c ( N s . Therefore, C imp and C exp scale as C imp $ OðN imp iter F imp iter cN s Þ and C exp $ OðmKN s Þ þ OðmcN s Þ. This suggests that C exp % C imp (depending on the exact values of m, K, c, N imp iter , and F imp iter ), which means that both solvers will have similar speeds for high frequency excitation (since C exp þ C rhs % C imp þ C rhs ). This is indeed demonstrated by the numerical results presented in Sec. III A.
Matrices Z imp hÀl ¼ 0 and Z exp hÀl ¼ 0 for h À l > D max = ðc 0 DtÞ þ P. Consequently, as the frequency decreases (large Dt; Dt % D max =c 0 ), the number of non-zero Z imp hÀl and Z exp hÀl decreases. However, these non-zero matrices (including Z imp 0 and Z exp 0 ) become denser (even full matrices). For example, for D max =ðc 0 DtÞ < 1, all non-zero Z imp hÀl and Z exp hÀl are completely full. This means that C rhs (cost of computing V imp h =V exp h ) becomes significantly smaller than C imp and C exp and the latter can be used to predict the total cost of time marching. As Z imp 0 and Z exp 0 become denser, c $ N s . As a result, C imp $ OðN imp iter F imp iter N 2 s Þ and C exp $ OðmN 2 s Þ. Since m is smaller than N imp iter F imp iter (as shown by the numerical results presented in Sec. III A), C imp < C exp . As a result, for large Dt, the explicit MOT solver is expected to be faster than its implicit counterpart. Indeed, this is demonstrated by the numerical results presented in Sec. III A.

F. Comments
Several comments about the surface integral equation formulation described in Sec. II A and its discretization using the methods described in Sec. II B, II C, and II D are in order.
(1) The TDKIE in Eq. (3) is formulated to analyze scattering from an acoustically rigid object. A perfectly rigid scatterer can be considered as a mathematical idealization since it refers to an object made of an infinitely stiff material. However, it can be considered as a good approximation in scattering scenarios where an object with a high bulk modulus resides in a background medium with a significantly lower bulk modulus. For example, a steel submarine embedded in water can be approximated as an acoustically rigid scatterer since the bulk modulus of steel is roughly 70 times that of water. For scenarios that do not fall in this category, where the object has to be considered as a penetrable scatterer, one can formulate a coupled system of surface integral equations by combining exterior and interior scattering integrals. 30 (2) The TDKIE in Eq. (3) suffers from the well-known internal resonance problem. 41 However, the amplitude of spurious modes can be kept under control (for only TDSIEs but not their frequency domain versions) by increasing the accuracy of the MOT scheme. 42,43 These modes can be completely eliminated from the solution if one solves a combination of the TDKIE in Eq. (3) and its normal derivative (Burton-Miller formulation). 10 However, this formulation calls for significantly more complicated singularity extraction/cancellation techniques (especially for higher-order surface discretizations like the one described in Sec. II B) since the order of singularity in the kernel of the integral equation is increased. 44 (3) The MOT schemes developed for solving the TDKIE in Eq. (3) can be extended to solve the coupled system of equations to account for penetrable scatterers or the Burton-Miller surface integral equation to obtain a more accurate solution. These extensions do not call for significant changes in either the method used for the temporal discretization or the implementation of the time marching. However, one needs to significantly modify the method used for the spatial discretization to account for different integral equation kernels. These modifications will be formulated and described in a future publication.

III. NUMERICAL RESULTS
In this section, the accuracy, efficiency, and applicability of the proposed explicit MOT are demonstrated through the numerical analysis of transient acoustic scattering from rigid objects residing in an unbounded background medium.
In all examples, the object/scatterer is excited by a plane wave with velocity potential u i ðr; tÞ ¼ u i 0 G mod ðt Àk i Ár=c 0 Þ: Here, u i 0 ¼ 1 m 2 =s is the amplitude,k i is the direction of propagation of the incident plane wave, G mod ðtÞ ¼ cos ½2pf 0 ðt À t 0 Þe ÀðtÀt 0 Þ 2 =2w 2 is a modulated Gaussian pulse, f 0 is the pulse center frequency, t 0 is the pulse delay, and w is a measure of the pulse width. For all the examples considered here, an effective bandwidth for G mod ðtÞ is defined as f bw ¼ 3=ð2pwÞ. This specific definition f bw ensures that 0.003% of G mod ðtÞ's total power is within the frequency band ½f 0 À f bw ; f 0 þ f bw . In all examples, it is assumed that the scatterer resides under water and therefore c 0 ¼ 1500 m/s. The spatial discretization uses the Nystr€ om method with second-order Lagrange polynomials, i.e., N u ¼ 6. Tolerance/ accuracy parameters for the implicit and explicit MOT schemes are the same: After the time-domain simulations are completed and the unknown coefficients I h ¼ IðhDtÞ; h ¼ 1; …; N t are solved for and stored, one can obtain the time-harmonic total velocity potential uðr; f Þ at frequency f using Iðf Þ È É ði;pÞ # À1 ðrÞL ði;pÞ ðrÞ: Here Here, u s ðr; f Þ is the time-harmonic scattered velocity potential at frequency f, k 0 ¼ 2pf =c 0 is the wavenumber in the background medium, andk s ¼x sin h cos / þŷ sin h sin / þẑ cos / is the direction along which rðh; /; f Þ is computed. Note that the traditional definitions of the angles h and / in spherical coordinate system are used here: h is measured from the z-axis and / is measured on the xy-plane from the x axis.

A. Accuracy and efficiency
For this example, transient acoustic scattering from a rigid unit sphere is analyzed using the implicit and explicit MOT schemes. For the first set of simulations, the sphere surface is discretized into N p ¼ 2560 curvilinear triangular patches resulting in N s ¼ 15 360 number of unknowns. The temporal interpolation function T(t) is the fourth-order Lagrange polynomial, i.e., P ¼ 4. For the explicit MOT solver, the SOR coefficient a ¼ 1 (i.e., no SOR is applied) and the coefficients p and c are computed using the numerical scheme that approximates the solution in terms of decaying and oscillating exponentials, 39  SCS of the sphere is also obtained at the same frequency and angles from the (analytical) Mie series solution. 46 Let this value be represented by r Mie ðh; 0 ; f 0 Þ. Figure 3 where type 2 fimp; exp g; Dh ¼ 0:5 ; and / ¼ 0 . Figure  3(b) plots r imp err ðf Þ and r exp err ðf Þ versus f changing from 600 to 1150 Hz. It is clearly seen that both solvers provide accurate results within the effective band of the excitation.
For the next set of simulations, the sphere surface is discretized into N p ¼ 116, N p ¼ 396, and N p ¼ 1126 curvilinear triangular patches, resulting in N s ¼ 696, N s ¼ 2376, and N s ¼ 6756 numbers of unknowns, respectively. The temporal interpolation function T(t) is the third-order Lagrange polynomial, i.e., P ¼ 3. For the explicit MOT solver, the SOR coefficient is a ¼ 0.6 and the coefficients p and c are computed using the numerical scheme that approximates the solution in terms of decaying and oscillating exponentials, 39 resulting in a PE(CE) m scheme with K ¼ 22. For each mesh, four different simulations with different f 0 and w are executed using the explicit and implicit MOT solvers (a total of 24 solver simulations). For these sets of simulations,k i ¼ẑ; f 0 ¼ 20; 65; f 220; 350g Hz; f bw ¼ 0:3125f 0 ; w ¼ 3=ð0:625pf 0 Þ, and t 0 ¼ 6w; and the simulations are carried out for N t ¼ 600 time steps with Dt ¼ 1=ð32:5f 0 Þ. Different Dt results in a different value for the denseness of Z imp 0 and Z exp 0 : c ¼ f696; 106; 30; 22g; c ¼ f2376; 309; 58; 34g, and c ¼ f6756; 857; 130; 61g for N s ¼ 696, N s ¼ 2376, and N s ¼ 6756, respectively. Table I   As expected, T imp fix and T exp fix are almost the same. At low frequencies (i.e., large Dt), the T exp For these values of Dt; T exp PECE is almost one third of T imp TFQMR ; and since T exp PECE ) T exp fix and T imp TFQMR ) T imp fix , the explicit MOT solver is more than two times faster than the implicit MOT solver. For high frequency excitations (i.e., small Dt), both solvers require almost the same time. The measurement results provided in Table I verify the statements about computational complexities of the solvers given in Sec. II E. Table II compares r imp err ðf 0 Þ and r exp err ðf 0 Þ, and it shows that the explicit MOT solver is slightly more accurate than its implicit counterpart (which uses the same Dt) in all the simulations. The table also shows that when the number of unknowns is increased from N s ¼ 696 to N s ¼ 2376, the accuracy increases for the same f 0 (and also the same Dt). However, when the number of unknowns is increased from N s ¼ 2376 to N s ¼ 6756, there is no change in the accuracy.
This can simply be explained by the fact that for these spatial discretization levels, the accuracy is bounded by the temporal discretization and can, for example, be increased by reducing Dt and/or increasing P.

B. Submarine
For this example, transient acoustic scattering from a rigid submarine is analyzed using the implicit and explicit  MOT solvers. The submarine fits in a box of dimensions 32.1 Â 3.6 Â 6.5 m as shown in Fig. 4, Its surface is discretized into N p ¼ 3422 curvilinear triangular patches resulting in N s ¼ 20 532 number of unknowns. The temporal interpolation function T(t) is the third-order Lagrange polynomial, i.e., P ¼ 3. For the explicit MOT solver, the SOR coefficient a ¼ 1 (i.e., no SOR is applied) and p stores the coefficients of the Adams-Bashforth method and c stores the coefficients of the backward difference method resulting in a PE(CE) m scheme with K ¼ 4. 38

IV. CONCLUSION
A fully explicit MOT scheme to solve the TDKIE for analyzing transient acoustic scattering from rigid objects is developed. The velocity potential on the surface of the scatterer is expanded using a higher-order polynomials in space. Inserting this expansion in the TDKIE and point-testing the resulting equation (i.e., Nystr€ om method in space) yield a system of ODEs. This system is integrated in time for the unknown expansion coefficients using a PE(CE) m method. The resulting time marching scheme does not call for matrix equation solution since the Gram matrix resulting from the point-testing of the polynomials is diagonal. Numerical results demonstrate that the proposed explicit MOT scheme is significantly faster than its implicit counterpart under low frequency excitation (i.e., for large time steps) while maintaining the accuracy and the stability of the solution.
Extension of the proposed explicit solver, which could be used in the solution of resonance-free Burton-Miller surface integral equation, is underway. thank the King Abdullah University of Science and Technology Supercomputing Laboratory (KSL) for providing the required computational resources.

APPENDIX: NUMERICAL EVALUATION OF MATRIX ELEMENTS
The evaluation of the first term in the right hand side of Eq. (8) is rather straightforward and its value depends only on # À1 ðr ðj;qÞ Þ and @ t TðtÞj t¼ðhÀlÞDt .
The second term in the right hand side of Eq. (8) is equal to Z exp hÀl È É ðj;qÞði;pÞ and is evaluated as described next. When patches q and p are not within close proximity of each other (e.g., they are more than two patches apart), the evaluation of Z where R ¼ r ðj;qÞ À r ði;pÞ , R ¼ jRj, and x ði;pÞ is the quadrature weight at the interpolation point r ði;pÞ . Note that r ði;pÞ exactly coincides with the Gaussian quadrature points. 33 Also, note that r ði;pÞ (and Gaussian quadrature points) are generated on the rightangle triangle in the (n, g) space and mapped back to Cartesian coordinate system where the curvilinear triangular patch resides. When patches q and p are not within close proximity of each other but do not overlap, a higher degree Gaussian quadrature is used to evaluate the integral, Z exp hÀl È É ðj;qÞði;pÞ ¼ À x L ði;pÞ ðr Þn 0 ðr Þ Á R 8pR 3 where R ¼ r ðj;qÞ À r ; R ¼ jRj; x is the quadrature weight at the quadrature point r , and N is the number of Gaussian quadrature points. Note that the term L ði;pÞ ðr Þ is needed since r and r ði;pÞ do not overlap. Also note that, like before, these quadrature points are generated on the right-angle triangle in the (n, g) space and mapped back to Cartesian coordinate system as r . When patches q and p overlap, the integrand becomes singular. In this case, the integral is evaluated in two parts, Â @ t TðtÞ ½ t¼ðhÀlÞDtÀR=c 0 ds 0 ¼ I 1 þ I 2 : (A3) The first integral I 1 has a singularity of 1/R and is evaluated after Duffy transformation is applied. 47 The second integral I 2 has a singularity of 1=R 2 . To account for this singularity, the integral in the (n, g) space is expressed in terms of q and h, variables of the local polar coordinate system centered at r ðn;gÞ [r ðj;qÞ transformed into the (n,g) space] (Fig. 8 Here, ds 0 ¼ #ðr 0 Þdndg; dndg ¼ qdqdh; qðhÞ is the distance between r ðn;gÞ and the boundary of the right-angle triangle along the direction of h [ Fig. 8(b)]. F(q, h) in Eq. (A4) has a singularity of 1/q, which is "extracted" as described below. 48 F(q, h) can be expressed as where F -1 (h) is a real function of h, which satisfies Ð 2p 0 F À1 ðhÞdh ¼ 0. The first term of the right hand side of Eq. (A5) is subtracted and added (after being integrated) back to I 2 (Ref. 48), All the integrands in Eq. (A6) are regular and are numerically evaluated using a Gauss-Legendre quadrature.