Pulsed Electromagnetic Field Transmission through a Small Rectangular Aperture: A Solution Based on the Cagniard–DeHoop Method of Moments

: Pulsed electromagnetic (EM) ﬁeld transmission through a relatively small rectangular aperture is analyzed with the aid of the Cagniard–deHoop method of moments (CdH-MoM). The classic EM scattering problem is formulated using the EM reciprocity theorem of the time-convolution type. The resulting TD reciprocity relation is then, under the assumption of piecewise-linear, space–time magnetic-current distribution over the aperture, cast analytically into the form of discrete time-convolution equations. The latter equations are subsequently solved via a stable marching-on-in-time scheme. Illustrative examples are presented and validated using a 3D numerical EM tool.


Introduction
Apertures in conducting planes are frequently encountered in the form of windows, cracks around doors, coupling slots in microwave devices, or imperfect seams between two metallic plates (e.g., ( [1], Chapter 4) and ( [2], Sections 6.7 and 7.9)). In order to quantify both intentional and undesired effects of such apertures, their EM scattering is a major concern in applied electromagnetics.
Wavefield penetration through an aperture in a thin screen is a classical problem of wavefield physics, solutions of which can be traced back to the work of Lord Rayleigh (see [3] and ( [4], § 6.1)). His solution applying to the plane-wave diffraction by a relatively small (with respect to the wavelength) aperture was later extended by Bouwkamp and Van Bladel [5][6][7]. Without restricting himself to the plane-wave excitation, Bethe [8] demonstrated that the EM scattering by a small aperture can be, in certain circumstances (see ([5], § 9)), attributed to the action of equivalent electric and dipole moments [9,10]. Whenever the maximum linear dimension of the aperture is not sufficiently small with respect to the operating wavelength, however, the approximate aperture-diffraction models due to Rayleigh and Bethe are no longer applicable and one has to resort to more general formulations. A way out of the difficulty has been offered by Levine and Schwinger through their variational formulation [11,12], thus extending the steady-state solution to a wide range of frequencies.
The available rigorous solutions of the aperture diffraction problem were derived under the assumption of sinusoidally in time-varying wavefields. However, owing to the widespread use of communication and radar systems relying on the transmission, detection, and subsequent interpretation of digital signals, the steady-state assumption may no longer be computationally efficient and/or physically legitimate. Therefore, in the present work, we remove the restriction and solve the diffraction wavefield problem in its (original) spacetime domain. The vast majority of previous works aiming at TD solutions rely more or less on the pertaining steady-state results-exact, purely TD analytical solutions of the aperture diffraction are lacking in the literature on the subject. Indeed, the pertinent solutions can be, broadly speaking, divided into two categories. The first line of reasoning employes purely numerically approaches such as the inverse FFT algorithm [13,14] or the finite-difference TD technique [15,16]. The second approach yields TD approximate field expressions based on the use of Bethe's model [17] or Kirchhoff's approximation ( [18], Section 5.5.1). To deal with the issue, the classic CdH technique [19] (see also [20][21][22][23], for example) and the TD EM reciprocity theorem of the time-convolution type ( [24], Section 28.2) (see also ([25], Section 5.2), ( [26], Section 1.4.1), and [27]) are in the present work combined to introduce a novel, rigorous TD integral-equation approach to analyzing the pulsed EM scattering by a relatively small rectangular aperture in a PEC screen. The new solution strategy, to be referred to as the CdH-MoM, has been previously applied to TD performance studies of cylindrical and planar antennas [28] and, more recently, to the TD EM scattering analysis of transmission lines in the presence of thin sheets [29] and EM metasurfaces [30,31].

