Approximation of High Intensity Radiated Field by Direct Current Injection using matrix methods based on Characteristic Mode Analysis

. This contribution discusses the approximation of radiated by conducted immunity tests by the example of High Intensity Radiated Field (HIRF) and Direct Current Injection (DCI) based on a surface current analysis. For this purpose, Characteristic Mode Analysis (CMA) is applied to provide basis functions for a surface current expansion in Characteristic Modes. Via a matrix-based basis transformation algo-rithm involving Characteristic Mode data of both HIRF and DCI test setups, suitable DCI surface currents are derived. The approximation of HIRF surface currents by the computed DCI surface currents is analyzed for exemplary DUTs over a broad frequency range. Within this frequency range, those DCI frequencies leading to an optimal approximation of the HIRF current are determined. Concerning practical issues in DCI testing, the inﬂuence of DCI adapter parameters on the surface current approximation is elucidated. The numerical results show that DCI can approximate HIRF at low frequencies largely independent from the DCI adapter setting, whereas at high frequencies an approximation is difﬁ-cult to realize.


Introduction
In Electromagnetic Compatibility (EMC) testing of aircrafts, radiated immunity tests at high field levels are customarily referred to as High Intensity Radiated Field (HIRF) tests (EUROCAE, 2010).This procedure requires an irradiation of the Device Under Test (DUT) from multiple angles with different polarizations, making HIRF rather time-consuming.It is moreover expensive to implement at low test frequencies, because large antennas and powerful amplifiers are necessary.To remedy, inter alia, these deficiencies, the conducted test method of Direct Current Injection (DCI) was introduced in avionics (Carter and Willis, 1991).It employs the injection of currents onto the conducting surface of a DUT using voltage or current sources as injection adapters, whereas additional termination adapters realize a current return path via a ground plane (Leat, 2007).Instead of a ground plane as return conductor, a coaxial return structure enclosing the DUT is likewise possible (Pérez et al., 2019).Despite the introduction of DCI several decades ago, its correlation to other EMC test methods like HIRF is still subject to contemporary research, see Rothenhäusler et al. (2019Rothenhäusler et al. ( , 2023)), Wang et al. (2020), Ückerseifer et al. (2019, 2021, 2023).
In order to investigate whether DCI can generate EMC test conditions similar to HIRF, the surface equivalence principle is applied.In its generalized form, it states that the EM-sources in a domain can be replaced by equivalent electric and magnetic surface currents (Huygens sources) at the boundary of a closed surface separating an outer domain I from an inner domain II (Harrington, 1961).The fields in the outer (E I , H I ) or inner (E II , H II ) domain are considered equivalent, whereas the fields in the other domain are usually set to zero (Love's equivalence principle) as most simple solution of Maxwell's equations satisfying the boundary conditions on the surface.In the special case of a scattering problem in which the closed surface is identical to the DUT's perfectly conducting surface, the equivalent surface current corresponds to the electric surface current J S excited by the EM-sources (Jin, 2010), see Fig. 1.Consequently, if the surface current of HIRF caused by irradiation and the injected current of DCI have similar structures, both methods can be considered equivalent from an EMC point of view with respect to the outer domain I.
As mathematical instrument to investigate surface currents, Characteristic Mode Analysis (CMA) (Harrington and Mautz, 1971;Cabedo-Fabres et al., 2007;Chen and Wang, 2015) is applied.A CMA yields orthogonal basis functions, called Characteristic Modes, in which surface currents can be expanded.Characteristic Modes only depend on geometric and material properties of the test setup and are thus different for HIRF and DCI because of distinct boundary conditions introduced by the DCI adapters (Rothenhäusler and Gronwald, 2017).
Based on Characteristic Mode data of HIRF and DCI setups, a basis transformation of the HIRF surface current from HIRF modes to DCI modes is performed, by which expansion coefficients for a DCI surface current expansion follow.Its approximation to the HIRF current is analyzed by means of a relative error measure for complex-valued vector quantities.Unlike in Ückerseifer et al. (2019), the basis transformation algorithm employs only vector-matrix operations and is thus easier to automate for a broad frequency range and arbitrary DUTs including spherical bodies.Apart from the discussed surface current transformations, practical issues concerning the influence of DCI adapter impedance as well as injection and termination adapter positions on the surface current approximation is discussed.
Section 2 starts with a mathematical overview over CMA and its specific modal quantities.They constitute the foundation for a basis transformation algorithm subsequently demonstrated yielding DCI surface currents as approximation to HIRF surface currents.In order to apply this algorithm, Sect. 3 introduces HIRF and DCI test setups for three different canonical DUTs, whose modal properties are determined by a CMA.Finally, Sect. 4 presents computed DCI currents for these DUTs and compares them to HIRF currents in a broad frequency spectrum for different physical configurations of the DCI setup.

Mathematical background
The mathematical notation throughout this paper indicates complex quantities with an underscore.In contrast to scalar quantities, vectors and matrices are symbolized by bold variables.

Characteristic Mode Analysis
Characteristic Mode Analysis for perfect electric conductors (PEC) is commonly based on the Electric Field Integral Equation (EFIE) operator Z = R+jX (impedance matrix) set up in the Method of Moments (Harrington, 1967).In a CMA, the generalized eigenvalue problem at a frequency f with real eigenvalues λ n (f ) ∈ (−∞, ∞) and eigenvectors J n (f ) ∈ R called Characteristic Modes, both quantities of order n ∈ N, is solved.The eigenvalues are commonly expressed in a more convenient form with limited value range as the modal significance to a surface current distribution.As an example, Fig. 2 shows amplitude-normalized Characteristic Modes J H n , J D n of a straight wire acting as DUT in HIRF and DCI configuration.The boundary conditions for the HIRF modes J H n manifest as vanishing currents at both wire ends, whereas for the DCI modes J D n non-vanishing boundary currents occur because of DCI adapters attached to each wire end.Similar differences between HIRF and DCI modes occur for more complex bodies as well.
As the J n are linearly independent, they form a basis in which the surface current J HIRF on a DUT resulting from a HIRF test can be expanded.Assuming a test frequency f H , an expansion in M Characteristic Modes is set up according to Here, J H n denote Characteristic Modes of the HIRF setup and the expansion coefficients β n are called modal weighting coefficients.Both quantities are provided by numerical software tools.To investigate whether a DCI surface current J DCI can theoretically approximate a HIRF current distribution, the latter needs to be expandable in Characteristic Modes of the DCI setup, J D n , and corresponding modal weighting coefficients γ n like where f D is the frequency, at which a CMA of the DCI setup is carried out.Appropriate γ n for an adequate approximation of J HIRF by J D can be obtained via a basis transformation algorithm described in Knabner and Barth (2018).It is implemented in MATLAB (MathWorks, Version R2021b) using CMA data generated by the integral equation solver of CST Microwave Studio (Dassault Systèmes, Version 2021.05) for DUTs with N DUT mesh elements.Due to its enhanced mesh capabilities, all meshes are created in the numerical software FEKO (Altair Engineering, Version 2022.0.1) with a standard level of discretization and imported to CST in NAS-TRAN format.The DCI mesh with N DCI = N DUT + N ADAP mesh elements serves as reference, from which the HIRF mesh with N HIRF = N DUT elements is obtained by omitting the adapter mesh elements N ADAP .This ensures equal positions of all mesh elements for the HIRF and DCI setups and thus enables an element-wise surface current comparison of both test methods.

One-dimensional DUT
A basis transformation matrix T with elements t nm for the transformation of one-dimensional current distributions is deduced by expanding each J H n as linear combination of all The approximation sign is due to finite summation over M Characteristic Modes as well as the different boundary conditions for HIRF and DCI modes, see Fig. 2. In Eq. ( 5), the J D n only include the N DUT mesh elements belonging to the DUT, since the basis transformation can only be applied to these current elements.
In order to transform Eq. ( 5) from an index representation to a matrix one, both sets of Characteristic Modes J H n (i), J D n (i), evaluated at the ith mesh element, are collected as column vectors in two matrices yielding At a given simulation frequency f , these matrices contain the most significant modes, i.e. s 1 (f ) > s 2 (f ) > . . .> s M (f ), which are found by deactivating mode tracking (Ludick et al., 2014) in the CMA solver of CST.With Eq. ( 6), the equivalent matrix form of Eq. ( 5) reads as The desired transformation matrix T is hence and (V D ) −1 is in general, i.e. for non-quadratic matrices V D (M = N), evaluated as Moore-Penrose inverse (pseudoinverse), see Campbell and Meyer (2009).
With the help of T, the vectorized DCI weighting coefficients γ = [γ 1 , . .., γ M ] T in Eq. ( 4) follow from the HIRF weighting coefficients β = [β 1 , . .., β M ] T via a coordinate transformation.It is deduced, in a first step, by inserting Eq. ( 5) in Eq. ( 3) with the result and a subsequent comparison of coefficients with Eq. (4) as or in matrix-vector notation Having determined γ , expansion (4) can be set up and compared to the current J HIRF .

Two-or three-dimensional DUT
In two or three dimensions the HIRF weighting coefficients are transformed into separate DCI weighting coefficients γ x , γ y , γ z for each cartesian current component.For this purpose, the matrices (6) are split into V H x , V H y , V H z and V D x , V D y , V D z containing only one single cartesian component of J H n , J D n .This yields transformation matrices and weighting coefficients by which the cartesian current components are expanded as They can be combined to T , which equals the seeked approximation (4) of the surface current J HIRF .
Figure 3 summarizes the overall basis transformation algorithm implemented in MATLAB based on Characteristic Mode data of HIRF and DCI setups provided by CST.Having determined the DCI expansion J D , its goodness of approximation to the reference current J HIRF can be analyzed in a last step.

Test setups and DUT properties
The outlined basis transformation is applied to three canonical DUTs (straight wire, rectangular plate, circular cylinder) modelled as PEC in HIRF and DCI configuration for f H = f D , see Fig. 4.Each HIRF setup consists of an infinite PEC ground plane placed 1 m beneath the DUT.It is exposed to a plane wave with an electric field amplitude of |E i | = 1 V m −1 and a polarization angle of 45 • at two exemplary frequencies f H 1 = 50 MHz and f H 2 = 100 MHz for which the DUT appears electrically large (resonance region), such that a sufficient number of Characteristic Modes is excited (Blevins, 2006).The corresponding DCI setup vertically connects the DUT to this ground plane by one injection and one termination adapter, either having a characteristic impedance of Z inj = Z term = 50 .
The number of Characteristic Modes N (f ) present on all three test objects in free space is illustrated in Fig. 5 for modes fulfilling s n ≥ 0.1.Due to low-frequency instabilities inherent to RWG-based (Rao-Wilton-Glisson) EFIE formulations of CMA (Qian and Chew, 2008), (Dai et al., 2017), especially for those containing line sources like DCI excitations (Hofmann et al., 2022), the lower frequency is chosen sufficiently high in order to obtain a substantial number of modes.Based on the polynomially fitted number of modes in Fig. 5, the mode density D(f ) = dN (f )/df in Fig. 6 is calculated.It shows a constant mode density for the wire, a linearly rising one for the plate and a parabolic one in accordance with results for ordinary and generalized eigenvalue problems in Blevins (2006), Courant and Hilbert (1924).

Current distributions for HIRF and DCI
Figure 7 illustrates the currents J HIRF , J D on a straight wire (length l = 5 m, radius r = 1 mm) with N DUT = 58 mesh segments for M = 5.Both curves show a good agreement at f H 1 , but at f H 2 slight amplitude deviations and, more importantly, non-vanishing boundary currents due to current injec-  tion by the DCI adapters occur.Both effects can be mitigated by increasing M, as can be seen in Fig. 8.
The surface currents J HIRF , J D on a rectangular plate (5 m × 2 m) with N DUT = 1464 mesh triangles are plotted in Fig. 9 (top view).Even though the number of modes is increased to M = 10, certain amplitude deviations at f H 2 remain visible.This is explainable with an increased mode density for two-dimensional DUTs, as emphasized in Fig. 6.
Finally, the surface currents J HIRF , J D on a cylinder (l = 5 m, r = 1 m) meshed with N DUT = 5692 triangles are shown in Fig. 10 for M = 10.In this scenario, the approximation of the HIRF current at both frequencies is insufficient owing to the even higher mode density of volumetric bodies, again conveyed by Fig. 6.Since it is known that inner resonances of volumetric bodies deteriorate the accuracy of EFIE-based CMA formulations because of an ill-conditioned  impedance matrix (Dai et al., 2015), the more stable Combined Field Integral Equation (CFIE) formulation of CMA is used for the cylinder (Mautz and Harrington, 1977).

Broadband comparison of HIRF and DCI
The approximation of the HIRF current J HIRF by the DCI expansion J D can be quantified over a broad frequency range for all three DUTs.A measure for the relative error of two complex-valued vector quantities (Mittra, 2014) whose result is shown in Fig. 11 from 1 to 200 MHz in steps of 1 MHz and different M for wire, plate and cylinder.In accordance with the previous results, it is visible that more modes are required at higher frequencies due to a higher mode density for geometric complex DUTs.In contrast, a relatively small error of approximately 20 % can be observed for the wire and plate up to frequencies of around 100 MHz.
Steep ascents and descents in ε between neighbouring frequencies indicate changes in the sorting of Characteristic Modes deriving from deactivated mode tracking in CST (Ludick et al., 2014).In Eq. ( 15), all DUT mesh elements are evaluated, although for EMC testing an approximation of interference maxima is most important.A criterion to include only mesh elements of relevant current amplitude is their energy contribution to the total current distribution.Since the PEC boundary condition on a boundary described by its unit normal vector n enforces identical amplitudes of surface current and magnetic field H , each current element can be assigned a magnetic energy (Pozar, 2011) within a local mesh element volume V i close to the DUT surface.In Eq. ( 17), the complex conjugate H * (i) of H (i) is indicated by an asterisk.Calculating all W i and sorting them in descending order (W 1 > W 2 > . . .> W N DUT ), these elements of highest energy can be summed up and related to the total energy present in the accumulated volume V = N DUT i=1 V i according to until at a certain number of considered elements N lim a predefined percentual limit ξ lim is exceeded.For several values ξ lim , Fig. 12 shows the relative error ε with N DUT in Eq. ( 15) replaced by N lim regarding the cylinder.Apparently, mesh elements with high energy dominate the overall surface current approximation, as adding low-energy elements does not significantly affect the error at most frequencies.Figure 13, in addition, reveals the required number of elements N lim referenced to all N DUT elements to reach the limits ξ lim .It is remarkable, comparing to Fig. 12, that a very limited number of mesh elements is responsible for the surface current approximation.
Hitherto, the frequencies f H and f D were identical.Since the DCI adapters change the electrical properties of the DUT, an investigation, whether CMA frequencies f D = f H lead to a better approximation of the HIRF current at f H 1 and f H 2 ,   As additional degree of freedom, the termination adapter impedance Z term in Fig. 4 can be altered from a matched condition (Z term = 50 ) to the extrema of a short circuit (Z term = 1µ ) and an open circuit (Z term = 1 M ), see Fig. 16.In general, the curves show a good correspondence despite the significant impedance range.Consequently, the impedance can be chosen to fulfil practical criteria like impedance matching to achieve sufficient test levels.
A further study involves the placement of injection and termination adapters (Z inj = Z term = 50 ), which are indicated by numbered dots in Fig. 17.In addition to the excitation scenario 1 → 2 of Fig. 4, the cases 3 → 4 and 1 → 5 are ana-  lyzed with the outcome in Fig. 18.All adapter configurations result in similar errors, such that for practical considerations those adapter positions enabling an independent excitation of Characteristic Modes should be selected (Peitzmeier andManteuffel, 2018, 2019).

Conclusions
The presented results show that, in theory, radiated immunity tests can be approximated well by conducted immunity tests at relatively low test frequencies for simple DUTs.As a means of analysis, CMA and the outlined basis transformation algorithm provide a systematic approach to determine suitable DCI excitations for a large variety of DCI configurations.Especially at low frequencies with a small mode density a good agreement between the HIRF currents and their approximation in DCI modes is achieved.The approximation is primarily determined by a restricted amount of relevant mesh elements.When the mode density increases simultaneously to the frequency, a sound approximation is difficult to realize for complex DUTs and a large number of modes is required.Moreover, the algorithm is computationally expensive, since for each physical change of the DCI setup, a new CMA in the entire frequency range is necessary.

Figure 1 .
Figure 1.Surface equivalence principle for scattering by a perfect electric conductor (PEC) (a) Original problem.(b) Equivalent problem.

Figure 2 .
Figure 2. Normalized Characteristic Modes J H n , J D n of a straight wire (l = 5 m) in HIRF and DCI configuration for n = 1, 2, 3.

Figure 3 .
Figure 3. Determination of suitable DCI surface currents from Characteristic Mode data of HIRF and DCI setups.

Figure 4 .
Figure 4. Canonical test objects in HIRF and DCI configuration.

Figure 5 .
Figure 5. Unfitted and fitted number of Characteristic Modes N(f ) for wire, plate and cylinder.

Figure 6 .
Figure 6.Mode density D(f ) for wire, plate and cylinder.

Figure 11 .
Figure 11.Relative error ε for f H = f D .

Figure 16 .
Figure 16.Relative error ε for three termination impedances Z term and M = 10.

Figure 18 .
Figure 18.Relative error ε for different adapter positions and M = 10.