A New Perspective on the Reactive Electromagnetic Energies and Q factors of Antennas

It is required to calculate the stored reactive energy of an antenna in order to evaluate its Q factor. Although it has been investigated for a long time, some issues still need further clarification. The main difficulty involved is that the reactive energy of an antenna tends to become infinitely large when integrating the conventionally defined energy density in the whole space outside a small sphere containing the antenna. The reactive energy is usually made to be bounded by subtracting an additional term associated with the radiation fields. However, there exists no well-accepted accurate definition for this additional term that is valid for all cases. By re-checking the definition of reactive energies, an alternative formulation is proposed in this paper which can separate the reactive energy and the radiation energy explicitly based on source-potentials. The clearly defined reactive energy is bounded without necessity to subtract that additional term, and the resultant formulae are easy to implement.


I. INTRODUCTION
The stored reactive electromagnetic energy of an antenna is a very important quantity. For example, it can be used to evaluate the Q factor of the antenna and predict its bandwidth. The reactive energy has been investigated by many researchers, and the calculation methods proposed so far can be roughly divided into two categories: (1) methods in early stage based on spherical mode expansion technique [1]- [3], and (2) methods based on determining the fields using computational electromagnetic methods or simulators [4]- [7]. In 1948, Chu in his paper [8] discussed the radiation problem associated with electrically small antennas, and derived ladder type equivalent circuits for TM n0 /TE n0 spherical waves, by which the upper bound on their Q factors can be predicted. The reactive energy only includes those stored in the reactive elements in the equivalent circuit, hence is bounded and can be accurately evaluated. Collin and Rothschild [1] calculated the reactive energies strictly with fields obtained The associate editor coordinating the review of this manuscript and approving it for publication was Su Yan .
using mode decomposition method [9], [10], where the reactive energies of spherical modes and cylindrical modes are obtained by directly integrating the term ε 0 4 E · E * and µ 0 4 H · H * in the whole space outside a sphere with a small radius. Since the integration is infinite, the energy density associated with the radiation fields has to be subtracted from the integrand. Fante [2] extended the results of Collin, and McLean re-examined the case of small antennas and calculated the Q factors of TM 10 and TE 10 mode [3]. For small antennas, spherical mode expansion solution for reactive energy is a good approximation, and can provide satisfactory upper bound for Q factors. It has been extended to analyzing antennas with larger sizes [11], where the radiation fields by a current distribution are expanded with spherical modes. However, it is not efficient because many modes may need to be taken into account for antennas with large size and complicated structures. Furthermore, the fields inside the sphere enclosing the antenna cannot be addressed accurately. Therefore, it is more natural to use numerical methods to calculate the reactive energies, as has been investigated by many researchers [5], [12]- [14]. Basically, a typical numerical procedure to calculate the reactive energies can be carried out in two steps. The first step is to determine all currents involved in the system, such as the induced currents on metal surfaces, the polarization currents in dielectrics, the magnetization currents in magnetic materials, and the currents in lossy media. This can be done using common methods that have been developed in the computational electromagnetic community, such as methods based on electric field integral equations (EFIEs) [15]- [18], volume integral equations (VIEs) [19], [20], and methods based on equivalence principles [21]- [24]. The second step is to calculate the radiation electromagnetic fields from these current sources, and then the stored reactive energies by integrating the energy densities in the whole space. However, as pointed out by many researchers, the stored reactive energy obtained in this way is infinitely large if the conventionally defined electric energy density and magnetic energy density are used. Subtracting an additional term of energy density associated with the radiation fields from the integrand has become a common strategy. Unfortunately, this strategy is not a rigorous solution to reactive electromagnetic energy, and at least three drawbacks need to be addressed. (i) The additional energy density is ultimately an ambiguous concept without a rigorous definition. It is usually determined based on spherical waves with their wave centers located at the origin. Accordingly, their results are coordinate-dependent [4], [25]. (ii) It might be acceptable for small antennas, but may become less accurate for large antennas, in which the propagation pattern is quite different from spherical waves, especially in the region near the radiators. (iii) Theoretically, if we divide the electric field into two parts, for example, E = E rad + E react , with E rad accounting for radiation electric energy and E react for the stored reactive electric energy, then the total electric energy density should include three terms, namely, 0.5ε 0 E rad · E rad , ε 0 E rad · E react , and 0.5ε 0 E react · E react . Obviously, subtracting 0.5ε 0 E rad · E rad from the total energy may not necessarily result in the reactive electric energy because we cannot prove that the two parts of the electric fields are orthogonal to ensure that ε 0 E rad · E react will vanish. This inference is true also for magnetic energy.
Yaghjian and Best have adopted this method to calculate the Q factor of antennas [4]. In order to give a clear explanation, the core idea of their formulation is revisited, but only the radiation problems in free space is considered so that we need not take into account of the effect of materials and can focus on the key issue. In this case, (A.8) in [4], which is the starting point of Yaghjian-Best formulation, can be simplified as where I 0 is the excitation current at the feeding port of the antenna. The primes stand for derivatives w.r.t ω, and X 0 is the input reactance at the feeding port when it is tuned with a series positive inductance or capacitance. V is the spherical domain with radius r. Equation (1) is described as reactive theorem by Rhodes in [26] (the case in free space). The term V Re B · H * + D * · E dV is related to the conventionally defined electromagnetic energy. In Yaghjian-Best formulation (1) is transformed into where ε and η are respectively the permittivity and intrinsic impedance in free space, and F is the far electric field. The first two terms in RHS of (2) are considered as the reactive electromagnetic energy in Yaghjian-Best formulation, which is denoted as Vandenbosch [25] proposed a set of formulae for calculating the reactive energies, which are expressed in closed form of integrations with respect to the current densities in the antenna structure. Again, consider the free space situation, the equation (27) in [25], which is the starting equation of Vandenbosch formulation, can be simplified as Recall that |I 0 | 2 X 0 = −Im V s E · J * dV , and assume J = 0, then (4) is the same as (1). Therefore, the theory bases for Vandenbosch formulation and Yaghjian-Best formulation are almost the same, with a minor difference that the postulation of J = 0 in Vandenbosch formulation is a little bit stronger than the postulation of I 0 = 0 in Yaghjian-Best formulation. The main difference lies in the way to calculate the reactive energies. In Yaghjian-Best formulation, the reactive energies are calculated with W F . In Vandenbosch formulation, −0.5W rad,G is used as the additional term in [25] to replace the term associated with the radiation power. With this modification, the reactive energy can be directly computed with a set of closed-form expressions that are coordinateindependent. Gustafsson and Jonsson evaluated W F analytically to get [27], where W van = W m vac + W e vac + W rad,G is the total reactive energy in Vandenbosch formulation, and W F 2 is a VOLUME 8, 2020 coordinate-dependent term. If the origin shifts within a small sphere containing the antenna, the variation of W F 2 is small.
Yaghjian-Best formulation and Vandenbosch formulation have attracted many attentions from researchers [28]- [33]. They have been successfully applied to analysis and optimization of small antennas [34]- [39], and can be applied for highly dispersive lossy media [40], [41]. The Vandenbosch formulation has been extended to time domain [42], [43]. However, it is reported that the Vandenbosch formulation can produce negative values of stored energy for electrically large structures [38]. Furthermore, the formulation in time domain may give results that are a little bit different from those obtained with the formulation in frequency domain [43].
There are other methods to calculate the reactive energies [44]. A comprehensive comparison of them can be found in [32]. However, for those methods that require to evaluate the reactive energies, W F and W van are the most popular choices. Unfortunately, it lacks rigorous proof to show either W van or W F is the correct expression for the reactive energy. Therefore, the issue has not yet been solved completely. There is still no unique definition of stored energy that is widely accepted in literatures and valid in all cases [45].
A careful re-examination on this issue reveals that the difficulty involved in reactive energy can be traced back to an old classical problem: for a given time-varying current distribution J ( r, t) in domain V a , how to determine the reactive electromagnetic energy stored in the whole free space. Note that for static fields, the electric energy can be computed in different ways. The most commonly used ones are the following two approaches It is easy to check that the second equation in (7) does not hold true for electro-dynamic fields. Only one integral in (7) is correct for expressing the electric energy in time varying fields. Although V ∞ 0.5 E · Dd r 1 is commonly chosen to be the proper one, there exists no rigorous proving or experimental validation. Feynman considered that it was so chosen probably because that the simplest choice is usually the best choice [46]. The above discussed difficulty in calculating the reactive electromagnetic energy of antennas may root in this well known paradox in electromagnetics, which has unfortunately not yet been solved satisfactorily. Based on this observation, the paradox and the basic concept of the reactive energy are revisited in this paper, and a new formulation for calculating the reactive energies of radiators is derived. The reasons for selecting source-potential combinations as the appropriate expressions for reactive electromagnetic energies are provided in Section II. The detailed deduction of an alternative formulation is described in Section III, and validated with a canonical example in Section IV. The expressions for calculating reactive energies and Q factors of antennas are shown in Section V, numerical examples are provided in Section VI, with a brief summary in Section VII.