Problem Definition
The problem configuration under analysis is shown in Figure 1. The position in this configuration is specified by the coordinates {x, y, z} with respect to an orthogonal, Cartesian reference frame with the origin O and the three base vectors, {i x , i y , i z }, forming in the indicated order a right-handed triad. The position vector is then written as r = xi x + yi y + zi z . The time coordinate is t, and the continuous time-convolution operator is represented by * . The Heaviside unit-step function is denoted by H(t), and the Dirac-delta distribution is δ(t). Partial differentiation is denoted by ∂, which is supplied with the pertaining subscript. For instance, to differentiate with respect to x, we use ∂ x . We shall analyze the EM-field penetration through a bounded aperture in an unbounded PEC screen. The PEC plane lies in z = 0, and the aperture occupies a rectangular In the present analysis, it is assumed that the maximum dimension of the aperture is relatively small with respect to the spatial support of the exciting EM pulse. The surrounding medium is assumed to be linear, isotropic, and loss-free in its EM properties. It is described by (real-valued, positive scalar) electric permittivity 0 and magnetic permeability µ 0 . The corresponding EM wave speed is c 0 = ( 0 µ 0 ) −1/2 > 0, and the wave admittance is denoted by The aperture is supposed to be irradiated by the (causal) incident EM wave field (denoted by superscript i ), which is specified by its electric-field strength, E i (r, t), and magneticfield strength, H i (r, t). The incident EM field is not necessarily a plane-wave. To account for the presence of the aperture, we next define the scattered EM wave field (denoted by superscript s ) as the difference between the total fields, {E, H}(r, t), and excitation EM wave fields (denoted by superscript e ), viz.
for all r ∈ R 3 and t > 0, where the excitation fields have the meaning of total fields in the absence of the aperture (i.e., with the "short-circuited aperture"). Since i z × E e (r, t) = 0 over the entire z = 0 plane for all t > 0, the use of the explicit-type boundary condition applying to the total EM field on the PEC screen in Equation (1) implies that i z × E s (r, t) = 0 on the PEC screen for all t > 0.

Time Domain Problem Formulation
The EM scattering problem under consideration is formulated using the TD Lorentz reciprocity theorem ( [24], Section 28.2). To that end, the reciprocity theorem is applied to the scattered EM wave fields and the (causal) testing EM wave fields (denoted by superscript T ), which are generated by the testing magnetic-current surface density, ∂K T (r, t), whose spatial support is the surface of the aperture, supp(∂K T ) = S. Then, upon applying the TD reciprocity theorem to the upper half-space, z > 0, while enforcing the surface boundary condition applying to (the tangential component of) the scattered field, ∂K s = E s × i z , on the PEC plane, we arrive at where · represents the standard inner product of two vectorial quantities ( [26], Equation (1.2)), * denotes the continuous time-convolution operator ( [26], Equation (1.11)), and the superscript + indicates that we approach the surface of the aperture, S, from above as z ↓ 0. Since both states are causal and each other's adjoint, the contribution from the bounding sphere "at infinity" vanishes ( [26], Section 1.4.3). In the second step, the TD reciprocity theorem is applied to the lower half-space, z < 0, which leads to a similar relation: where − indicates that we approach the surface of the aperture from below as z ↑ 0.
Recalling the definition of the excitation field, we may use Equation (1) with applying to all r ∈ S and t > 0, to combine relations (2) and (3). This yields the desired TD reciprocity relation: The final TD relation (4) will next be solved for the equivalent magnetic-current spacetime distribution, ∂K s (x, y, t), as induced in the aperture. Owing to the relatively small dimension of the aperture, we may assume that the magnetic-current spatial distribution has the following form: where k x (t) and k y (t) are the magnetic-current pulse shapes (to be computed), is the triangular function, and represents the rectangular function. The spatial distribution of the magnetic-current as represented by Equations (5) and (6) is chosen such that the pertinent end conditions, ∂K s x = 0 and ∂K s y = 0 at x = ±∆ x and y = ±∆ y , respectively, are satisfied for all t > 0. Along their transverse direction, the magnetic-current components are assumed to have the uniform spatial distribution. Consequently, the magnetic-current flowing in parallel to the edge does not exhibit the inverse-square-root singularity as required by the pertaining edge condition ([1], Chapter 4). It can be expected, however, that this choice will have a negligible impact on the estimation of far-field characteristics [32]. Finally, it remains to specify the testing currents on the right-hand side of the TD reciprocity relation. To this end, we may choose the "razor-type" testing functions It is noted that the thus-chosen testing currents satisfy the end conditions, ∂K T x = 0 and ∂K T y = 0 at x = ±∆ x and y = ±∆ y , respectively.

