Cagniard-DeHoop Technique-Based Computation of Retarded Partial Coefficients: The Coplanar Case

Efficient computation of partial elements plays a key role in the Partial Element Equivalent Circuit (PEEC) method. A novel analytical method for computing retarded partial coefficients based on the Cagniard-DeHoop (CdH) technique is proposed. The methodology is first theoretically developed and then illustrated on the computation of a surface retarded partial coefficient pertaining to two coplanar rectangular surface elements. An efficient way for incorporating loss mechanisms in the time domain (TD) via the Schouten-Van der Pol theorem is proposed. Illustrative numerical examples demonstrating the validity of the introduced solution are given.


I. INTRODUCTION
The PEEC method is an integral-equation technique (e.g. [1], [2]) that is capable of translating an EM scattering problem to an equivalent-circuit representation [3]. Electric and magnetic field couplings are represented separately with partial potential coefficients (P) [4], [5], partial inductances (Lp) [6], [7], as well as resistances. The PEEC method is well suited for both high and low frequency problems down to dc [8], [9]. If a volume-equivalence formulation is assumed [10], the free-space full-wave Green's function (e −ikR /R with k = ω/c 0 , with c 0 the speed of light in the free space) is to be considered in the computation of partial elements. If the surface-equivalence principle is adopted, the full-wave Green's function within homogeneous lossy dielectric/magnetic materials is to be used along with its curl [11] in the computation of interaction integrals describing the electric and magnetic field coupling. Potential integrals over triangular domains in the static limit (k → 0) have recently been analyzed in [12].
When the volume-equivalence is used to derive the PEEC formulation, the exponential term in the full-wave Green's The associate editor coordinating the review of this manuscript and approving it for publication was Su Yan . function corresponds to a retardation. In this case, the exponential term was moved outside of the integrals by using the distance between the center of the source and the observation point [13]. Such an approximation is possible for moderate frequency problems [14]- [16] and allows to avoid to integrate the frequency dependent matrix elements at each frequency. However, this does not lead to highly accurate results for very high frequencies. A better approximation can be obtained by resorting to a Taylor series expansion of the Green's function [17]. Although more accurate to catch the frequency-dependence of the Green's function and more effective than the brute-force numerical integrationbased approach, the resulting partial elements are easily computed in the frequency domain but their TD counterpart is not straightforward and would require auxiliary equations to implement higher order terms of the Taylor's expansion.
In this paper, we shall provide a fundamentally novel analytical approach for computing retarded coefficients of a PEEC model based on the CdH technique [18]. The methodology is illustrated on the computation of a retarded partial coefficient pertaining to two coplanar rectangular surface elements. It demonstrated that in such a case the surface retarded potential coefficients can be derived analytically in the TD in terms of elementary functions only. Moreover, a general procedure based on the Schouten-Van der Pol theorem [19, p. 1056] enabling the inclusion of general loss mechanisms is proposed.

II. PROBLEM FORMULATION
A PEEC model is represented through a set of partial elements, the value of which is found upon evaluating spatial integrals over the domains of interacting discretization elements. For the sake of simplicity, we shall limit ourselves in this work to the interaction of two coplanar rectangular surface elements (see Fig. 1). Other configurations such as rectangular surfaces on distinct parallel or orthogonal planes or parallelepipeds will be considered in future reports. To localize the position in the problem configuration, we employ coordinates {x, y, z} with respect to an orthogonal Cartesian reference frame with the origin O and the standard basis {i x , i y , i z }. Consequently, the position vector is r = xi x + yi y + zi z . The time coordinate is denoted by t. Finally, * t denotes the time-convolution operator, H(x) represents the Heaviside unit-step function and δ(x) denotes the 1-D Diracdelta distribution.
With reference to [3, eq. (D.1)], we shall study a surface retarded partial coefficient expressed through a 2-D integral , where x > 0 and y > 0 denote the spatial discretization steps in the x-and y-direction, respectively. Furthermore, s is the (real-valued and positive) Laplace-transform parameter, S m,n = x y is the surface area of domain A m,n , and is the free-space Green's function pertaining to a homogeneous, isotropic and loss-free medium described by its (scalar and real-valued) electric permittivity and magnetic permeability µ with the corresponding EM wave speed c = ( µ) −1/2 > 0. As demonstrated in Sec. IV, conductive electric and linear hysteresis magnetic losses can subsequently be incorporated via the Schouten-van-der-Pol theorem [19, p. 1056] by replacing s with [(s + α) (s + β)] 1/2 , where α and β are non-negative constants (see [19, sec. 26.5]).