II. ALTERNATIVE EXPRESSIONS FOR REACTIVE ENERGIES
In order to give a clear description on the formulation, the introducing of the energy densities for static fields are reexamined at first and the derivations are described in detail even though some of them are quite fundamental.
Consider a static charge ρ ( r) existing in a bounded region V a in free space. A popular method to obtain the total energy associated with the charge is to assume that all charges are moved piece by piece from infinite to their current positions. According to energy conservation law, it can be deduced that the total electrostatic energy of the whole system is equal to the work done to the charges, which is derived to be where φ ( r) is the scalar electric potential with its zero reference point located at infinity, and E = −∇φ for static fields. Applying Gauss' Law,∇ · D = ρ, (8) can be cast into For the sake of simplicity, hereafter, S ∞ and V ∞ are used to denote the spherical surface and the space when r → ∞.
Since lim r→∞ D ·r ∼ O 1 r 2 , lim r→∞ φ ∼ O 1 r , the first term at the RHS of (9) approaches zero, and it seems natural to define the electric energy density as The magnetic energy associated with steady-state current distribution J ( r, t), r ∈ V a , can be expressed by where A ( r) is the vector magnetic potential relating to the magnetic flux density B with B = ∇ × A. The zero reference point of A ( r) also locates at infinity. (11) can again be derived based on energy conservation law by considering the process of exciting a loop with its current increases from zero to a certain value i 0 . Applying Ampere's Law in static case, J = ∇ × H , (11) can be transformed to Since lim r→∞ H × A ·r ∼ O 1 r 3 for magneto-static fields, the surface integration in (12) also vanishes. Hence, the conventional magnetic energy density is then defined as However, for time-varying electric fields, E = −∇φ − ∂ A ∂t, the electric energy (8) has to be changed tõ The upper script '∼' is used to indicate time varying variables. The surface integral still vanishes as φ D ·r ∼ O 1 r 3 . We have therefore In time varying case, the Ampere's Law must include the displacement current.
Since H × A ·r ∼ O 1 r 2 , the surface integral in (16) is not zero but a bounded value. The term of the surface integral is the flux passing through S ∞ and can be interpreted as the energy stored beyond the surface S ∞ , which is related to the radiation energy. Apparently, it is reasonable to consider the volume integral in the RHS of (16) as the stored reactive energy associated with the current source, which generally does not equal to V ∞ 0.5 B · H d r 1 for time varying fields.
It is quite clear that, for time varying fields, we have to decide which expression is correct for the electromagnetic energies, or even worse, both of them are not the strictly correct expressions. Contrast to the conventional choice, it is proposed in this paper to select the source-potential combinations as the correct expressions for electromagnetic energies in time varying situations, as some other researchers have proposed [47], [48]. For the sake of convenience, only electric energy is taken into account in the following discussions to support this choice.
Firstly, the choice of source-potential expression for reactive energy is bounded if the sources are bounded in region V a . With this choice, the corresponding reactive electric energy can be calculated with volume integral in terms of fields and potentials as Secondly, there is no direct and rigorous proof or experimental validation to support that V ∞ 0.5 E · Dd r 1 is valid for time varying fields. Seen from (9) and (17), it is more natural to consider V ∞ 0.5 E · Dd r 1 as a secondary quantity relating to 0.5 V a ρ ( r 1 ) φ ( r 1 )d r 1 . Moreover, V ∞ 0.5 E · Dd r 1 can be taken as a special case of the modified electric energy (17) for static fields with ∂ A ∂t · D = 0, which also to some extent justifies the new choice.
Thirdly, as pointed out previously, for time varying fields V ∞ 0.5 E · Dd r 1 is infinitely large. Consider the following scenario. The electric energy associated with a static charge in a bounded region is finite in the whole space outside a small sphere containing the charge. However, with the conventional definition of the reactive electric energy, it will abruptly become infinitely large when the charge begins to vary with time, no matter how slowly it varies. This is quite unnatural, and no method has been found to remedy this issue satisfactorily. Obviously, there is no such difficulty with the source-potential formulation.
Fourthly, the most possible evidence to justify 0.5 E · D as the correct expression of electric energy density for time varying fields may come from the Poynting theorem, which in free space can be expressed as where S is the Poynting vector, and the problem region V a is enclosed by surface S a . The Poynting vector is often regarded as the power flux density. The Poynting relation describes the power balance of the system. Its common interpretation states that the sum of the power flowing out of the surface S a and the power dissipated in the region should equal to the decreasing rate of the total energy in the region. Therefore, it seems quite natural to define the integral at the RHS of (18) as the total electromagnetic energy and the integrand as the energy density. However, the interpretation of Poynting theorem has always been controversial [49], [50]. As a matter of fact, Poynting theorem basically describes the balance between the varying rates of the energies in the system instead of the stored energies themselves. There is no doubt that Poynting equation itself is correct because it is derived directly from the Maxwell equation. However, when we trace back to Poynting's work [51], there was no proof to show that 0.5 E · D is exactly the correct expression for electric energy density of time varying electric fields.
The above deduction can also be applied to magnetic fields. According to (16), the magnetic energy in the new formulation can be calculated with volume integral in terms of fields and potentials over the whole space as For static fields, ∂ D ∂t = 0, (19) yields the conventional expression for reactive magnetic energy. VOLUME 8, 2020 Finally, a case-study with the Hertzian dipole provides a very important support to the new formulation. The analytical results of the electric and magnetic energies obtained using the new formulation exactly agree with those by using Chu's equivalent circuit model, without need to introduce any additional terms like that suggested in [3].