Problem Solution
The problem will be solved via the CdH joint transform technique [19]. Accordingly, we combine a unilateral Laplace transformation, i.e., for {s ∈ R; s > 0}, by virtue of Lerch's uniqueness theorem ( [26], Appendix), with the wave slowness representation taken along the surface parallel to the PEC screen: where κ and σ are (imaginary-valued) slowness parameters in the xand y-direction, respectively. As a matter of fact, Equation (12) is a two-dimensional Fourier inversion integral, where the (real-valued and positive) Laplace-transform parameter, s, plays the role of a scaling parameter. This representation entails the properties ∂ x → −sκ and ∂ y → −sσ. Under the wave slowness representation, the TD reciprocity relation (4) can be written as The (tangential components of the) transform domain testing fields as generated by the equivalent magnetic-current surface density, ∂K T , can be determined from the transform domain counterparts of standard source-type representations for the electric-and magnetic-field strengths in a loss-free medium (cf. ( [24], Equations (26.10-14) and (26.10-15) where Ω 2 0 (κ) = 1/c 2 0 − κ 2 , and the transform domain Green's function,G, follows as the bounded solution of Using the limit from Equation (18) in Equations (14) and (15), it can be easily verified that the transform domain testing fields satisfy the excitation conditionsẼ T The (bounded) transform domain Green's function can be written asG for (19) can be subsequently used in the transform domain expressions (16) and (17) to determine the testing fields on the left-hand side of the reciprocity relation (13). The latter relation is subsequently solved in the TD. To that end, we first expand the unknown pulses k x,y (t) (see Equations (5) and (6)) in a piecewise linear manner: where v k;x and v k;y are (yet unknown) coefficients and Λ k (t) represents the temporal triangular function: along the discretized time axis {t k = k∆ t ; ∆ t > 0, k = 1, 2, · · · , M}. Finally, upon substituting the transform domain images of Equations (5) and (6) with (20) and (21), (9) and (10) in the reciprocity relation (13), we end up with a system of equations in the s domain. Its constituent can be transformed to the original TD analytically via the CdH technique, which yields the following system of discrete time-convolution equations: where Y k = Y (t k ) represents a [2 × 2] TD admittance array that is specified in Appendix A. Furthermore, V k = [v k;x , v k;y ] T is the [2 × 1] array of the unknown coefficients (see Equations (20) and (21)), and finally, I k = I(t k ) = [I x (t k ), I y (t k )] T denotes a TD excitation [2 × 1] array, which, as a matter of fact, corresponds to the interaction term on the right-hand side of the TD reciprocity relation (4). Its specific form will be given in Section 5. Once the admittance and excitation arrays are specified, Equation (23) can be readily solved via the marching-on-in-time technique. This way leads to the following step-by-step updating scheme: which yields the desired coefficients for all m = {1, 2, · · · , M}. Illustrative numerical examples validating the TD solution (24) are presented in the following section.