III. PROBLEM SOLUTION
The retarded partial coefficient as expressed through Eq. (1) will be next reformulated via a spatial Fourier representation. To that end, the Green's function in the s-domain is represented via the wave slowness representation aŝ has the meaning of the wave slowness normal to surfaces A m,n . Making use of the representation (3) in Eq. (1), we get where i 0 (x) denotes the modified spherical Bessel function of the first kind. Upon expanding the product of the modified spherical Bessel functions into their exponential factors it is found thatP mn (s) can be expressed in terms of the generic integral analyzed in Appendix VI. Accordingly, we may writê Transforming the latter to the TD, we end up with where the TD function I (x, y, t) is specified in the Appendix via Eq. (21). It is worth mentioning that Eqs. (6) and (7) are valid for arbitrary placed coplanar regions and, thus, they hold not only for mutual partial elements but also for the self-interaction term which is typically the most challenging part to be computed because of the singularity. Other geometrical configurations will be considered in forthcoming works.

Loss effects can be incorporated via the Schouten-Van der
Pol theorem [19, p. 1056] for the replacement of s with (s + α) 1/2 (s + β) 1/2 as described in the present section. This can be seen upon observing that losses can accounted for by substituting in Eq. (1) the free-space Green's function of the dissipative 3-D scalar, modified Helmholtz equation (cf. Eq. (2)) where α = σ/ε, β = /µ are non-negative constants representing the effect of conductive electric and linear hysteresis magnetic losses, respectively. Consequently, assuming that P mn (t) is available (see Eq. (7)), the retarded coefficient affected by the losses, say P † mn (t), can be represented viâ (9) Since the original of exp[−(s + α) 1/2 (s + β) 1/2 τ ] is known in closed form, the desired TD function immediately follows as where I 1 (x) is the modified Bessel function of the first kind. Apparently, if α = β, then P † mn (t) = P mn (t) exp(−αt), as expected. Thus, in line with the result represented by Eq. (10), the inclusion of rather general losses amounts merely to evaluating an additional integral with the bounded range of integration.

V. NUMERICAL EXAMPLE
A popular approximation in PEEC implementations is the midpoint one that amounts to replacing |r − r | in Eq. (1) by This simplification makes it possible to obtain an approximation of P mn (t) at once. We next discuss the corresponding results concerning both loss-free and dissipative Green's functions as given by Eqs. (2) and (8), respectively.

A. LOSS-FREE COEFFICIENT
Application of the midpoint approximation to Eq. (1) with (2) yields a simple s-domain expression that can be transformed to the TD as Equation (11) will next be compared with the exact solution represented by Eq. (7). For this purpose, we define the following parameter which, under the midpoint approximation, takes the value one.
In the following examples, we choose x = y = 1.0 m, which is supposed to be relatively small (with respect to c× (incident wave pulse time width) or the wavelength at the frequency of analysis). Without any loss of generality, the center of A m is located at the origin, that is, x m = 0 and y m = 0. In the first example, we take x n = y n = 5 x , so that surfaces A m and A n are relatively far apart, namely, r mn = 5 √ 2 x 7.07 x . The corresponding coefficient (scaled by 4π r mn /c) is shown in Fig. 2a. Apparently, the resulting pulse is relatively narrow, thus resembling the Dirac-delta impulse. Also, the corresponding value of I as found using the trapezoidal rule is approximately 1.0017, which is very close to the value one pertaining to the midpoint approximation. In the second case, we choose x n = y n = x , so that the surfaces are relatively near each other, namely, r mn = √ 2 x 1.41 x . In such a case, the computed pulse is not relatively narrow anymore, which indicates inadequacy of the midpoint approximation. Finally note that the corresponding value of parameter I is 1.0592.