III. REACTIVE ENERGIES AND Q factors
For the sake of convenience, we denotẽ The stored reactive electric energy can be computed by integratingw e over the whole space, It can be checked from (15) thatW E (t) =W ρ (t). So the electric energy can also be calculated with integration over the source distribution region. Similarly, the stored reactive magnetic energy can be computed by integratingw m over the whole space, Making use of (16), it can also be calculated with integration over the source distribution region subtracting a surface integral, that is, For pulse sources, the surface integral in (24) is zero since the fields will never reach S ∞ during a limited time period. For harmonic waves, it will be shown later that the time average of the surface integral is also zero. Therefore, in most situations we haveW M (t) =W J (t), without necessity to evaluate the surface integral at S ∞ . For sinusoidal time varying electromagnetic fields with time dependence of exp (jωt), we havẽ For the sake of simplicity, the same symbols are used for phasors in the expressions. It can be verified that Poynting relation can be written as where the term jω A · D * 4 cancels out.
The time-averaged stored reactive energies of an antenna in free space can then be calculated with The real part of the surface integral in (29) is found to be zero (See Appendix A), so it can be simplified as The radiation power is calculated in terms of far fields, Making use of the Poynting theorem in free space, it can also be calculated in terms of the current source and electric field as With the new formulation for stored reactive energies, the Q factor of a radiator is defined by Similar to the method proposed by Vandenbosch, (33) can be evaluated at a single frequency alone.
Yaghjian and Best have proposed in [4] two Q factors, namely, Q z and Q FBW . Both of them are defined with the antenna being tuned at ω 0 with a positive series inductance or capacitance. Q z is derived from the derivative of an input impedance of the antenna, Although requiring the derivatives, it is significant that (34) has avoided evaluating the stored energy. Q FBW is calculated with the fractional bandwidth, where s is the VSWR used to define the fractional bandwidth FBW s at the tuned angular frequency ω 0 . The bandwidth for a given VSWR at a tuned ω 0 can be determined by searching the reflection coefficient based on the input impedance in both sides around ω 0 . Q FBW depends on the choice of VSWR and it can be seen that Q FBW ≈ Q Z when s → 1.
It is worthwhile to note that Yaghjian-Best formulation is specialized for antennas being tuned at ω 0 with a positive series inductance or capacitance. A more general way is to use (67) in [25] to calculate Q z . The resultant Q factor is independent on the tuning structure, so it is possible to design a wideband matching network much better than a single series inductance or capacitance.