Illustrative Examples
Throughout this section, it is assumed that the rectangular aperture is irradiated by a uniform E-polarized TD plane-wave: where the unit vector in the direction of polarization, α, can be expressed via the azimuthal angle, φ, as α = sin(φ)i x − cos(φ)i y . The pertaining slowness parameters are then described by κ 0 = cos(φ) sin(θ)/c 0 , σ 0 = sin(φ) sin(θ)/c 0 and γ 0 = cos(θ)/c 0 , where θ denotes the polar angle. The (causal) plane-wave signature, e i (t), has the property e i (t) = 0 for all t < 0. With reference to the right-hand side of the TD reciprocity relation (4), the corresponding tangential components of the magnetic-field strength at the plane of aperture, z = 0, then follow as To determine the corresponding TD excitation array I(t) (see Equation (23)), we may use either Equations (9) and (10) with (26) and (27) on the right-hand side of Equation (4) or, equivalently, evaluate the right-hand side of the transform domain reciprocity relation (13). Both ways yield the following elements of the excitation array: The limits pertaining to the normal incidence characterized by θ = 0, implying κ 0 = σ 0 = 0, follow as For the presented example, we take 2∆ x = 0.10 m and ∆ y = ∆ x /2 with θ = 0 and φ = π/4. Furthermore, the plane-wave signature has the shape of a bipolar triangle, which can be described by where we take e m = 1.0 kV/m and c 0 t w = 1.0 m (see Figure 2). Consequently, c 0 t w = 10 max(2∆ x , 2∆ y ), which implies that the aperture is relatively small, as assumed in our analysis. In the first step, the marching-on-in-time solution (24) is applied to calculate the magnetic-current coefficients, v k;x and v k;y , which, through Equations (20) and (21) with (5) and (6), determine the space-time distribution of the magnetic-current surface density induced throughout the aperture. Its tangential components can be expressed using the actual electric-field strengths in the aperture as ∂K s x (x, y, t) = E y (x, y, 0, t) and ∂K s y (x, y, t) = −E x (x, y, 0, t). The electric-field pulses at the center of the aperture as calculated with the aid of the CdH-MoM method and via the finite-integration technique (FIT) as implemented in CST Microwave Studio are shown in Figure 3a,b. As can be seen, the results agree very well. To quantify the deviation between the signals, we evaluated the normalized root-mean-squared error according to ERR = ∑ M k=1 (f k − f k ) 2 /M/∆ f , wheref k and f k represent the time samples of the FIT and CdH-MoM signals, respectively, and ∆ f = f max − f min . The error corresponding to the E x -field component (see Figure 3a) is about ERR 1.5%, while ERR 3.2% for the E y -field component (see Figure 3b). The higher error in the E y -field component can be attributed to the postulated magneticcurrent space-time distribution (see Equations (5) and (6)), which is less accurate along the relatively longer side of the aperture.
Typical computational times to obtain the TD responses shown in Figure 3 are about 20 s using a non-optimized MATLAB CdH-MoM implementation and about 30 min using CST Microwave Studio . An even more striking difference is in the number of required unknowns and accompanying memory requirements. While our dedicated CdH-MoM computational model solves the system for two unknown quantities only (see Equations (5) and (6)), the FIT model consists of about 600 thousand discretization elements. The calculations were conducted on a standard laptop with Intel(R) Core(TM) i7-10510U CPU @ 1.80 GHz.
To indicate the range of applicability of the presented computational model, we increased the aperture's length along the x-direction such that 2∆ x is now a fifth of the excitation pulse's spatial support c 0 t w , i.e., 2∆ x = 0.20 m = c 0 t w /5. In addition, we took ∆ y = ∆ x /10, so that the model represents a relatively narrow slit in the PEC screen. The remaining parameters were kept the same. Figure 4 shows the corresponding pulse shapes of the E y -field at the center of the aperture. As the length of the slit is no longer sufficiently small with respect to c 0 t w , discrepancies between the signals, quantified by ERR 5.8% are apparent. While this result can still be useful for initial estimates, improving the accuracy through the incorporation of the piecewise-linear spatial expansion of the axial current density (e.g., ( [28], Section 14.3)) is advisable in this case.    Once the space-time distribution of the equivalent magnetic-current surface density in the aperture is known, it can be used to evaluate the scattered and, hence, total EM-fields via Equation (1). To that end, one may employ the pertaining EM wave field representations (e.g., ([24], Section 28.12)). For example, the TD scattered far-field amplitudes, {E s;∞ , H s;∞ }, defined via ( [24], Section 26.12) as |r| → ∞, with ξ = r/|r| = cos(φ s ) sin(θ s )i x + sin(φ s ) sin(θ s )i y + cos(θ s )i z being the unit vector in the direction of observation, can be approximately evaluated from where k x (t) and k y (t) directly result from the TD solution (24) via Equations (20) and (21). The corresponding magnetic-type far-field amplitudes simply follow from H s;∞ (ξ, t) = Y 0 ξ × E s;∞ (ξ, t). It is further worth noting that Equations (34)-(36) are, as a matter of fact, exact in the directions normal to the plane of aperture, θ s = 0 and θ s = π, and that the scattered wave fields in z < 0 are equal to the total fields transmitted through the aperture (cf. Equation (1)). For the sake of validation, the z-component of the far-field electric-type amplitude at φ s = π/4 and θ s = 5π/6 was evaluated via Equation (36) using the CdH-MoM solution (24) and with the help of a "far-field probe" in CST Microwave Studio . Figure 5 shows the resulting pulse shapes. Again, the results correlate well with ERR 6.2%.