B. DISSIPATIVE COEFFICIENT
The midpoint rule applied to Eq. (1) with the dissipative Green's function (8) yields the approximate coefficient which consists of the impulsive (instantaneously-reacting) part and the relaxational one forming the tail of the TD coefficient. Apparently, if both α = 0 and β = 0, which represents the loss-free case, then Eq. (13) boils down to (11). Again, for the sake of comparison, we define (cf. Eq. (12)) which, under the midpoint approximation, can be expressed as The integral in Eq. (15) can be readily evaluated numerically. Again, it is noted that I † takes the value one for the loss-free case.
The dissipative coefficient will be evaluated for two distances r mn = 5 √ 2 x and r mn = √ 2 x , with x = y = 1.0 m, again. First, we consider relatively weak electric-type dissipation represented by r mn α/c = 0.5 and r mn β/c = 0. In this case, the midpoint approximation of parameter I † as given by Eq.  Finally, Fig. 4 shows the computed TD coefficients for relatively strong dissipation described by r mn α/c = 2.5 and r mn β/c = 0.1. Apart from the significant attenuation of the pulse amplitudes, we may observe the change of the resulting pulse shapes, in particular, in their late-time parts. The slow decline of the pulse tails is attributed to the relaxational part of the TD coefficients that is proportional to |β − α| (see Eqs. (10) and (13)). Note that the value of I † as predicted by the midpoint approximation via Eq. (15) is I † 0.40467 in this case.

C. CAPACITIVE COUPLING IN THE FREE SPACE
In the last example, the capacitive coupling between regions m and n is considered. They are assumed as squares of size 1 cm on the same plane. The coefficient of potential describing their full-wave electric field coupling in the free space is computed as which corresponds to the CdH -based exact solution. It admits three approximations:   • Center-to-center (CC): P FW mn (t) 1/(4π r mn ) δ(t − r mn /c).
where P QS mn is the static coefficient of potential [3]. Region n is assumed to be excited by a charge pulse q n (t) of amplitude 1 pC with rise and fall times of 6.67 ps. The self induced potentials are computed using the proposed CdHbased approach, the QS and the QS − delayed approximations through convolution mn (t) = P FW mn (t) * t q n (t), n = 1, m = 1, 2. For the QS − delayed approximation, it has been assumed r mn = x = y = 1 cm. Figure 5 shows the self induced potential on region m. It is worth  to observe that the causality is correctly reproduced, while the QS and QS − delayed are inaccurate. Then, region n is moved along the diagonal at a distance r mn = 1.41 x and r mn = 14.142 x in order to compute the near-field and far-field coupling. Figures 6 and 7 show the mutual induced potentials. As expected, when the two regions are close, the effects of the finite size of the two regions are clearly visible and captured by the proposed CdH -based approach, while the approximated solutions QS, QS − delayed and CC fail. When the two regions are far apart, the QS −delayed and CC approximated solutions are quite close to the proposed CdH -based one because the geometrical dimensions are less important and the propagation delay is well approximated by using the center-to-center distance. On the contrary, the QS approach results to be non-causal since this approximation neglects the propagation delay.

VI. CONCLUSION
A fundamentally new methodology based on the Cagniard-DeHoop technique for calculating TD retarded coefficients VOLUME 8, 2020 has been proposed. For the case of two coplanar discretization surfaces, it is demonstrated that this technique yields an exact, closed-form analytical TD expression for the TD surface retarded coefficient of potential. The inclusion of dissipation mechanisms is accomplished with the aid of the Schouten-Van der Pol theorem. Future research directions will address an extension of the methodology to more general problem configurations including noncoplanar and orthogonal elementary surfaces.