IV. VALIDATION WITH HERTZIAN DIPOLE
In order to compare the proposed formulation with conventional methods, the stored energy and Q factor of a Hertzian dipole is analyzed, which has been calculated by Mclean in [3]. It's known that a Hertzian dipole generates TM 10 spherical mode in free space (denoted with TM 01 in [3]). The fields of a Hertzian dipole are symmetrical about the z-axis, the components of which can be readily derived as below The amplitude of the Hertzian dipole is assumed to be j4π k so that these expressions are exactly the same as those shown in [3]. However, the vector potential given in [3] is not appropriate because it is difficult to find a scalar potential to satisfy the Lorentz Gauge. A proper vector potential can be obtained directly from the dipole current as follows, The terms relating to reactive energy densities are −Im As can be checked straightforwardly that the integration of w e and w m in the space outside a small sphere with radius a is infinite due to the contribution of the 1 r 2 terms. These terms are canceled in [3] by subtracting the energy density associated with the radiation fields, which is With the new definition, the integrand of the stored reactive electric energy and magnetic energy are derived to bẽ w e = µ 4 sin 2 θ 1 k 4 r 6 +4 cos 2 θ 1 2k 2 r 4 + 1 k 4 r 6 (44) Obviously, there is no 1 r 2 term. Integrating them in the space outside a sphere with radius a, the total reactive energies can be obtained The radiation power is found to be Hence, the Q factor for the TM 10 mode is The Q factor is exactly the same as that obtained by Mclean [3]. It can be further verified that, by multiplying a proper scale factor, the reactive electric energy and the reactive magnetic energy are also exactly the same as those stored in the capacitor and the inductor in the equivalent circuit proposed by Chu [8]. Note that the propagating wave of this point source behaves much like a spherical wave. Therefore, subtracting the energy density associated with the radiation power happens to give good results for the total reactive energy in this case. Fig. 1(a). An excitation current J ex exists on the antenna port S p . There is a dielectric with permittivity ε 1 and permeability µ in region V d , together with a PEC conductor with surface S c . The radiation problem can be solved based on equivalence principle. Assume that there is an induced surface current J s on S C , a polarization current J pol in region V d . These equivalent sources are then all placed in free space and are used to account for the effect of the conductor and the dielectric, as shown in Fig.1 (b). Denote the electric field generated by the excitation current as input field to the conductor and the dielectric, which is denoted by

Consider a typical radiation problem shown in
The operator L X ( r 1 ); r 1 ∈ is defined as where g ( r, r 1 ) = e −jkR (4π R) is the scalar Green's function and k is the wave number in free space. R = | r − r 1 |. The tangential component of the electric field vanishes on the PEC surface, so we have the electric field integral equation, In the dielectric, the total electric field includes two parts, where the polarization current relates to the total electric field, Inserting (53) cos kr 12 r 12 × d r 2 d r 1 (54) where r 12 = | r 1 − r 2 |. The integration region is the combination of all source domains, which is denoted by V = V d ∪ S C ∪ S P . The total current is denoted by J = J pol ∪ J S ∪ J ex . Note thatW E andW M are respectively only the first part of W e vac andW m vac given in Vandenbosch formulation. The Poynting theorem in this case can be expressed as where the electric field includes contributions from all sources, From (56), we have the power balance equation, P in = P rad + P lc + P ld (59) Therefore, the Q factor of the antenna can be denoted as All Q factors are defined in the form of (33), with the radiation power therein being replaced by P in , P rad , P lc and P ld , respectively. For PEC conductors, P lc = 0 since the tangential component of the total electric field on the surface is zero. The current on the PEC surface contributes to the storage of reactive energies, but does not contribute to the radiation power directly. However, the surface current will generate electric field on the antenna, hence, indirectly affect the radiation power.
For dielectrics, The power loss is P ld = −0.5ωε 1 E · E * . The imaginary part of 0.5 E · J * pol represents reactive energy stored in the dielectric, which has been already addressed in the reactive energies related to the corresponding polarization current. The polarization currents also influence the radiation power by affecting the electric field distribution at the antenna port.