Conclusions
Using the EM reciprocity theorem of the time-convolution type and the CdH technique, we introduced a fundamentally new solution to the transient EM-field penetration through a relatively small aperture of a rectangular shape. It was demonstrated that for the postulated piecewise-linear space-time distribution of the equivalent magnetic-current surface density induced in the aperture, the pertaining TD admittance array can be derived analytically in terms of elementary functions only. The presented TD solution can be hence viewed as an exact "weak" solution of the EM aperture problem. Subsequently, the CdH-MoM solution was implemented and successfully verified using a commercial 3D EM computational tool. It turns out that the CdH-MoM solution is easy to implement and introduces huge computational savings of several orders of magnitude with respect to general-purpose, TD differential numerical approaches.

Informed Consent Statement: Not applicable.
Data Availability Statement: The presented data are available upon request from the author.

Conflicts of Interest:
The author declares no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Time Domain Admittance Array
Since the xand y-components of the induced magnetic-current surface density are in a small rectangular aperture mutually uncoupled, the TD admittance array has only non-zero diagonal terms, viz.
The elements of the admittance array can be for the piecewise-linear space-time distribution of the induced magnetic-current density (see Equations (5) and (6) with (20) and (21)) obtained analytically in the following form: and Y yy (t) is given by a similar TD expression that can be obtained upon replacing x with y (and vice versa) in Equation (A2). The TD function, Ψ(x, y, t), follows as the inverse of the complex-slowness integral (cf. ( [28], (G.3))): for {s ∈ R; s > 0}, {x ∈ R; x = 0}, {y ∈ R; y = 0}, and recall that γ(κ, σ) = (1/c 2 0 − κ 2 − σ 2 ) 1/2 . Furthermore, K 0 and S 0 are integration paths that run along the imaginary axes in the complex-slowness κ-and σ-planes, respectively, and that are indented to the right around their origins with small semi-circular arcs of vanishingly small radii (cf. ( [28], Figure G.1)). The inversion ofΨ(x, y, s) can be carried out using the CdH procedure as closely described in ( [28], Appendix G). First, the integrand in the integral with respect to σ is continued analytically into the complex σplane, while keeping Re(γ) ≥ 0. This implies the horizontal branch cuts extending along {Im(σ) = 0, Ω 0 (κ) < |Re(σ)| < ∞}. Consequently, using Jordan's lemma and Cauchy's theorem [24] (p. 1054), path S 0 is replaced with the loop around the branch cut, a parametric form of which is σ(u) = −uΩ 0 (κ)sgn(y) ± i0 for all {1 ≤ u < ∞} with sgn(y) = |y|/y. In addition, the contribution from the simple pole at σ = 0 must be accounted for when y ≥ 0. In the integral around the branch cut, we introduce the new variable of integration, u, while the integration around the pole is readily carried out analytically. The thustransformed inner integral with respect to σ is next substituted back in Equation (A3). Proceeding further in a similar manner with the integration in the complex κ-plane, we, after some lengthy, yet straightforward algebra, arrive at integral expressions that can be evaluated analytically. Transforming finally the result of integration to the TD via standard tables ( [33], Section 29), we end up with the TD original of Equation (A3) in the following form ( [28], (cf. Equations (G. 24 +2 v 2 y 2 x 2 y 2 − 3 + x 2 y 2 3 x 2 y 2 − 2 + 3 + 3 4 v x 2 1 (v 2 /y 2 −1) 3/2 v 2 y 2 + x 2 y 2 − 1 . (A5) The calculation of the convolution integral in Equation (A4) is computationally the most exacting task. An efficient way to remedy the issue is the recursive-convolution technique ( [28], Appendix H). Alternatively, one may apply integration by parts and write where ∂ −n v denotes the n-th integration with respect to v. Thanks to the relatively simple form of f (x, y, v) (see Equation (A5)), the required integrals can be readily found analytically, thus enabling exact and fast computation of the TD admittance array elements via Equations (A1), (A2), and (A4).