APPENDIX THE ON-PLANE CASE
The TD counterpart of the retarded coefficient under consideration (see Eq. (1)) is given by Eq. (7) in terms of TD functions I (x, y, t). This calls for the inversion of the following integral representation for x ∈ R, y ∈ R and {s ∈ R; s > 0}, where K 0 and S 0 are the integration paths extending along Re(κ) = 0 and Re(σ ) = 0, respectively, that are indented to the right with semi-circular arcs with centers at the origins and vanishingly small radii (see Fig. 8). To obtain the TD original of Eq. (17), we shall pursue the CdH approach closely described in [20]. Accordingly, we first deform S 0 into the CdH-path, say L ∪ L * (here * denotes the complex conjugate), where for {1 ≤ u < ∞} with (κ) = (1/c 2 − κ 2 ) 1/2 , thus representing the loop around the branch cuts extending along { (κ) < |Re(σ )| < ∞, Im(σ ) = 0}. The deformation is permissible by virtue of Jordan's lemma and Cauchy's theorem. Consequently, introducing u as the new variable of integration, the inner integral with respect to σ can be expressed as where we have combined the contributions from L and L * and accounted for the non-zero contribution of the pole singularity at σ = 0. Substitution of Eq. (19) back in (17) leads tô whereÛ (x, y, s) andV (x, y, s) are transformed to the TD in the following subsections. Employing these results, we finally end up with which is the main result of the present work, with r = (x 2 + y 2 ) 1/2 > 0. This result can be used in Eq. (7) to obtain the desired retarded PEEC coefficient. Its limiting cases follow as and a similar expression applies to I (0, y, t). Finally, upon taking x → 0, we get which is useful to computing the coefficient pertaining to overlapping surface elements (cf. Eq. (7)).

GENERIC INTEGRALÛ
The first integral to be handled has the following form To obtain the TD counterpart ofÛ (x, y, s), we deform K 0 to the CdH path that is defined via where τ is the real-valued and positive (time) parameter. The deformation is permissible, again, as the condition for Jordan's lemma to apply is safely met. Upon solving (25) for κ, we find that the equality is met along the hyperbolic arcs, say G ∪ G * , parametrized via for all τ ≥ r(u)/c, where r(u) = (x 2 + u 2 y 2 ) 1/2 > 0. Further, introducing τ as the variable of integration and combining the contributions from G and G * , we arrive at where in the integrand we take the values along G and U 0 (x, y, s) denotes the contribution from the pole at κ = 0. The latter is given bŷ Next, in the second term of Eq. (27), we change the order of integration according to in which r = r(1) and U (τ ) = (c 2 τ 2 − x 2 ) 1/2 /|y|. Subsequently, the second integral with respect to u is carried out analytically, which after lengthy but straightforward algebra yieldŝ where the first integral on the right-hand side of Eq. (30) is, in fact, the pole contribution (28), where we substituted τ = u|y|/c for all τ ≥ |y|/c. The integrals in (30) have the form of Laplace integrals, which makes it possible to employ Lerch's uniqueness theorem [21, appendix] and get the original U (x, y, t) at once. Owing to the factors s −2 and s −3 in front of the integrals, the thus constructed U (x, y, t) is composed of convolution-type integrals. Since they are amenable to analytical solution, we finally end up with which is used to derive the main result represented by Eq. (21).

GENERIC INTEGRALV
The second integral to be transformed to the TD has the following form In the first step, again, the indented integration contour K 0 is deformed to G ∪G * , which now represents the loop encircling the branch cuts {1/c < |Re(κ)| < ∞, Im(κ) = 0} in the complex κ-plane (see Fig. 8b), that is (cf. Eq. where the first term on the left-hand side represents the contribution from the pole singularity at κ = 0. The inversion of the first term is straightforward, while the transform of the second term yields a convolution-type integral, again. Carrying out the latter analytically, we arrive at which is finally used to derive the main result represented by Eq. (21).