VI. NUMERICAL EXAMPLES
Six examples are analyzed to show the difference between the new formulation and the method proposed in [4] and [25]. In order to get rid of possible ambiguity, the feeding structures are clearly illustrated in the examples. All feeding currents have amplitude of I 0 = 1A and distribute uniformly along the reference direction on the feeding patch.
In the examples, the plates in the antennas are assumed to be perfectly electrically conducting (PEC), and are meshed into conforming triangular meshes with meshing software. In order to guarantee the numerical accuracy, the average edge length of the mesh is about 1/10 of the wavelength of the highest frequency considered. All surface currents are expanded with Rao-Wilton-Glisson (RWG) basis functions [15] and are calculated by solving the EFIE (49) on the antennas. Standard Galerkin scheme is adopted in the computations, and singularity-subtraction techniques have been used to accurately evaluate the entries of the impedance matrix.
Basically, five Q factors are compared. Q po , Q Z , and Q FBW are respectively calculated with (33), (34) and (35). Q F , Q van are calculated in the way similar to Q po , but with the reactive energies replaced by W F and W van , respectively. Unless specified differently, the VSWR used to determine the bandwidth is 1.5 and the corresponding Q factor is denoted by Q FBW −1.5 . Q F is calculated with the origin located at the center of the structure.
A fractional discrepancy between the Q factors obtained with the new formulation and that with Vandenbosch formulation is defined as, A mathematical uniform surface current ring is analyzed firstly. Similar to that in [52], the current is not associated with a real metal plate, but is a purely feeding source distributed over a circular ring with radius of 15mm and width of 0.5mm. The phase of the current is e −j2ϕ , which varies linearly along the circle. The calculated energies and Q factors are shown in Fig.2(a). It can be seen that all the reactive electric energies (WE van , WE F ) and reactive magnetic energies (WM van , WM F ) calculated with Vandenbosch formulation and formulae of W F are negative near 29GHz, where the reactive energies (W ρϕ , W AJ ) calculated with source-potential formulation are all positive. The Q factors are shown in Fig.2(b). Q F is calculated with the origin located at the center of the ring, and almost coincides with Q van in this example. Owing to the negative energies, Q van and Q F are negative near 29GHz. It can be noted that there is a jump in Q FBW −1.5 at 25.6GHz, which is caused by a neighboring local resonance. Q po is always positive. We have re-calculated the energies for a larger ring with a radius of 150mm and strip width of 0.5mm. The phase of the current also varies linearly along the circle. The numerical results are not shown here. From the frequency of 1GHz to 15GHz, the reactive energies (WE van,F ) have 10 negative regions, however, the reactive energies (W ρϕ , W AJ ) are positive in the whole frequency range.
Next, a plate dipole is analyzed. It consists of two PEC plates with size of 500mm × 2mm. A feeding patch with size of 2mm × 2mm is placed between the two plates. The surface current on the antenna is calculated by solving the corresponding EFIF with Galerkin method. The results of the input resistance R and reactance X are shown in Fig.3(a), and the calculated Q factors are plotted in Fig.3(b). It can be seen that all the five Q factors are roughly close to each other in the examined frequency band. The relative discrepancy Q is less than 16% in this example. Although not plotted in the figure, it has been checked that the effect of the origin is quite small in this case due to the symmetrical property of the fields. It has also been checked that when the VSWR is set to be 2.0, 1.5, and 1.05, the corresponding Q factors, denoted respectively by Q FBW −2.0 , Q FBW −1.5 , and Q FBW −1.05 , are all close to Q Z (only Q FBW −1.5 is plotted.) The third example is a square loop antenna with edge length of 30mm, as shown in Fig. 4(a), excluding the PEC ground plate. The width of the PEC strip is 0.5mm. A 0.5mm × 0.5mm feeding patch is put at the center of one segment of the square. The results of Q factors are plotted in Fig. 4(b). The five Q factors reveal similar behavior but with larger discrepancies than that in the second example. In the fourth example, a PEC plate with size of 35mm × 35mm is placed under the loop antenna, as shown with dot lines in Fig. 4(a). When the ground plane is 2mm away, the calculated Q factors are plotted in Fig.5(a). It can be seen that Q po , Q van and Q F are very close. Q Z and Q FBW agree well with each other, and also roughly agree with the other three Q factors. However, Q FBW has spurs near the natural resonances, where the input impedance varies sharply, as can be seen from Fig. 5(b). With the increasing of h, the distance between the loop and ground plate, the Q factors gradually decease and approach to those of the loop without ground, as shown in Fig. 5(c). The relative discrepancies are shown  in Fig. 5(d). The discrepancy tends to become smaller when the ground plate approaches the loop.
The fifth example is a PEC strip grid antenna consisting of uniform strips with width of 0.8mm. The antenna is symmetrical, but the feeding point is usually offset to get optimal matching performance. The calculated Q factors are shown in Fig.6, from which it can be seen that Q F and Q van almost coincide, and Q po is close to them, with discrepancy illustrated in Fig. 6(b). Due to the effect of local resonances, Q FBW-1.5 does not agree with Q Z at frequency range of about 8∼9 GHz. The influence of the choice of VSWR on Q FBW is illustrated in Fig.6(c), where Q FBW-1.05 almost coincides with Q Z , but Q FBW-1.5 and Q FBW-2.0 may have jumps in them. Although larger VSWR causes larger bandwidth, the corresponding Q FBW is not necessary to be smaller because its definition includes the VSWR (s), as can be seen from Fig.6(c).
The effect of the choice of origin is also illustrated in Fig.6(c). The origin is chosen to be located uniformly at 20 points on the yellow dot-lined circle shown in Fig.6(a). The resulted Q F s are shown with grey lines in Fig.6(c). Although the effect is bounded, it is not always neglectable.
The last example is a Vivaldi antenna with structure shown in Fig.7(a). The opening rate is R a = 0.0458 mm [53]. The calculated input resistance and reactance are shown in Fig.7(b). Since its resistance is near 50ohm in a very wide frequency range, it is easy to realize a wide band antenna with further optimization. The Q factors are shown in Fig.7(c). Similar to the results of the strip grid antenna, Q po , Q van and Q F (with origin at the center of the feeding patch) are close to each other with a relative discrepancy at most 18%. However, Q FBW and Q Z are obviously away from them. The effect of the choice of VSWR is illustrated in Fig.7(d).
Since the Vivaldi antenna is basically a wideband antenna, the variation of Q FBW −2.0 is not quite smooth because many local resonances may fall into the pass band at a specified tuned frequency. By changing the origin to 20 points located on the yellow dot-lined circle in Fig. 7(a), the resultant Q F s are shown with the grey line in Fig. 7(d).
All these examples show that Q po computed using the new formulation are always close to Q van with relatively small discrepancies. This is not strange since the main body of the expressions for the reactive energies of the two formulations is identical. Q F usually depends on the choice of the origin, and the variation range is proportional to the offset distance. For antennas with symmetrical radiation fields, the effect may be very small.
The numerical results also show that Q Z and Q FBW calculated using Yaghjian-Best formulation are also close to Q po and Q van for narrow band radiators, such as the dipoles, but may differ largely from them for wide band radiators, as shown in the last two examples. This implies that the input impedance of a radiator may not necessarily completely describe the field-based energy processes or distributions, which is also observed in [54]. The definitions of Q Z and Q FBW are different from those directly defined based on reactive energies. So their results are not necessarily to be equal. They may have different advantages in different applications.

