Padé approximation for time-domain plane wave reflected from a lossy earth

To deal with the behaviour of wire structures subjected to electromagnetic waves with the finite-difference time-domain (FDTD) method, it is necessary to establish the external electromagnetic excitation in time-domain. For a homogeneous lossy half-space problem, the Fresnel reflection coefficients for arbitrary polarization plane wave are usually multiplied by the pulse spectrum in the frequency domain and the result is transformed into the time-domain. In this paper, an accurate method based on the Padé approximation theory in conjunction with the partial fraction expansion is used, and leads to expressions of the field above a homogeneous lossy earth that can analytically be transformed into the time-domain. The influence of the common degree of the numerator and denominator of the Padé approximant on the accuracy of the proposed method is discussed.


Introduction
In many electromagnetic coupling to transmission lines problems, an accurate knowledge of electromagnetic fields of external excitation is required to calculate induced voltages in the line by using field-to-line coupling equations. When a ground is included as a homogeneous lossy half-space, the incident, reflected and transmitted fields are usually pre-calculated in the frequency domain where it is relatively simple to apply Fresnel coefficients to the incident pulse spectrum. Time-domain fields are typically obtained through the use of the inverse fast Fourier transform (IFFT), but it can be computationally time demanding and prone to aliasing effects. This conventional frequency domain approach has been widely discussed in literature [1][2][3][4][5][6][7][8].
For some instances, especially the analysis of the reflection of an electromagnetic plane wave incident upon a dielectric slab [9], the study of electromagnetic lateral pulses generated by vertical and horizontal electric dipoles [1] (Chapters 13 and 14), and finite-difference time-domain analysis of overhead lines within a circuit simulator [10][11][12], it would be useful to have an expression for the transient reflected field from the boundary between two electrically different half-spaces. To do this, it is appropriate to find an adequate representation of the Fresnel coefficients in the time domain. Several studies on this topic have developed a method or approximations to deal with reflection coefficients directly in the time domain. Barnes and Tesche [13] developed approximate analytical expressions for the transient earth-reflected fields by considering an impulsive plane wave for both vertical and horizontal polarizations. The impulse response for a lossy ground is obtained directly by analytical inversion of the Laplace transform which has given rise to approximations and thus a limitation to grazing angles.
Other methods based on the frequency to time-domain transformation technique generated time-domain approximations of Fresnel coefficients. E J Rothwell et al [14] applied the inverse Fourier transform to the frequency-domain horizontal reflection coefficient, and split its time-domain expression into a pulse term representing the lossless case, and another characterizing the dispersive nature of the conducting half-space. The dispersive term involves an integral of Bessel functions which has been approximated by an infinite sum. A series form for a time-domain vertical Fresnel coefficient, which require numerous Bessel function evaluations, has been introduced in [15]. The procedure is similar to that developed in [14] and it is based on inverse Fourier transform. To achieve reasonable results, a corrective term was introduced to the series initially valid for angles of incidence less than the Brewster angle.
To circumvent the use of the Bessel functions and their evaluations in a repetitive way when calculating the time-domain vertical Fresnel coefficient, S Antonijevic et al [16] introduced a method which is valid in case where the relative electrical permittivity is not too small and the angle of incidence not too large. The technique is based on the decomposition of the frequency-domain vertical Fresnel coefficient into a space-dependent and a space-frequency-dependent (SFD) functions. The SFD function is then derived using Gaver-Stehfest algorithm for numerical inverse Laplace transform and a simple convergence acceleration algorithm.
In this paper, the proposed technique for achieving the impulsive ground response is more efficient and easier than those mentioned above. It does not require the recursive evaluation of Bessel function or the use of the numerical inverse Laplace transform. Time-domain reflection coefficients are obtained using the Padé approximation in conjunction with the partial fraction expansion.

