Instruction FACTA UNIVERSITATIS Series: Electronics and Energetics Vol. 27, N o 2, June 2014, pp. 221 - 234 DOI: 10.2298/FUEE1402221D CAUSAL MODELS OF ELECTRICALLY LARGE AND LOSSY DIELECTRIC BODIES  Antonije Djordjević 1 , Dragan Olćan 1 , Mirjana Stojilović 2 , Miloš Pavlović 3 , Branko Kolundžija 1 , Dejan Tošić 1 1 University of Belgrade – School of Electrical Engineering, Belgrade, Serbia 2 University of Applied Sciences Western Switzerland, Yverdon-les-Bains, Switzerland 3 WIPL-D d.o.o., Belgrade, Serbia Abstract. This paper presents a novel formula for the complex permittivity of lossy dielectrics, which is valid in a broad frequency range and is ensuring a causal impulse response in the time domain. The application of this formula is demonstrated through the analysis of wet soil, where the coefficients of the formula are tuned to match the measured data from the literature. Additionally, an analytical expression for the impulse response of the relative permittivity is derived. The influence of the frequency dependence of the complex permittivity on the causality of responses is illustrated through the analysis of 1-D, 2-D, and 3-D electromagnetic systems. Being the most complex, the 3-D system is also used as a test bed for comparing the computational limitations of two commercially available solvers, CST and WIPL-D. Key words: causal response, complex permittivity, wet soil, impulse response, CST, WIPL-D. 1. INTRODUCTION Contemporary software tools for electromagnetic (EM) simulation can be efficiently used for modeling and analysis of various complex systems. Yet, when a system comprises electrically large but highly-detailed objects filled with lossy dielectrics, as is usually the case for simulations at microwave and millimeter-wave frequencies, software limitations can easily be reached. An example of a large and complex system is a human body, which is often modeled when considering body-area networks. In order to find the response of such a system in the time domain, there are two general simulation strategies: to apply a time-domain solver, or to apply a frequency-domain solver and then use the inverse discrete Fourier transform.  Received January 23, 2014 Corresponding author: Antonije Djordjević University of Belgrade – School of Electrical Engineering, Bulevar kralja Aleksandra 73, P.O. Box 35-54, 11120 Belgrade, Serbia (e-mail: edjordja@etf.bg.ac.rs) mailto:edjordja@etf.bg.ac.rs 222 A. DJORDJEVIĆ, D. OLĆAN, M. STOJILOVIĆ, M. PAVLOVIĆ, B. KOLUNDŢIJA, D. TOŠIĆ If a frequency-domain solver is used, the system needs to be analyzed at a large number of equispaced frequencies. At the lower end of the frequency range, the electrical dimensions of the body are small compared with the wavelength, so that integral-equation solvers are the preferred choice [1]. At the upper end of the frequency range, various asymptotic techniques may be used, provided that the shape of the body is simple (e.g., a sphere). If, however, the shape of the body is highly-detailed, then an integral-equation solver (IS) is preferred again. However, the computational resources required for IS (the memory and processor requirements) quickly increase with increasing the simulation frequency. That is why the contemporary commercial integral-equation solvers often fail to analyze objects whose overall dimensions exceed several tens or hundreds of wavelengths. In this paper, we will challenge the limits of two available commercial solvers implemented in WIPL-D [2] and CST [3] electromagnetic simulation tools. The time-domain response of any real (physical) system is causal: the response cannot start before the excitation. All time-domain solvers implement a time-stepping procedure and naturally incorporate the causality feature, as in [3]. However, this is not the case with frequency-domain solvers. To ensure causal response, one must not forget to properly model the dielectric relative permittivity. If the dielectric is lossless, the relative permittivity can be independent of frequency. Yet, if the dielectric is lossy, the variations with frequency must be described by an appropriate function of frequency, as discussed in Section 2. Otherwise, a noncausal, and thus unrealistic, response will be obtained. In order to illustrate the causality issues, we use wet soil as an example. The parameters of the soil are evaluated based on measured data available in the literature, presented in Section 3.1. These data are fitted in a very broad frequency range, as described in Section 3.2. The fitting function involves the broadband term from [4]. In the time-domain analysis, a dispersive relative permittivity is described by the corresponding impulse response. Hence, for the broadband term from [4], we reveal in Section 3.3 the corresponding impulse response, which is not available in the literature. In Section 4 we present three examples to demonstrate differences between a causal and a noncausal model of the soil. The examples are sorted by computational complexity. The first example is a one-dimensional (1-D) electromagnetic problem (plane-wave propagation). The second example is a two-dimensional (2-D) problem (a cylindrical dielectric scatterer). Finally, the third example is the most complex one: a three- dimensional (3-D) problem consisting of a large dielectric cube and two dipole antennas. 2. CAUSALITY ISSUES WITH FREQUENCY-DEPENDENT PERMITTIVITY Parameters of lossy media are frequency-dependent. Thus, a linear nonmagnetic medium is characterized by the complex relative permittivity )(j)()j( rrr  , where f 2 is the angular frequency and f is the frequency. Here, the negative of the imaginary part, )( r  , takes into account both conductive and dielectric losses. A physical system is causal: the response cannot occur before the excitation. Consequently, )( r  and r ( )  are not mutually independent. Under certain conditions [5], they are related by the Hilbert transform, or, equivalently, the Kramer- Kronig relations (dispersion relations). Causal Models of Electrically Large and Lossy Dielectric Bodies 223 Time-domain solvers that are capable of dealing with dispersive parameters inherently use causal models of media. Usually, they fit the relative permittivity in the frequency domain using simple terms, most often Debye terms [6], so that the impulse response is obtained analytically using the Fourier transform. The impulse response is afterwards used in convolution integrals in the time domain. Frequency-domain solvers can deal with any kind of frequency dependence of )j( r  , because they are not bound by causality issues. However, if we use the results of the frequency-domain analysis to compute the time-domain response, )j( r  should be such as to provide a causal response. Otherwise, the response in the time domain would have a non-physical behavior. For example, it could start before the excitation, thus violating strict causality [7], or the speed of propagation of electromagnetic fields could exceed the speed of light in a vacuum, violating Einstein’s causality. Although it is often stated in the literature that )( r  can be evaluated from )( r  [8], and vice versa, the required numerical integration is not easy, and sometimes even not doable. The reasons for that lie in singular, highly oscillatory, or even diverging integrands, and in infinite integration limits in the Hilbert transform. Even analytically, the integrals cannot be evaluated in many important practical cases because the integrals are divergent or undefined [5]. As a simple example, let us consider a leaky dielectric characterized by const)( r  (equal to the electrostatic relative permittivity) and by a constant conductivity const)(  , which is independent from r  . The equivalent (complex) permittivity of the material is r r 0 ( j ) j( / )       , so that r 0 ( ) /     . Clearly, if only r  is known, it is impossible to find )( r  without an additional piece of information, and vice versa. In consequence, the data for )j( r  that satisfy causality conditions can be most reliably and easily supplied in terms of an analytic function of the complex frequency s (  js on the imaginary axis). This function, )( r s , cannot have poles in the right half-plane. It can have only simple poles on the imaginary axis, where it must possess conjugate symmetry: )j(*)j( rr  . Hence, )(**)( rr ss  . The function )( r s can be supplied directly by the user, in an analytic form. Alternatively, the user tabulates the frequency-dependent data for )( r  and )( r  , and the solver evaluates an appropriate interpolation formula as in [3]. For direct analysis in the time-domain, the impulse response of r  is needed. It is convolved with the vector )( 0 tE to obtain )(tD . This convolution is the time-domain counterpart of the relation ED 0r  in the frequency domain. 3. COMPLEX RELATIVE PERMITTIVITY OF WET SOIL In order to clearly demonstrate the difference between causal and noncausal models, it is preferable to have a medium with relatively high losses. In that case, the causality issues can be noted even after an EM wave propagates along a short distance. We have selected wet soil as an example of a dispersive medium throughout the remainder of this paper. We have characterized the soil based on experimental data presented in Subsection 3.1. The analytic approximation for )( r s is given in Subsection 3.2. 224 A. DJORDJEVIĆ, D. OLĆAN, M. STOJILOVIĆ, M. PAVLOVIĆ, B. KOLUNDŢIJA, D. TOŠIĆ 3.1. Experimental data on relative permittivity of water and soil In this paper we use two sets of measured relative permittivity values: those of water and those of soil. We combine them to estimate the relative permittivity of wet soil. Measurement results and analytical model for the relative permittivity of water are given in [9]. The measured complex permittivity is fitted by a constant ( r ), two Debye (relaxation) terms, and a frequency-independent conductivity (  ), as 0 2w r2 lw r1 rr 11 )(             sss s . (1) In this model, conductive losses are attributed to  and polarization losses to the two Debye terms. The measured relative permittivity of soil is given in [10], and comprise data for r  , r  , and  . Although it is not sufficiently clear from the report, r and  are not independent; they are related by r 0 /( )    , which can be verified from the numerical results presented in the paper. In other words, two distinct descriptions of losses are used in reference [10]. In the first description, all losses (conductive and polarization) are attributed to r  . In the second description, all losses are attributed to  . In this paper we use results from the middle row of Fig. 36 in [10]. 3.2. Broadband approximation of relative permittivity of wet soil In [4], an approximation for frequency-dependent complex relative permittivity of lossy dielectrics is proposed, covering a very wide frequency range. It uses a logarithmic function, which provides practically constant )( r  , while )( r  slowly decays with frequency. This approximation yields a causal response in the time domain. The complete expression for the relative permittivity is given by o 1 2 12 r j 10ln j j ln ' ')j(          mm , (2) where [rad/s]1101 log m and [rad/s]2102 log m . The first term is the relative permittivity at very high frequencies, the second term is the broadband logarithmic term, and the third term comes from the conductivity, which is assumed to be independent of frequency. For the frequency range where 21  , the real part of the logarithmic term is 10ln ln ' 10ln j j ln ' 10ln j j ln ' Re 2 12 1 2 12 1 2 12                             mmmmmm , (3) and it linearly decays with the logarithm of the frequency, for 2 1 ' ( )m m  per decade. In the same frequency band, the imaginary part of the integral, Causal Models of Electrically Large and Lossy Dielectric Bodies 225 10ln 2' 10ln j j arg ' 10ln j j ln ' Im 12 1 2 12 1 2 12                                   mmmmmm , (4) is practically constant. For angular frequencies below 1  or above 2  , the imaginary part of the logarithmic term tends to zero, while the real part tends to be constant. This logarithmic function can replace several Debye terms in a wide frequency range. The formula (2) is often quoted as “Djordjevic-Sarkar” model, and it has been built into Ansoft [11], Agilent [12], Simberian [13], and other software. We make an approximation of the parameters of wet soil by combining the soil parameters from [10], the logarithmic term from [4], and the approximation for pure water at 25ºC, based on data from [9] and [14]. The measured data for the soil are in the frequency range from 0.1 GHz to 3 GHz, whereas the approximation is valid outside of this frequency range as well. Our approximation for the permittivity of wet soil reads:                          s s mm ps 1 2 12 r r ln 10ln 11)(   0 2w rw w lw rw w 11 1                        ss p s pp , (5) where  11.0p is the relative contribution of water,  11 2 f , where MHz 1 1 f is the lower cutoff frequency of the broadband term,  22 2 f , where THz 100 2 f is the upper cutoff frequency of the broadband term,  [rad/s]1101 log m ,  [rad/s]2102 log m ,  r rd 2 1 ( )m m    is the total variation of the real part of the broadband term, where 61.1 rd  is the slope per decade,  S/m025.0 is the constant conductivity,  1w1w 2 f , where GHz 25 1w f is the location of the first Debye term for water,  2w2w 2 f , where GHz 200 2w f is the location of the second Debye term for water,  5.76 rw  is the total variation of the real part of the permittivity of water, and  065.0 w p is the relative contribution of the second Debye term. 226 A. DJORDJEVIĆ, D. OLĆAN, M. STOJILOVIĆ, M. PAVLOVIĆ, B. KOLUNDŢIJA, D. TOŠIĆ Fig. 1 compares the measured data from [10] with the results obtained from our approximation formula in the frequency range from 0.1 GHz to 3 GHz. The measured data exhibit stochastic errors, but the achieved agreement can be considered to be good. (a) (b) Fig. 1 Comparison of measured and analytically calculated parameters of wet soil: (a) relative permittivity and (b) conductivity. The results show good agreement. 3.3. Impulse response For time-domain solvers, we need the impulse response of )( r s from equation (5). This response can be evaluated by summing the responses for all terms. The first term is a constant (unity), the second term is a broadband term, two Debye terms follow, and the last term has the form s/1 . The impulse responses for all terms, except the broadband term, are elementary and can be found in standard tables of the inverse Laplace transform. However, the response for the broadband term is not available in the literature. Hence, we have evaluated it analytically using the inverse Laplace transform, so that the impulse response corresponding to (5) reads: l 2 r r 2 1 e e ( ) ( ) (1 ) h( ) ( ) ln10 t t t t p t m m t                    lw 2w w rw lw w rw 2w 0 ((1 ) e e ) h( ) h( ) t t p p p t t             , (6) where )(t is the Dirac (delta) function and )(h t is the Heaviside (step) function. 4. EXAMPLES OF (NON)CAUSAL RESPONSE In this section, we present three examples to demonstrate differences in the time- domain response when using a causal and when using a noncausal model of wet soil. The examples are ordered according to the complexity and dimensionality of the analyzed electromagnetic problems, from the simplest to the most complex ones. Causal Models of Electrically Large and Lossy Dielectric Bodies 227 4.1. Plane wave We consider a uniform plane wave that propagates through a homogeneous nonmagnetic medium. This is a one-dimensional (1-D) electromagnetic problem. The excited wave is described by a delta-function. The distance of wave propagation is m 5.0d . We analyze the propagation in the frequency domain at 4096 frequency points, starting from 0, with a step of 1 MHz. Thereafter, we use the inverse discrete Fourier transform to obtain the response in the time domain. We consider two models. First, when the complex relative permittivity is given by (5). Second, when the complex relative permittivity is independent of frequency and equal to 36402j725116 r . .  (which is an estimated mean value of the permittivity of wet soil in the first model). The results are shown in Fig. 2. For the first model, a causal response is obtained. It has a crisp start at 6 ns. For the second model, a noncausal response is obtained. It is characterized by a premature and “lazy” leading edge of the pulse. Fig. 2 Time-domain response when a plane wave is propagating through a causal and a noncausal medium. The causal response has a crisp start, while the noncausal response has an early start and slow leading edge. 4.2. Dielectric cylinder We consider an infinitely long cylinder of a square cross-section, whose side length is 0.5 m. The cross-section of the cylinder is shown as an inset in Fig. 3. The axis of the cylinder coincides with the z-axis of the Cartesian coordinate system. Hence, we deal here with a two-dimensional (2-D) EM system. 228 A. DJORDJEVIĆ, D. OLĆAN, M. STOJILOVIĆ, M. PAVLOVIĆ, B. KOLUNDŢIJA, D. TOŠIĆ Fig. 3 Electric field z-component at the axis of a dielectric cylinder for the cases when the permittivity is frequency-constant (noncausal) and frequency-dependent (causal). A uniform plane electromagnetic wave illuminates the cylinder. The wave propagates along the x-axis, in the opposite direction of the x-axis. The electric-field vector of the wave is a Gaussian pulse in the time domain, defined as 2 0 2 ( ) 2 0 ( ) e t t z t E   E i , (7) where, V/m 1 0 E , ns 3 0 t , ns 1.0 , and z i is the unit vector in the z-direction. These data are for the transversal plane at m 1x . Since the vector of the electric field is parallel to the cylinder axis, the electromagnetic field in this system is a 2-D transversal magnetic-field, usually termed as TM mode. The numerical analysis is performed using the method of moments (MoM) with the surface integral-equation formulation, piecewise constant approximation for electric and magnetic surface currents, PMCHWT formulation, and point-matching testing procedure [1], [15]. As the response, we calculate the z-component of the electric field at the cylinder axis (i.e., at point O at the inset of Fig. 3). The frequency-domain analysis is done from 9.99512 MHz to 10.235 GHz over 1024 equidistant frequency samples. The total number of unknowns increases with the increase of the analysis frequency, but it does not exceed 1600. The total analysis time is 815 s on a desktop computer with Intel i7 CPU and 32 GB of DDR3 RAM. The time-domain response is calculated for the time interval from 0 to 100 ns over 2048 equidistant time samples. The results for constant permittivity in the whole frequency range, 3j16 r  , which is a noncausal model, and for r ( )s given by (5), which is a causal model, are shown in Causal Models of Electrically Large and Lossy Dielectric Bodies 229 Fig. 3 for the first 15 ns. As in Fig. 2, the response obtained by the first model shows a premature beginning of the leading edge. In contrast, the causal model has a clear start, which allows for much more precise timing when evaluating the beginning of the impulse response. 4.3. Dielectric cube As the most resource-demanding problem, we consider the three-dimensional (3-D) system shown in Fig. 4. It consists of a cube, made of a lossy dielectric, and two symmetrical dipoles. The side of the cube is 2c, where we take two values for c: mm 100c for a smaller cube and mm 250c for a larger cube. Inside the cube, at its center (which coincides with the coordinate origin O), one symmetrical dipole (dipole #1) is located. The length of one arm of the dipole is 5 mm (the overall dipole length is 10 mm). The wire radius is 0.1 mm (the diameter is 0.2 mm). Another dipole (dipole #2) is located outside the cube. The arm length of this dipole is 20 mm (40 mm overall) and the wire radius is 0.5 mm (the diameter is 1 mm). The dipoles are mutually parallel, and parallel to the height of the cube. The distance between the dipole centers (feeding points, ports) is 2c. Fig. 4 Two dipole antennas, one of which is inside a lossy dielectric cube. Not drawn to scale. Both cubes are analyzed at 1000 frequencies: 10 MHz, 20 MHz, ..., 10000 MHz (10 GHz) using the program WIPL-D [2]. For each frequency, the impedance and scattering parameters are computed. The nominal impedance for the scattering parameters is 50 Ω. Two models of the dielectric are used: a noncausal model, for which the complex relative permittivity is frequency independent, 3j16 r  (i.e., 16 r  and 3 r  ), and a causal model, for which the complex relative permittivity is evaluated from equation (5). For both cubes and for both dielectric models, the impulse responses for s21 and z21 are evaluated. The results are shown in Figs. 5–8. 230 A. DJORDJEVIĆ, D. OLĆAN, M. STOJILOVIĆ, M. PAVLOVIĆ, B. KOLUNDŢIJA, D. TOŠIĆ WIPL-D simulations are done using WIPL-D Pro version 11, on a desktop computer with Intel CPU core i7 3820 @3.60 GHz, 64 GB of (DDR3) RAM, Nvidia GeForce GTX590, and with Microsoft Windows 7 Pro 64-bit operating system. The EM system is modeled using three symmetry planes in order to maximally reduce the computer resources needed for the analysis. Although such symmetry introduces a parasitic image of the second (larger) dipole, the effect of the parasitic dipole is negligible as it is located far away. First, we analyze the smaller cube ( mm 100c ). In the case of the constant relative permittivity, 3j16r  , the analysis takes 26,931 s (i.e., approximately 7.5 hrs), and the total number of unknowns increases with frequency from 1,202 at the lowest frequency (10 MHz) up to 4,943 at the highest frequency (10 GHz). In the case when the relative permittivity is given by equation (5), the simulation takes 26,172 s (i.e., approximately 7.3 hrs). The total number of unknowns increases with frequency from 1,202 at the lowest frequency up to 4,803 at the highest frequency. Thereafter, we analyze the larger cube ( mm 250c ). The analysis for the constant permittivity lasts 113,753 s (i.e., approximately 31.6 hrs). The number of unknowns is 1,202 at the lowest frequency and rises with frequency up to 15,233 for the highest frequency. The analysis for the relative permittivity given by equation (5) lasts 105,484 s (i.e., 29.3 hrs), while the number of unknowns is in the range from 1,202 to 13,621. Note that in both cases, the analysis of the cube with frequency dependent permittivity lasts slightly less than the analysis with constant permittivity. This is due to the fact that WIPL-D allocates resources by taking into the account the electrical size of the structure at the operating frequency. At higher frequencies, the modulus of the frequency- dependent permittivity is smaller than the modulus of the constant permittivity, thus demanding fewer unknowns. The results for mm 100c obtained by WIPL-D are compared with the results obtained by program CST [3], which uses a time-domain solver (Fig. 5). CST simulations were performed using Microwave Studio software from the CST Studio 2013 package, on a Windows 7 64-bit server equipped with two Intel(R) Xeon(R) CPUs @2GHz and 192 GB RAM. Unfortunately, CST cannot analyze the case for mm 250c for the given hardware configuration. The CST time-domain solver uses the causal model of the relative permittivity. The hexahedral mesh size was about 7 million cells. The cube interior was filled with a frequency dispersive material, defined using the appropriate permittivity values given by equation (5) for each frequency point. The background was filled with a vacuum. The model boundaries were set to “open”. Two symmetry planes (one magnetic-field and one electric-field symmetry plane) were defined to reduce the total computational load. The solver accuracy was set to 80 dB. The excitation signal was of a 10 GHz CST-default Gaussian type. The simulation was run so to ultimately yield 1001 equidistant frequency points in the 0 to 10 GHz frequency range. The total simulation time is approximately 67.5 hours. The time-domain solver in CST evaluates only the scattering parameters. Thereby, only one port is excited, so that in one run of the program the parameters s11 and s21 are evaluated. In order to calculate the impedance parameter z21, the parameter s22 is needed as well. However, the computation of s22 requires another full-time run of CST, which was not performed to avoid the long run of the program. Consequently, only the impulse response for s21 is shown (Fig. 5). Causal Models of Electrically Large and Lossy Dielectric Bodies 231 The simulation for the noncausal response in CST can be performed using a frequency-domain solver only. The simulation parameters are as follows: integral solver (IS), a mesh with 1313 surfaces, a vacuum background, open boundaries, two symmetry planes, solver accuracy 1e3, s-parameters normalized to 50 Ω, 3rd order solver, and the solver type is MoM. The results are also shown in Fig. 5. The total simulation time is approximately 52 hours. The agreement between the results evaluated by WIPL-D and by CST is very good, both for the causal model and the noncausal model. Fig. 5 Impulse response for mm 100c , for s21, and zoom-in (inset), as computed by CST time-domain solver for the causal model, by CST frequency-domain solver for the noncausal model, and by WIPL-D for both the causal and noncausal models. Fig. 6 shows the impulse response for z21 for the smaller cube. Figs. 7 and 8 show the impulse response for the larger cube, for s21 and z21, respectively. These results were computed only by WIPL-D. The noncausal model of the dielectric yields a premature start of the response, which is more visible for z21 than for s21. The explanation is in the shape of the spectrum of these two parameters. The spectrum of the parameter z21 is wider than the spectrum of s21. Hence, the inadequate variations of the permittivity of the noncausal model have influence in a wider frequency range for z21 than for s21. 232 A. DJORDJEVIĆ, D. OLĆAN, M. STOJILOVIĆ, M. PAVLOVIĆ, B. KOLUNDŢIJA, D. TOŠIĆ Fig. 6 Impulse response for mm 100c , for z21, and zoom-in (inset), as computed by WIPL-D for the noncausal and causal models. Fig. 7 Impulse response for mm 250c , for s21, and zoom-in (inset), as computed by WIPL-D for the noncausal and causal models. Causal Models of Electrically Large and Lossy Dielectric Bodies 233 Fig. 8 Impulse response for mm 250c , for z21, and zoom-in (inset), as computed by WIPL-D for the noncausal and causal models. 7. CONCLUSION This paper presents an analytical expression of complex permittivity of wet soil, valid in a broad frequency range, which assures a causal response in the time domain. The parameters of the formula are tuned to fit the measured data for soil and water in a broad range of frequencies. The impulse response, needed for direct analysis in the time domain, is derived, too. The discrepancies between the causal and noncausal responses, and their relations with the complex permittivity of the material, are illustrated through several examples of different dimensionality and complexity. It is shown that in all cases the causal response has a crisp start, while the noncausal response has an early and slow leading edge. Additionally, a model of a 3-D EM system, being the most complex example, is used to test the present-day limits of some commercial EM solvers. Acknowledgement: The paper is a part of the research done within the project TR32005 of the Serbian Ministry of Education, Science, and Technological Development. REFERENCES [1] B. M. Kolundţija and A. R. Djordjević, Electromagnetic Modeling of Composite Metallic and Dielectric Structures, Boston: Artech House, 2002. [2] WIPL-D Pro 3-D. Available: http://www.wipl-d.com/ [3] CST, 3D Electromagnetic Simulation Software. Available: https://www.cst.com/ 234 A. DJORDJEVIĆ, D. OLĆAN, M. STOJILOVIĆ, M. PAVLOVIĆ, B. KOLUNDŢIJA, D. TOŠIĆ [4] A. R. Djordjević, R. M. Biljić, V. D. Likar-Smiljanić, and T. K. Sarkar, “Wideband frequency-domain characterization of FR-4 and time-domain causality”, IEEE Trans. Electromagn. Compat., vol. 43, no. 4, pp. 662–667, November 2001. [5] A. R. Đorđević and D. V. Tošić, “Causality of circuit and electromagnetic-field models”, in Proc. of 5th European Conference on Circuits and Systems for Communications (ECCSC'10), Belgrade, Serbia, 2010, pp. 12–21. [6] P. J. W. Debye, Polar Molecules, New York: The Chemical Catalog Company, 1929. [7] C. F. Bohren, “What did Kramers and Kronig do and how did they do it?”, Eur. J. Phys., vol. 31, pp. 573–577, 2010. [8] F. M. Tesche, “On the use of the Hilbert transform for processing measured CW data”, IEEE Trans. Electromagn. Compat., vol. 34, no. 3, August 1992, pp. 259–266. [9] T. Meissner and F. J. Wentz, “The complex dielectric constant of pure and sea water from microwave satellite observations”, IEEE Trans. Geoscience Remote Sens., vol. 42, no. 9, pp. 1836–1849, September 2004. [10] G. D. Smith and B. J. Stanton, Soil Parameters from Fort A.P. Hill Soil Permittivity and Conductivity Measurements for the Wide Area Airborne Minefield Detection Program, Army Research Laboratory, Adelphi, MD, ARL-TR-3049, Sept. 2003. Available: http://www.arl.army.mil/arlreports/2003/ARL-TR- 3049.pdf. [11] Ansys. (2012, May). Automating the SI Design Flow for HFSS. [Online]. Available: http://www.ansys.com/ staticassets/ANSYS/Conference/Confidence/Minneapolis/Downloads/automating-si-design-flow-for-ansys- hfss-1.pdf [12] Agilent Technologies. (2009). About Dielectric Loss Models. [Online]. Available: http://edocs.soco.agilent. com/display/ads2009/About+Dielectric+Loss+Models [13] Simberian Inc. (2008, Sept.). Modeling Frequency-Dependent Dielectric Loss and Dispersion for Multi- Gigabit Data Channels. [Online]. Available: [14] http://www.simberian.com/AppNotes/ModelingDielectrics_2008_06.pdf [15] J. Barthel, K. Bachhuber, R. Buchner, H. Hetzenauer, and M. Kleebauer, “A computer-controlled system of transmission lines for the determination of the complex permittivity of lossy liquids between 8.5 and 90 GHz”, Ber. Bunsenges. Phys. Chem., vol. 95, no. 8, pp. 853–859, 1991. [16] WIPL-D 2-D Solver. Available: http://www.wipl-d.com/