VII. CONCLUSION
A new formulation is proposed to calculate the stored reactive energy of a radiator based on sources and potentials for harmonic fields. Similar to the conventional definitions for the electromagnetic reactive energies, the new definition also lacks rigorous mathematical proof. However, some arguments are provided to support the new formulation. Particularly, the case study of Hertzian dipole has verified that the stored reactive energies evaluated by the new formulation exactly agree with the classical results obtained based on Chu's circuit model. On the other hand, by analyzing the energy balance relationship in a radiation problem, the new formulation can provide an insight into the composition of the energies radiated by an antenna. The expressions for the reactive energies are quite simple and easy to implement. Q factors can be evaluated at a single frequency point. At some situations in which the conventional methods may take negative reactive energies while the new formulation still provides positive results.
This work focuses on the reactive electromagnetic energies of radiators in free space or radiation problems that can be handled with equivalent models in free space. Further investigations are needed on systems consisting of isotropic, lossy or time varying media, as well as the mathematical aspect of the new formulation.

APPENDIX A
In order to show that the real part of S ∞ H × A * · d S in (29) approaches zero, we have to use the asymptotic behavior of the scalar Green's function in 3D free space. Since in this case the integration is performed for time harmonic fields, we can write where g 1,2 = g r, r 1,2 , J 1,2 = J r 1,2 . Furthermore, we have the asymptotic expressions lim r→∞ g 1,2 = e −jkr 4π r e −j k· r 1,2 (63) lim r→∞ ∇g 1,2 ∼ −jke −jkr 4π r e −j k· r 1,2r ≈ −jkg 1,2r (64) Making use of the relationship that lim r→∞ ∇g * 2 · J 1 J * 2 · ∇g 1 ≈ lim r→∞ ∇g 1 · J 1 J * 2 · ∇g * 2 (65) and V a ∇ 1,2 · g 1,2 J 1,2 d r 1,2 = 0 (66) It can be deduced that S ∞ V a V a ∇g 1 · J * 2 J 1 ·rg * 2 d r 2 d r 1 dS = 1 jk S ∞ V a V a ∇ 1 g 1 · J 1 J * 2 · ∇ 2 g * 2 d r 2 d r 1 dS Using (63) and performing the integration directly gives Hence we have (70), as shown at the top of the page. where it is understood that the two currents in the two-fold integral are the same one but at different locations in the same region V a . By exchanging r 1 and r 2 , together with the order of the inner-fold and the outer-fold integration, it is straightforward to check that which leads to Re S ∞ H × A * · d S = 0.

APPENDIX B
Inserting the scalar potential into (14) and applying the current continuity equation yields With the same analogy to prove (71), it is straightforward to check that the first two-fold integral at the RHS of the last equation is pure real, while the second one is pure imaginary. Hence (54) is proved. Inserting the vector potential into (16) and applying the current continuity equation yields (55) can be obtained in the same way to derive (54). It can also be verified that, if adopting the identity (A.1) in [25], (54) and (55) can be derived respectively from the direct integration of the energy densities (17)  His current research interests include computational electromagnetics and its application in scattering and radiation problems. VOLUME 8, 2020