Plane wave above ground
Consider a transient electromagnetic plane wave whose propagation direction S is −a r as shown in the spherical coordinate system (r, θ, f) (see figure 1 ). The wave is impinging on a finitely conducting half-space and propagates in free space with propagation velocity The region (y<0) is designated as a lossy homogeneous earth, characterized by a permittivity ò g =ò rg ò 0 , a permeability μ g =μ 0 , and a conductivity σ g . ò 0 , μ 0 are the electrical constants of the free space (y>0). The polarization of the electric field vector is described in terms of the unit vectors in the spherical coordinate system a θ and a f . γ is the polarization angle. The general expression of the electric field above the air-earth interface at some point, identified by the position vector r=xa x +ya y +za z , is where E 0 (s) is the complex amplitude, and s is the complex frequency Laplace transform variable. s x , s y , s z are the components of the unit propagation vector S along x−, y−, and z−axis s a cos sin 3 The time quantities t x , t y and t z are defined by The components of the polarization unit vector e=e x a x +e y a y +e z a z in the represent the components relating to vertical and horizontal polarization waves respectively. The reflected electric field is related to the incident field trough the Fresnel reflection coefficients, witch are determined by matching the appropriate electromagnetic field components across the air-earth interface [2]. For an arbitrary polarization, this reflected electric field is expressed as follows where Γ v and Γ h are the Fresnel reflection coefficients for the electric field parallel and perpendicular to the plane of incidence respectively. n denotes the complex refractive index of the assumed homogeneous earth. It is expressed as By rearranging the expression (6) so as to highlight the contribution for each polarization, and applying the inverse Laplace transform, we obtain the following expression of the transient reflected field t t t e t u t s t s t s t e P P , , 9 where  is the convolution product, δ(.) is the impulse function and u(.) the unit step function. γ v , γ h , and t e 0 ( ) denote inverse Laplace transforms of Γ v (s, θ), Γ h (s, θ), and E 0 (s) respectively. The vectors P v and P h witch characterize the vertical and horizontal polarizations are defined by The expression (9) shows that the transient reflection field is written as a convolution of the time-domain reflection coefficients and the waveform used for the excitation. It is therefore obvious that these reflection coefficients represent the impulse response for the conducting half space. For the remainder of this paper, we will focus on evaluating this impulse response efficiently and accurately. We will also discuss the comparison between the proposed method and the direct calculation from the spectrum by the standard inverse Fourier transform.
3. Padé approximation for impulse response evaluation 3.1. Transient response The impulse response in (9) requires time-domain reflection coefficients (RC). In order to obtain the simplest and accurate expressions for (RC), we use the Padé approximation [17]. The idea of this approximation is to expand Γ v (s, θ) and Γ h (s, θ) as a ratio of two power series and determining both the numerator and denominator coefficients using the coefficients of their Taylor series expansion in complex s-domain.
To begin, we expand Γ v,h (s, θ) into a Taylor series at (s=s 0 ) which can be written in the following compact form is the nth moment of Taylor expansion. This expansion is the fundamental starting point of any analysis using Padé approximation. A Padé approximation s, where a 0 , K, a L , b 1 , K, b M are the unknowns. b 0 is taken to be equal to 1, which eliminate an arbitrary multiplicative factor in (13). In our approach, the numerator and denominator degrees are chosen equal (L=M) to provide a good estimate of reflection coefficients in the whole frequency range. Matching η(s, θ), in terms of its moments, to the rational approximation given in (13) By matching the two sides of (14), the coefficients in the numerator and the denominator of the rational approximation ( The numerator coefficients can be found by equating the powers of (s−s 0 ) from s s 0 0 - Note that (13) can also formulated by partial fraction decomposition as where p î are approximate poles which can be easily obtained by applying Python [18] numpy.roots module on polynomial B(s). ĉ is equal to a M /b M , and k î are the residues. In order to solve for the residues k î , the calculated poles are used to perform the following limits By factoring the polynomial B s ( ) as a product of linear polynomials, it is easily to show that The form in (17) is more convenient for an inversion of Laplace transform into the time domain. The time domain impulse response of the conducting earth can then be given in closed form in terms of the poles and residues as x are the constant term, poles and corresponding residues for each wave polarization, i.e. vertical or horizontal.

Computation of Taylor expansion moments
It should be emphasized that the efficiency of the proposed method for time-domain reflection coefficients is essentially based on the accurate evaluation of the coefficients of Taylor series expansion. To perform the computation of the moments of Taylor expansion given in (11), we have adopted an algorithm, which was originally developed by Hosono [19], that involves neither manipulations nor numerical approximations to the derivatives. Thus, it does not suffer from the usual instabilities associated with the numerical differentiation formulas for calculating higher order derivatives; also, it does not involve the use of symbolic manipulation packages. The method relies on evaluating the kth expansion coefficient, for a given function h(z) holomorphic on the unit disk, by the following integral along the circumference on the unit circle (C), Substituting the exponential function k z exp( ) in (22) with the periodic function e a (k z) with respect to z defined as e k z a j k z a j n exp 2 . 1 0.5 , 23 a n and by closing the contour of (22) on the right half-plane z, the residue theorem leads to where a is a positive number which determines the degree of approximation and .
[ ] I represents the imaginary part. As mentioned in [19], the expression (24) is valid if h z h z * * = ( ) ( )( * is the complex conjugate sign) on the unit disk, and h(0)=0. In this work, the formula (24) can not be used directly for the reflection coefficients. Indeed, neither the restriction to the unit disk nor the null value in zero are satisfied. To deal with this case, let's consider where R is appropriate positive value which will be specified in the numerical results section. Using formula (24) for s, h q( ), we get the moments η [n] of the reflection coefficients at s=s 0 a n R a n j k n exp .
However, as will be shown in the next section, the accuracy of the expression (26) declines rapidly for the higher moments. To improve the accuracy, we introduced a second adjustment parameter b according to what was discussed in [19] so that the expression (26) becomes a b n n R a n b j k n exp .

Results and discussion
In this section, some numerical results will be presented to show the accuracy and efficiency of the proposed method for calculating the ground impulse response (GIR). To begin with, we will first check the accuracy of the equations (26) and (27).
In order to directly obtain the real-valued (GIR), the Taylor moments are computed in the vicinity of an optimum real value s 0 =10 6 , which provides a good precision for all the parameters involved in the reflection coefficients, namely the angle of incidence, the conductivity and the relative permittivity. Figure 2 shows the comparison, in terms of relative error with respect to Python results, between the relationships (26) and (27) where η [n],Py is the nth moment of the Taylor expansion of Γ v,h (s, θ) around s=s 0 =10 6 which is calculated using the Sympy.Series Python module [19]. The lossy earth is characterized by a relative permittivity ò rg =12. and a conductivity σ g =0.001. The incidence angle is fixed at θ=15°. The values of a, b and R in (27) are set at 7.5, 0.78 and 2s 0 respectively. These values are chosen so as to minimize the error in the whole frequency range for all electrical and geometric parameters. From figure 2, it appears that the formula (26) is not quite satisfactory, specially for the moments of higher orders, while equation (27) cancels this error by providing increasing accuracy as the order of the moment increases. Moreover, we note that the error is almost identical for both vertical and horizontal Fresnel coefficients. The first ten moments for the reflection coefficients used to calculate the error, for the same parameters mentioned above, are reported in table 1. Table 2 shows the poles and residues associated with the Fresnel reflection coefficients for both vertical and horizontal polarizations. It should be emphasized that the Padé approximation can result in the right half of the s-plane approximating poles for a given order which lead to a divergence of the impulse response. This is the case, for example, of the coefficient v G which exhibits a positive pole for M=10, ò rg =12, σ g =0.001S/m, and θ=15°(see table 2). In this work, positive poles are ignored. The results of the proposed approximation (dashed lines) for M=10 are in good agreement with the exact results (solid lines), and this for all the parameters involved in the reflection coefficients, as shown in figures 3, 4, and 6.
The curves in figure 3 represent the Fourier spectrum Γ ξ=v, h (−jω) and its approximation jw Gx ( ) as a function of frequency ω for different values of incidence angle θ. The imperfectly conducting earth is assumed to have a conductivity of σ g =0.001S/m and a relative dielectric constant ò rg =12. In figure 4, the angle of incidence is set at θ=45°and the dielectric constant ò rg =12. The real and imaginary parts of the exact and approximate spectra of Fresnel reflection coefficients are plotted for different values of the soil conductivity. From the examination of figures 3 and 4, it appears that the frequency behaviours of the vertical and horizontal reflection coefficients are similar, the behaviour of one at low frequency is exactly the same as that of the other at high frequency, and vice versa. The value of s 0 must therefore be judiciously chosen to allow a good precision for both v G and h G , and this for all the parameters on which the reflection coefficients depend. Figure 5 shows that the choice of the relatively low value s 0 (s 0 =10 5 ) leads to a loss of accuracy for high frequencies and conversely high values of s 0 (s 0 =5×10 7 ) distort the estimation of reflection coefficients for low frequencies. By testing several values of s 0 between 10 5 and 5×10 7 for different conductivities, permittivities and incidence angles, we retained the value s 0 =10 6 which produces more precision for all parameters. The effect of the dielectric constant on j figure 6 for θ=45°and σ g =0.001S/m. It should be noted that the proposed method presents a slight loss of accuracy, especially for the imaginary part of j v h , w G -x= ( ), for very low ground conductivities and grazing angles of incidence as shown in figures 3, 4 and 6. Figure 7 shows the effect of the denominator degree M of the Padé approximant on the accuracy of v G . It can be seen that the value M=10 is largely sufficient to accurately reproduce the frequency behaviour of the Fresnel reflection coefficients.
To deal with overhead lines subjected to a plane wave with an arbitrary polarization using finite-difference time-domain method and Agrawal coupling model [20], vertical and horizontal electric field components are required. A simple manipulation of (2) and (6) can be performed to yield the y-and z-directed components of the reflected electric field at a position denoted by (x, y, z) as In order to see the effects of the approximations used in developing the reflection coefficient expressions on the z-and y-directed transient reflected electric field components, we need to specify the time behaviour of the external source. The frequently used form of an electromagnetic pulse is In this work, the parameters assume the following values: E 0 =56.6 kV/m, α=4.086×10 6 (1/s), β=1.565×10 8 (1/s). This two-exponential shape is often used to simulate high-altitude electromagnetic pulses (HEMP). The Laplace transform of the double exponential is where M ξ=v,h is the effective order obtained from M by subtracting the number of positive poles, and t s t s t s t x x y y z z It should be noted that even though the excitation waveform is complex, the proposed method for reflection coefficients, which uses a sum of exponential functions, offers the possibility of using recursive convolution [21] to determine the reflected field for a specified waveform.
The reflected field with an arbitrary polarization, which is expressed using the original reflection coefficients, is calculated by converting the Fourier spectra  waveforms corresponding to the FFT and the direct calculation through the analytic formulation are superimposed for M=10.
Finally, it is interesting to compare the proposed method with those published in [13] and [16], in order to illustrate its effectiveness. To do this, consider the double exponential time-domain z-directed reflected field defined by the following convolution integral [13]:  ( ) I denote the real and imaginary parts.  where Here N T denotes the number of discrete samples of t, g q ¢ x ( )and Δt the time-domain resolution. The recursive formula (36) minimises the required mathematical operations involved with the numerical evaluation of the convolution integral. Indeed, only two additional past values of t, ), need to be retained at each step of the recurrence. As indicated in (36), one is used to update the function Ψ i and the other one to determine the value of e p t 1 , . It should be emphasized that recursive convolution algorithms, which have a much better numerical performance, can be used especially in the case where either of the functions to be convolved is a summation of exponential functions [21,22] . This is particularly the case when one has to compute the reflected electric field using the RC functions published in [13,16] by adopting an incident waveform that is not expressed as a summation of exponentials. Considering the form of RC approximations in (21), it is then obvious that the proposed method has the advantage of using recursive convolution algorithms for evaluating the reflected field.  The impulse responses of (21), (25) in [13] and (19) in [16] for both vertical and horizontal polarizations have been convolved numerically, by means of the above proposed recursive convolution algorithm, with the double exponential waveform defined in [13,16] by E 0 =52.5 kV/m, α=4×10 s 6 1 -, and β=4.76×10 8 s −1 . The resulting transient reflected electric fields are plotted in figures 9 and 10. The curves designated by 'FFT' correspond to the reflected electric field obtained by using inverse fast Fourier transform of Γ ξ (−jω) E 0 (−jω), (ξ=v, h). The curves marked 'approximate' result from only two-term approximation as defined in [13]. It should be noted that for the reflected field defined by (35), the propagation time shift t 0 has been omitted. Figure 9 shows the transient plane wave reflected from a conducting half-space of electrical parameters ò rg =10 and 0.01S m g s = .
For the incidence angle of θ=45°, the results of references [13] and [16] are exactly reproduced and compared to the current method for M=10 and M=4 as depicted in figure 9(a). It can be seen that for this angle of incidence which is within the range of angles that match the validity condition of the approximations in [13] and [16], namely cos 1 24.62 rg rg , the approximate solutions ( [13], equation (25)) and ( [16], equation (19), N=6) are in good agreement with the proposed solution for M=10. Similar responses are plotted in figure 9(b) for an incidence angle of θ=15°. As expected, a decrease in accuracy is observed for both methods published in [13] and [16] for grazing angle. In contrast, the current method for M=10 provides good accuracy. Figure 10(a) shows the reflected electric field for horizontal polarization which has been calculated using ( [13], equation (25)) and its two-term approximation. The result is in accordance with the one given in [13]. In addition, the present method provides good accuracy for grazing angles and with an order M even lower (M=4) when the conductivity is sufficiently large (here σ g =0.5 S/m ), as illustrated in figure 9(b).

Conclusion
This paper has discussed a new approach for calculating the electromagnetic field reflected by a homogeneous imperfectly conducting ground and emanating from an arbitrary polarization plane wave. The half-space impulse response has been described by a simple summation of real-argument exponentials which are obtained through the Padé approximation in conjunction with partial fraction expansion. In order to calculate the Taylor moments of the vertical and horizontal polarization reflection coefficients, we have adopted a technique that does not suffer from the usual instabilities associated with higher order numerical differentiation formulas and does not use symbolic manipulations. This technique is based essentially on the values of the reflection coefficients. The proposed approach has provided an analytical representation of the reflected transient field associated with a double exponential incidence frequently used in electromagnetic wave coupling studies. In addition, it makes possible the use of recursive convolution in the case where the incident waveform is quite complex. Finally, the proposed approach is essentially based on the two parameters s 0 and M which respectively represent the point around which Taylor series expansions of the reflection coefficients have been made and the order of the Padé approximant. From the study carried out on all the parameters, namely the conductivity, the permittivity and the incidence angle, it appears that the values s 0 =10 6 and M=10 provide good accuracy.  Finally, it is easy to show that (A.5) satisfies the following recursive relationship: