# Finite-Difference Time Domain Techniques Applied to Electromagnetic Wave Interactions with Inhomogeneous Plasma Structures

Authors: Julio L. Nicolini; José Ricardo Bergmann
Journal: International Journal of Antennas and Propagation
 (2018)

---

## Abstract

Motivated by the emerging field of plasma antennas, electromagnetic wave propagation in and scattering by inhomogeneous plasma structures are studied through finite-difference time domain (FDTD) techniques. These techniques have been widely used in the past to study propagation near or through the ionosphere, and their extension to plasma devices such as antenna elements is a natural development. Simulation results in this work are validated with comparisons to solutions obtained by eigenfunction expansion techniques well supported by the literature and are shown to have an excellent agreement. The advantages of using FDTD simulations for this type of investigation are also outlined; in particular, FDTD simulations allow for field solutions to be developed at lower computational cost and greater resolution than equivalent eigenfunction methods for inhomogeneous plasmas and are applicable to arbitrary plasma properties such as spatially or time-varying inhomogeneities and collision frequencies, as well as allowing transient effects to be studied as the field solutions are obtained in the time domain.


---

## Body

## 1. Introduction

A plasma dipole is an antenna with a radiating structure based on a plasma element instead of a metallic conductor [1–5]. The plasma, kept activated by an ionizing source, is the conducting material that acts as the source for electromagnetic fields, which can then be modulated and used to carry information in telecommunication links. This is achieved by applying a secondary signal source to the plasma, which will then reemit the signal as electromagnetic radiation. Plasma antennas have several technological advantages, such as their radiation properties being electronically controlled, which results in an antenna that is not restricted to its original fabrication characteristics and can be reconfigured with simplicity and in real time [6, 7]. If the ionizing source is turned off, the plasma antenna is deactivated, becoming an inert element that is invisible to electromagnetic radiation, eliminating coupling problems with other antenna elements.These advantages can be employed in solutions to one of the great challenges in the synthesis of antennas for telecommunication links, which is the lack of flexibility in changing parameters, like, for example, the radiation pattern, once the antenna is deployed. This problem can be avoided by application of antenna arrays, where control over the excitation of each element allows the radiation pattern to be conformed. These antenna arrays are particularly interesting for use in reconfigurable antennas in mobile communications and satellite link applications [8–10], but the presence of metallic elements creates additional complexity in the synthesis process due to parasitic interactions. Plasma antennas provide an alternative to metallic elements in these arrays, as they remain inert when deactivated and therefore parasitic interactions are minimized.Despite the advantages in utilizing plasmas as conducting elements in antennas, there are still several difficulties and obstacles for this technology. Obtaining the real characteristics of a plasma antenna, in a generic situation, requires the complete description of the plasma configuration, which creates theoretical and numerical difficulties when trying to analyze such systems, as outlined in [7]. In particular, the behavior of electromagnetic waves within a plasma structure is not always trivial to be analyzed.The theoretical investigation of electromagnetic field behavior within a cylindrical inhomogeneous plasma structure is usually carried out through eigenfunction expansions [11, 12], which consists of expanding the electromagnetic field in Bessel functions, or other eigenfunctions appropriate to the problem’s geometry, and then finding the unknown expansion coefficients by application of boundary conditions within the plasma and at the plasma container’s boundaries. There are limitations to this method, however, such as computational costs in developing the field solutions and the comparatively low spatial resolution obtained for the fields internal to the plasma structure or the requirement that the inhomogeneities be of special cases so that an analytical solution can be found.In this paper, a simulation scheme based on FDTD techniques is proposed to investigate the problem of electromagnetic propagation within arbitrary inhomogeneous plasma structures appropriate for telecommunication applications and the waves scattered by such structures. Such simulations have already been successfully applied to the study of electromagnetic propagation through or near the ionosphere [13–18]. Some aspects of plasma antennas have already been treated with FDTD techniques [19–22], but simulated results for electromagnetic fields pertaining to waves propagating within the plasma structure or scattered by it have not been presented to the authors’ knowledge, and aforementioned works explore only homogeneous `plasma.For the ease of comparison with results from the literature, the plasma under consideration is inhomogeneous, cold, unmagnetized, and collisional. The plasma is also assumed to be confined to a cylindrical structure and under illumination from a transverse magnetic (TM) plane electromagnetic wave with electric field parallel to the cylinder’s axis and in a situation of normal incidence, as shown in Figure1. The method described herein is more efficient than equivalent eigenfunction expansions and can also provide field solutions with greater resolution and accuracy.Figure 1
Plane wave incidence on the plasma cylinder.The remaining of this work is organized as follows. In Section2, the theory utilized to describe the plasma structure in both the eigenfunction solution and the FDTD simulation is briefly outlined. In Section 3, the particulars of simulating the plasma behavior with the FDTD algorithm are presented. Section 4 provides validating results for homogeneous plasmas through comparisons with the usual eigenfunction methods as well as results that explore the characteristics of electromagnetic propagation through the inhomogeneous plasma. Section 5 provides the concluding remarks.


## 2. Theory

The electromagnetic wave propagation in a cold, collisional, isotropic plasma with no background magnetic field can be characterized by the relative permittivity(1)εr=1−wp2w2−iwv,with
(2)wp=ne2ε0m,where ωp is the plasma frequency, ν is the electron collision frequency, ω is the angular frequency of the propagating electromagnetic wave, n is the electron density, e is the electron’s charge, m is its mass, and 0 is the permittivity of free space.Equating (1) with the relative permittivity of an arbitrary lossy material allows an effective conductivity σef and associated plasma current density in the frequency domain to be found as
(3)σef=ε0wp2v+iw,J~→p=σefE~→.Applying an inverse Fourier transform and some algebraic manipulations allows (3) to be written, in the time domain, as
(4)∂∂tJ→p=ε0wp2E→−vJ→p.The relevant equations that govern the evolution of the electromagnetic fields within the plasma are then given, in the time domain, by(5)∂∂tH→=−μ0∇×E→,∂∂tE→=1ε0∇×H→−J→p.


## 3. FDTD Technique

The FDTD technique consists of discretizing Maxwell’s equations through the use of finite difference approximations for the derivatives. This section will provide the basic characteristics of the FDTD method utilized in this work.


### 3.1. Computational Domain and Update Equations

The basic scheme laid out by Yee in his seminal paper [23] is used to create the two-dimensional computational domain with Nx × Ny spatial cells for a simulation run for Nt steps. Assuming the time step is given by ∆t and the spatial steps are given by ∆s = ∆x = ∆y, a function f of space and time evaluated at a discrete point in space-time (n∆t, i∆x, and j∆y), where n, i, and j are integers or half-integers, is denoted as
(6)fnΔt,iΔx,jΔy=fni,j.Electromagnetic fields are spatially discretized over a finite Cartesian grid as per Figure2, and the cylinder’s boundary is realized with a staircase Cartesian approximation. Of note is that electric and magnetic fields are shifted by half steps from each other both spatially and in time, that is, electric fields are stored at integer times and positions, while magnetic fields are stored at half-integer times and positions. This allows for Maxwell’s equations to be readily discretized as
(7)Hxn+1/2i,j+12=Hxn−1/2i,j+12−μ0ΔtΔsEzni,j+1−Ezni,j,Hyn+1/2i+12,j=Hyn−1/2i+12,j−μ0ΔtΔsEzni+1,j−Ezni,j,Ezn+1i,j=Ezni,j+Δtε0ΔsHyn+1/2i+12,j−Hyn+1/2i−12,j−Hxn+1/2i,j+12−Hxn+1/2i,j−12,which can be used to propagate electromagnetic waves through a discretized free space within the computational domain. The equations for propagation within the plasma remain unchanged for the magnetic fields, but the plasma current that arises in (3) has to be addressed in the update equation for the electromagnetic field.Figure 2
Generic Cartesian spatial cell (i,j) and surrounding cells in the 2D computational domain. Stored values for each cell are the z-component of the electric field at time n and position (i,j), the x-component of the magnetic field at time n + 1/2 and position (i,j + 1/2), and the y-component of the magnetic field at time n + 1/2 and position (i + 1/2,j).This term can be efficiently handled by means of an auxiliary differential equation (ADE) formulation [24]. The updated equation for the current term is given by
(8)Jzn+1i,j=Kpi,jJzni,j+βpi,jEzn+1i,j+Ezni,j,where
(9)Kpi,j=1−vi,jΔt/21+vi,jΔt/2,βpi,j=wpi,j2ε0Δt/21+vi,jΔt/2.The update equation for the electric field requires knowledge of the plasma current at time stepn + 1/2, which is obtained by a time average such that the update equation for the electric field within the plasma is given by
(10)Ezn+1i,j=2ε0−Δβtpi,j2ε0+Δβtpi,jEzni,j+2Δt2ε0+Δβtpi,j∇×H→n+1/2−1+Kpi,j2Jzni,j,where
(11)∇×H→n+1/2=Hyn+1/2i+12,j−Hyn+1/2i−12,j−Hxn+1/2i,j+12−Hxn+1/2i,j−12.These update equations naturally take into account any spatial variations in both the electron density and the collision frequency, as the values for each are different at each computational cell.


### 3.2. Incident Wave Generation

The incident wave under consideration for the FDTD technique is a plane wave propagating towards the right (positivex-axis) of the computational domain. It is generated by means of a total-field scattered-field (TFSF) formulation [25] as shown in Figure 3. This region separation will also be exploited to calculate the scattered far field from the structure, discussed in Section 3.5. The computational domain is divided into two regions, and an auxiliary one-dimensional FDTD simulation is concurrently run to account for the propagation of the incident wave. The field values from the incident wave are then directly added and subtracted in the field update equations in computational cells surrounding the region separation.Figure 3
Total-field scattered-field technique used to generate the incident fields and calculate the scattered far fields.
(a)
Region separation for the TFSF technique. Region 1 consists of total fields while region 2 consists only of scattered fields
(b)
Detailed view of the field components adjacent to the TFSF boundary. A computational cell at the corner of the TFSF boundary is shown


### 3.3. Grid Termination

Absorbing boundary conditions (ABCs) based on one-way wave equations [26, 27] are used to analytically absorb incoming waves at the computational boundaries and simulate propagation towards infinity to prevent unphysical backscattered waves from polluting the numerical results. Mur’s ABC is only valid in a vacuum, but since the TFSF technique is used to generate, the incident wave is implemented in vacuum cells surrounding the plasma cylinder that condition is automatically fulfilled.


### 3.4. Steady State Fields

The FDTD technique provides a time domain solution for the electromagnetic fields, but the frequency domain solution can be obtained by an efficient implementation of in-step Fourier transforms on the fields [25]. Specifically, any field component f can be written, in the frequency domain, as
(12)f~w=∫−∞∞fte−iwtdt=∫0NtΔtfte−iwtdt,which can be readily implemented concurrently to the leap-frogging algorithm at all grid points by means of
(13)f~i,j,w≈13f0i,j+2∑nevenfni,je−iwnΔt+4∑noddfni,jeiwnΔt+fNti,je−iwNtΔt,where Simpson’s rule was employed to numerically evaluate the Fourier integral [28].


### 3.5. Near-to-Far-Field Transformation

With the steady-state fields calculated as described in Section3.4 from the nodes in region 2 of the TFSF technique, a near-to-far-field transformation can be employed to numerically obtain the scattering amplitude A(φ) of the structure of interest. The strategy consists of constructing a fictitious boundary within region 2 and taking advantage of the equivalence principle as shown in Figure 4.Figure 4
Equivalence principle.
(a)
Fictitious boundary for the problem at hand. Fields inside and outside the boundary are continuous across it
(b)
Equivalent problem with fields inside the fictitious boundary set to zero and material properties set to that of free spaceWith the equivalent electric and magnetic currents at the fictitious boundary, the electromagnetic potentials for a two-dimensional problem can be written in terms of Green’s functions as(14)A→ρ→=−iμ4∮LJ→sρ→′H02kρ→−ρ→′dl′,F→ρ→=−iε4∮LM→sρ→′H02kρ→−ρ→′dl′,where primed coordinates denote points upon the fictitious boundary and the closed path integral is calculated over the boundary. With the electromagnetic potentials, the scattered fields can be found by means of
(15)E→=−1ε∇×F→+iwA→+iwμε∇∇⋅A→,H→=1μ∇×A→+iwF→+iwμε∇∇⋅F→.Additionally, by considering an observation point in the far field, the asymptotic expression for the Hankel function can be used, resulting in, for thez-component of the electric field,
(16)Ezsρ→=−i8πke−ikρρ∮Lwμ0J→sρ→′⋅ϕ^eikρ′cosψdl′,where the angle ψ is given by cosψ=ρˆ·ρˆ0. The numerical scattering amplitude is then given by
(17)Aϕ=−i8πk∫Lwμ0J→sρ→′⋅z^−kM→sρ→′⋅ϕ^eikp′cosψdl′,which can be readily calculated with the application of Simpson’s rule, as all necessary quantities are known at the end of the FDTD run.


## 3.1. Computational Domain and Update Equations

The basic scheme laid out by Yee in his seminal paper [23] is used to create the two-dimensional computational domain with Nx × Ny spatial cells for a simulation run for Nt steps. Assuming the time step is given by ∆t and the spatial steps are given by ∆s = ∆x = ∆y, a function f of space and time evaluated at a discrete point in space-time (n∆t, i∆x, and j∆y), where n, i, and j are integers or half-integers, is denoted as
(6)fnΔt,iΔx,jΔy=fni,j.Electromagnetic fields are spatially discretized over a finite Cartesian grid as per Figure2, and the cylinder’s boundary is realized with a staircase Cartesian approximation. Of note is that electric and magnetic fields are shifted by half steps from each other both spatially and in time, that is, electric fields are stored at integer times and positions, while magnetic fields are stored at half-integer times and positions. This allows for Maxwell’s equations to be readily discretized as
(7)Hxn+1/2i,j+12=Hxn−1/2i,j+12−μ0ΔtΔsEzni,j+1−Ezni,j,Hyn+1/2i+12,j=Hyn−1/2i+12,j−μ0ΔtΔsEzni+1,j−Ezni,j,Ezn+1i,j=Ezni,j+Δtε0ΔsHyn+1/2i+12,j−Hyn+1/2i−12,j−Hxn+1/2i,j+12−Hxn+1/2i,j−12,which can be used to propagate electromagnetic waves through a discretized free space within the computational domain. The equations for propagation within the plasma remain unchanged for the magnetic fields, but the plasma current that arises in (3) has to be addressed in the update equation for the electromagnetic field.Figure 2
Generic Cartesian spatial cell (i,j) and surrounding cells in the 2D computational domain. Stored values for each cell are the z-component of the electric field at time n and position (i,j), the x-component of the magnetic field at time n + 1/2 and position (i,j + 1/2), and the y-component of the magnetic field at time n + 1/2 and position (i + 1/2,j).This term can be efficiently handled by means of an auxiliary differential equation (ADE) formulation [24]. The updated equation for the current term is given by
(8)Jzn+1i,j=Kpi,jJzni,j+βpi,jEzn+1i,j+Ezni,j,where
(9)Kpi,j=1−vi,jΔt/21+vi,jΔt/2,βpi,j=wpi,j2ε0Δt/21+vi,jΔt/2.The update equation for the electric field requires knowledge of the plasma current at time stepn + 1/2, which is obtained by a time average such that the update equation for the electric field within the plasma is given by
(10)Ezn+1i,j=2ε0−Δβtpi,j2ε0+Δβtpi,jEzni,j+2Δt2ε0+Δβtpi,j∇×H→n+1/2−1+Kpi,j2Jzni,j,where
(11)∇×H→n+1/2=Hyn+1/2i+12,j−Hyn+1/2i−12,j−Hxn+1/2i,j+12−Hxn+1/2i,j−12.These update equations naturally take into account any spatial variations in both the electron density and the collision frequency, as the values for each are different at each computational cell.


## 3.2. Incident Wave Generation

The incident wave under consideration for the FDTD technique is a plane wave propagating towards the right (positivex-axis) of the computational domain. It is generated by means of a total-field scattered-field (TFSF) formulation [25] as shown in Figure 3. This region separation will also be exploited to calculate the scattered far field from the structure, discussed in Section 3.5. The computational domain is divided into two regions, and an auxiliary one-dimensional FDTD simulation is concurrently run to account for the propagation of the incident wave. The field values from the incident wave are then directly added and subtracted in the field update equations in computational cells surrounding the region separation.Figure 3
Total-field scattered-field technique used to generate the incident fields and calculate the scattered far fields.
(a)
Region separation for the TFSF technique. Region 1 consists of total fields while region 2 consists only of scattered fields
(b)
Detailed view of the field components adjacent to the TFSF boundary. A computational cell at the corner of the TFSF boundary is shown


## 3.3. Grid Termination

Absorbing boundary conditions (ABCs) based on one-way wave equations [26, 27] are used to analytically absorb incoming waves at the computational boundaries and simulate propagation towards infinity to prevent unphysical backscattered waves from polluting the numerical results. Mur’s ABC is only valid in a vacuum, but since the TFSF technique is used to generate, the incident wave is implemented in vacuum cells surrounding the plasma cylinder that condition is automatically fulfilled.


## 3.4. Steady State Fields

The FDTD technique provides a time domain solution for the electromagnetic fields, but the frequency domain solution can be obtained by an efficient implementation of in-step Fourier transforms on the fields [25]. Specifically, any field component f can be written, in the frequency domain, as
(12)f~w=∫−∞∞fte−iwtdt=∫0NtΔtfte−iwtdt,which can be readily implemented concurrently to the leap-frogging algorithm at all grid points by means of
(13)f~i,j,w≈13f0i,j+2∑nevenfni,je−iwnΔt+4∑noddfni,jeiwnΔt+fNti,je−iwNtΔt,where Simpson’s rule was employed to numerically evaluate the Fourier integral [28].


## 3.5. Near-to-Far-Field Transformation

With the steady-state fields calculated as described in Section3.4 from the nodes in region 2 of the TFSF technique, a near-to-far-field transformation can be employed to numerically obtain the scattering amplitude A(φ) of the structure of interest. The strategy consists of constructing a fictitious boundary within region 2 and taking advantage of the equivalence principle as shown in Figure 4.Figure 4
Equivalence principle.
(a)
Fictitious boundary for the problem at hand. Fields inside and outside the boundary are continuous across it
(b)
Equivalent problem with fields inside the fictitious boundary set to zero and material properties set to that of free spaceWith the equivalent electric and magnetic currents at the fictitious boundary, the electromagnetic potentials for a two-dimensional problem can be written in terms of Green’s functions as(14)A→ρ→=−iμ4∮LJ→sρ→′H02kρ→−ρ→′dl′,F→ρ→=−iε4∮LM→sρ→′H02kρ→−ρ→′dl′,where primed coordinates denote points upon the fictitious boundary and the closed path integral is calculated over the boundary. With the electromagnetic potentials, the scattered fields can be found by means of
(15)E→=−1ε∇×F→+iwA→+iwμε∇∇⋅A→,H→=1μ∇×A→+iwF→+iwμε∇∇⋅F→.Additionally, by considering an observation point in the far field, the asymptotic expression for the Hankel function can be used, resulting in, for thez-component of the electric field,
(16)Ezsρ→=−i8πke−ikρρ∮Lwμ0J→sρ→′⋅ϕ^eikρ′cosψdl′,where the angle ψ is given by cosψ=ρˆ·ρˆ0. The numerical scattering amplitude is then given by
(17)Aϕ=−i8πk∫Lwμ0J→sρ→′⋅z^−kM→sρ→′⋅ϕ^eikp′cosψdl′,which can be readily calculated with the application of Simpson’s rule, as all necessary quantities are known at the end of the FDTD run.


## 4. Results

### 4.1. Validation Results with Homogeneous Plasma

To validate the FDTD technique, it is applied to a homogeneous plasma and the solution compared with solutions obtained from eigenfunction expansions, which can be efficiently calculated for homogeneous cases. The simulation parameters are as follows: time discretization ∆t = 2.35702 × 10−12 seconds, spatial discretization ∆s = 1 × 10−3 meters, maximum temporal step Nt = 1200, and number of spatial cells N = Nx × Ny = 63001.Figure5 shows the time evolution of the intensity of the electric field across the computational domain for a homogeneous plasma, showing the transient effects of the plasma on the electromagnetic wave such as refraction at the cylinder boundaries as well as internal reflections within the plasma.Figure 5
Time evolution of the intensity of the electric field, in volts/meter, across the computational domain, for the homogeneous case where incident wave frequency is set tofin = 10 GHz, electron density is set to n = 5 × 1017 m−3, and electron collision frequency is set to ν = 500 × 106 Hz.
(a)
n = 0
(b)
n = 200
(c)
n = 400
(d)
n = 600Figures6–8 show comparisons between the magnitudes of the electric field within the homogeneous plasma cylinder obtained by the FDTD simulation and by the eigenfunction method for different incident frequencies, electron densities, and collision frequencies, respectively. Figures 9–10 show the same, but for linear cuts along the x- and y-axes, which allows for more precise comparisons between the two methods than the two-dimensional colour plots.Figure 6
Comparison between the magnitude of the electric field, in volts/meter, within the cylinder obtained by the eigenfunction method and the FDTD simulation for different incident frequencies; electron density is set ton0 = 5 × 1017 m−3, and electron collision frequency is set to ν = 500 × 106 Hz.
(a)
Eigenfunction,fin = 5 GHz
(b)
FDTD,fin = 5 GHz
(c)
Eigenfunction,fin = 10 GHz
(d)
FDTD,fin = 10 GHz
(e)
Eigenfunction,fin = 15 GHz
(f)
FDTD,fin = 15 GHzFigure 7
Comparison between the magnitude of the electric field, in volts/meter, within the cylinder obtained by the eigenfunction method and the FDTD simulation for different electron densities; incident wave frequency is set tofin = 10 GHz, and electron collision frequency is set to ν = 500 × 106 Hz.
(a)
Eigenfunction,n0 = 5 × 1017 m − 3
(b)
FDTD,n0 = 5 × 1017 m−3
(c)
Eigenfunction,n0 = 1 × 1018 m−3
(d)
FDTD,n0 = 1 × 1018 m−3
(e)
Eigenfunction,n0 = 5 × 1018 m−3
(f)
FDTD,n0 = 5 × 1018 m−3Figure 8
Comparison between the magnitude of the electric field, in volts/meter, within the cylinder obtained by the eigenfunction method and the FDTD simulation for different collision frequencies; incident wave frequency is set tofin = 10 GHz, and electron density is set to n0 = 5 × 1017 m−3.
(a)
Eigenfunction,ν = 500 × 106 Hz
(b)
FDTD,ν = 500 × 106 Hz
(c)
Eigenfunction,ν = 500 × 107 Hz
(d)
FDTD,ν = 500 × 107 Hz
(e)
Eigenfunction,ν = 500 × 108 Hz
(f)
FDTD,ν = 500 × 108 HzFigure 9
Comparison between the magnitude of the electric field, in volts/meter, inside the plasma for linear cuts through the cylinder obtained by the eigenfunction method and the FDTD simulation for different incident frequencies; electron density is set ton0 = 5 × 1017 m−3, and electron collision frequency is set to ν = 500 × 106 Hz.
(a)
Variablesx and y = 0
(b)
Variablesy and x = 0Figure 10
Comparison between the magnitude of the electric field, in volts/meter, inside the plasma for linear cuts through the cylinder obtained by the eigenfunction method and the FDTD simulation for different electron densities; incident wave frequency is set tofin = 10 GHz, and electron collision frequency is set to ν = 500 × 106 Hz.
(a)
Variablesx and y = 0
(b)
Variablesy and x = 0Figure 11
Comparison between the magnitude of the electric field, in volts/meter, inside the plasma for linear cuts through the cylinder obtained by the eigenfunction method and the FDTD simulation for different collision frequencies; incident wave frequency is set tofin = 10 GHz, and electron density is set to n0 = 5 × 1017 m−3.
(a)
Variablesx and y = 0
(b)
Variablesy and x = 0Figure12 shows a comparison between the resulting scattering amplitude for the electric field obtained by the eigenfunction expansion method and the FDTD simulation for varying values of the homogeneous plasma parameters. Excellent agreement is found between the two approaches.Figure 12
Comparison between the scattering amplitudes obtained by the eigenfunction method and the FDTD simulation for different variable parameters of the homogeneous plasma.
(a)
Variable incident frequencyfin
(b)
Variable collision frequencyν
(c)
Variable electron densitynThese results show that the FDTD simulations provide solutions that differ less than 1% from those obtained by the eigenfunction method that is widely used in the literature in most cases and less than 5% in the worst case. It can also be seen that the FDTD technique allows the study of the time domain evolution of the wave propagation, so transient effects can be analyzed. Additionally, it can be seen that the algorithm is still precise even for the plasma that is overly dense in relation to the incident frequency (Figures6(b) and 7(f)) or that displays large losses (Figure 8(f)), so that spatial accuracy is only dependent on the mesh being fine enough to resolve the wave propagation, which is a well-known condition in finite-difference algorithms.


### 4.2. Investigation of Inhomogeneous Plasma

An inhomogeneous plasma is now investigated, with the same simulation parameters as the homogeneous case. The plasma inhomogeneity is defined by the quadratic density profile(18)nρ=n01−ρr2,where n0 is the plasma density at the center of the plasma, ρ is the distance from the center, and r is the cylinder radius.Figure13 shows the time evolution of the intensity of the electric field across the computational domain for an inhomogeneous plasma, showing the transient effects of the plasma on the electromagnetic wave such as the deflection in the propagation direction of the wave.Figure 13
Time evolution of the intensity of the electric field, in volts/meter, across the computational domain, for the inhomogeneous case where incident wave frequency is set tofin = 10 GHz, central electron density is set to n0 = 5 × 1017 m−3, and electron collision frequency is set to ν = 500 × 106 Hz.
(a)
n = 0
(b)
n = 200
(c)
n = 400
(d)
n = 600Figures14–16 show the magnitude of the electric field within the inhomogeneous plasma cylinder obtained by the FDTD simulation for different incident frequencies, central plasma densities, and electron collision frequencies, respectively. Figures 17–19 show the same, but for linear cuts along the x- and y-axes. Figure 20 shows the resulting scattering amplitude for the electric field obtained by the FDTD simulation for varying values of the inhomogeneous plasma parameters.Figure 14
Magnitude of the electric field, in volts/meter, within the inhomogeneous plasma cylinder obtained by the FDTD simulation for different incident frequenciesfin. Central electron density is set to n0 = 5 × 1017 m−3, and electron collision frequency is set to ν = 500 × 106 Hz.
(a)
fin = 5 GHz
(b)
fin = 10 GHz
(c)
fin = 15 GHzFigure 15
Magnitude of the electric field, in volts/meter, within the inhomogeneous plasma cylinder obtained by the FDTD simulation for different central electron densitiesn0. Incident wave frequency is set to fin = 10 GHz, and electron collision frequency is set to ν = 500 × 106 Hz.
(a)
n0 = 1 × 1017 m−3
(b)
n0 = 5 × 1017 m−3
(c)
n0 = 1 × 1018 m−3
(d)
n0 = 1.5 × 1018 m−3
(e)
n0 = 2 × 1018 m−3
(f)
n0 = 5 × 1018 m−3Figure 16
Magnitude of the electric field, in volts/meter, within the inhomogeneous plasma cylinder obtained by the FDTD simulation for different collision frequenciesν. Incident wave frequency is set to fin = 10 GHz, and central electron density is set to n0 = 5 × 1017 m−3.
(a)
ν = 500 × 106 Hz
(b)
ν = 500 × 107 Hz
(c)
ν = 500 × 108 HzFigure 17
Magnitude of the electric field, in volts/meter, inside the plasma for linear cuts through the cylinder obtained by the FDTD simulation for different incident frequencies. Central electron density is set ton0 = 5 × 1017 m−3, and electron collision frequency is set to ν = 500 × 106 Hz.
(a)
Variablesx and y = 0
(b)
Variablesy and x = 0Figure 18
Magnitude of the electric field, in volts/meter, inside the plasma for linear cuts through the cylinder obtained by the FDTD simulation for different central electron densitiesn0. Incident wave frequency is set to fin = 10 GHz, and electron collision frequency is set to ν = 500 × 106 Hz.
(a)
Variablesx and y = 0
(b)
Variablesy and x = 0Figure 19
Magnitude of the electric field, in volts/meter, inside the plasma for linear cuts through the cylinder obtained by the FDTD simulation for different collision frequenciesν. Incident wave frequency is set to fin = 10 GHz, and central electron density is set to n0 = 5 × 1017 m−3.
(a)
Variablesx and y = 0
(b)
Variablesy and x = 0Figure 20
Comparison between the scattering amplitudes obtained by the eigenfunction method and the FDTD simulation for different variable parameters of the inhomogeneous plasma.
(a)
Variable incident frequencyfin
(b)
Variable collision frequencyν
(c)
Variable central electron densityn0One qualitative analysis is straightforward: as expected from inspecting (1), variations in the plasma frequency ωp (which depends on the plasma density) change the behavior of electromagnetic propagation inversely to variations in both wave frequency ω and electron collision frequency ν. This effect can be seen, for example, by comparing the results for fin = 5 GHz in Figure 14 with the results for n0 = 5 × 1018 m−3 in Figure 15.Additionally, another qualitative analysis of the results shows a phenomenon of wave path deflection when the propagation is through inhomogeneous plasmas. This is due to the spatial variation in the electron density, which in turn causes a spatial variation in the refraction index of the plasma medium. Continuous spatial variations in refraction indexes, in turn, are well-known to cause ray deflection.In broad terms, two different behaviors can be observed from the presented results: (1) electromagnetic waves penetrating the plasma and propagating while being conditioned by the plasma, that is, suffering dispersion, attenuation, and deflection, when appropriate to each situation’s characteristics, and (2) electromagnetic waves being reflected from the plasma and exhibiting very low penetration (or, for the inhomogeneous cases, very low penetration after a certain point in the inhomogeneous cylinder). These two different types of behavior are related to the real part of the plasma’s dielectric permittivity, with penetration possible forℜεr>0 and reflection occurring for ℜεr<0. This behavior shift is shown in Figure 21.Figure 21
Behavior shift for the plasma dielectric permittivity as a function of the local plasma densityn, the plasma collision frequency ν, and the incident wave frequency fin. Red region represents ℜεr>0, and blue region represents ℜεr<0.
(a)
n and ν
(b)
n and fin
(c)
ν and finThese results are consistent with previous qualitative results in the literature obtained by the application of the eigenfunction expansion technique, but solutions are obtained at lower computational cost and with greater resolution when applying the FDTD technique.A limitation that comes with the FDTD technique for inhomogeneous plasmas, however, is the staircasing effect on the local plasma densityn. Even with an overall spatially dependent profile, the density is taken to be constant within each grid cell, so there is an additional constraint that the spatial step must be small enough to properly approximate the desired profile up to a tolerance threshold that will depend on the application and desired error measure for the solution.For the cases presented herein, the difference between the actual plasma density function and its staircase approximation in the FDTD grid is shown in Figure22. The maximum absolute point-wise error between the approximation and the density function is less than 2.5%.Figure 22
Comparison between the actual inhomogeneous density function for the plasmas considered herein and its staircase approximation in the FDTD grid.


## 4.1. Validation Results with Homogeneous Plasma

To validate the FDTD technique, it is applied to a homogeneous plasma and the solution compared with solutions obtained from eigenfunction expansions, which can be efficiently calculated for homogeneous cases. The simulation parameters are as follows: time discretization ∆t = 2.35702 × 10−12 seconds, spatial discretization ∆s = 1 × 10−3 meters, maximum temporal step Nt = 1200, and number of spatial cells N = Nx × Ny = 63001.Figure5 shows the time evolution of the intensity of the electric field across the computational domain for a homogeneous plasma, showing the transient effects of the plasma on the electromagnetic wave such as refraction at the cylinder boundaries as well as internal reflections within the plasma.Figure 5
Time evolution of the intensity of the electric field, in volts/meter, across the computational domain, for the homogeneous case where incident wave frequency is set tofin = 10 GHz, electron density is set to n = 5 × 1017 m−3, and electron collision frequency is set to ν = 500 × 106 Hz.
(a)
n = 0
(b)
n = 200
(c)
n = 400
(d)
n = 600Figures6–8 show comparisons between the magnitudes of the electric field within the homogeneous plasma cylinder obtained by the FDTD simulation and by the eigenfunction method for different incident frequencies, electron densities, and collision frequencies, respectively. Figures 9–10 show the same, but for linear cuts along the x- and y-axes, which allows for more precise comparisons between the two methods than the two-dimensional colour plots.Figure 6
Comparison between the magnitude of the electric field, in volts/meter, within the cylinder obtained by the eigenfunction method and the FDTD simulation for different incident frequencies; electron density is set ton0 = 5 × 1017 m−3, and electron collision frequency is set to ν = 500 × 106 Hz.
(a)
Eigenfunction,fin = 5 GHz
(b)
FDTD,fin = 5 GHz
(c)
Eigenfunction,fin = 10 GHz
(d)
FDTD,fin = 10 GHz
(e)
Eigenfunction,fin = 15 GHz
(f)
FDTD,fin = 15 GHzFigure 7
Comparison between the magnitude of the electric field, in volts/meter, within the cylinder obtained by the eigenfunction method and the FDTD simulation for different electron densities; incident wave frequency is set tofin = 10 GHz, and electron collision frequency is set to ν = 500 × 106 Hz.
(a)
Eigenfunction,n0 = 5 × 1017 m − 3
(b)
FDTD,n0 = 5 × 1017 m−3
(c)
Eigenfunction,n0 = 1 × 1018 m−3
(d)
FDTD,n0 = 1 × 1018 m−3
(e)
Eigenfunction,n0 = 5 × 1018 m−3
(f)
FDTD,n0 = 5 × 1018 m−3Figure 8
Comparison between the magnitude of the electric field, in volts/meter, within the cylinder obtained by the eigenfunction method and the FDTD simulation for different collision frequencies; incident wave frequency is set tofin = 10 GHz, and electron density is set to n0 = 5 × 1017 m−3.
(a)
Eigenfunction,ν = 500 × 106 Hz
(b)
FDTD,ν = 500 × 106 Hz
(c)
Eigenfunction,ν = 500 × 107 Hz
(d)
FDTD,ν = 500 × 107 Hz
(e)
Eigenfunction,ν = 500 × 108 Hz
(f)
FDTD,ν = 500 × 108 HzFigure 9
Comparison between the magnitude of the electric field, in volts/meter, inside the plasma for linear cuts through the cylinder obtained by the eigenfunction method and the FDTD simulation for different incident frequencies; electron density is set ton0 = 5 × 1017 m−3, and electron collision frequency is set to ν = 500 × 106 Hz.
(a)
Variablesx and y = 0
(b)
Variablesy and x = 0Figure 10
Comparison between the magnitude of the electric field, in volts/meter, inside the plasma for linear cuts through the cylinder obtained by the eigenfunction method and the FDTD simulation for different electron densities; incident wave frequency is set tofin = 10 GHz, and electron collision frequency is set to ν = 500 × 106 Hz.
(a)
Variablesx and y = 0
(b)
Variablesy and x = 0Figure 11
Comparison between the magnitude of the electric field, in volts/meter, inside the plasma for linear cuts through the cylinder obtained by the eigenfunction method and the FDTD simulation for different collision frequencies; incident wave frequency is set tofin = 10 GHz, and electron density is set to n0 = 5 × 1017 m−3.
(a)
Variablesx and y = 0
(b)
Variablesy and x = 0Figure12 shows a comparison between the resulting scattering amplitude for the electric field obtained by the eigenfunction expansion method and the FDTD simulation for varying values of the homogeneous plasma parameters. Excellent agreement is found between the two approaches.Figure 12
Comparison between the scattering amplitudes obtained by the eigenfunction method and the FDTD simulation for different variable parameters of the homogeneous plasma.
(a)
Variable incident frequencyfin
(b)
Variable collision frequencyν
(c)
Variable electron densitynThese results show that the FDTD simulations provide solutions that differ less than 1% from those obtained by the eigenfunction method that is widely used in the literature in most cases and less than 5% in the worst case. It can also be seen that the FDTD technique allows the study of the time domain evolution of the wave propagation, so transient effects can be analyzed. Additionally, it can be seen that the algorithm is still precise even for the plasma that is overly dense in relation to the incident frequency (Figures6(b) and 7(f)) or that displays large losses (Figure 8(f)), so that spatial accuracy is only dependent on the mesh being fine enough to resolve the wave propagation, which is a well-known condition in finite-difference algorithms.


## 4.2. Investigation of Inhomogeneous Plasma

An inhomogeneous plasma is now investigated, with the same simulation parameters as the homogeneous case. The plasma inhomogeneity is defined by the quadratic density profile(18)nρ=n01−ρr2,where n0 is the plasma density at the center of the plasma, ρ is the distance from the center, and r is the cylinder radius.Figure13 shows the time evolution of the intensity of the electric field across the computational domain for an inhomogeneous plasma, showing the transient effects of the plasma on the electromagnetic wave such as the deflection in the propagation direction of the wave.Figure 13
Time evolution of the intensity of the electric field, in volts/meter, across the computational domain, for the inhomogeneous case where incident wave frequency is set tofin = 10 GHz, central electron density is set to n0 = 5 × 1017 m−3, and electron collision frequency is set to ν = 500 × 106 Hz.
(a)
n = 0
(b)
n = 200
(c)
n = 400
(d)
n = 600Figures14–16 show the magnitude of the electric field within the inhomogeneous plasma cylinder obtained by the FDTD simulation for different incident frequencies, central plasma densities, and electron collision frequencies, respectively. Figures 17–19 show the same, but for linear cuts along the x- and y-axes. Figure 20 shows the resulting scattering amplitude for the electric field obtained by the FDTD simulation for varying values of the inhomogeneous plasma parameters.Figure 14
Magnitude of the electric field, in volts/meter, within the inhomogeneous plasma cylinder obtained by the FDTD simulation for different incident frequenciesfin. Central electron density is set to n0 = 5 × 1017 m−3, and electron collision frequency is set to ν = 500 × 106 Hz.
(a)
fin = 5 GHz
(b)
fin = 10 GHz
(c)
fin = 15 GHzFigure 15
Magnitude of the electric field, in volts/meter, within the inhomogeneous plasma cylinder obtained by the FDTD simulation for different central electron densitiesn0. Incident wave frequency is set to fin = 10 GHz, and electron collision frequency is set to ν = 500 × 106 Hz.
(a)
n0 = 1 × 1017 m−3
(b)
n0 = 5 × 1017 m−3
(c)
n0 = 1 × 1018 m−3
(d)
n0 = 1.5 × 1018 m−3
(e)
n0 = 2 × 1018 m−3
(f)
n0 = 5 × 1018 m−3Figure 16
Magnitude of the electric field, in volts/meter, within the inhomogeneous plasma cylinder obtained by the FDTD simulation for different collision frequenciesν. Incident wave frequency is set to fin = 10 GHz, and central electron density is set to n0 = 5 × 1017 m−3.
(a)
ν = 500 × 106 Hz
(b)
ν = 500 × 107 Hz
(c)
ν = 500 × 108 HzFigure 17
Magnitude of the electric field, in volts/meter, inside the plasma for linear cuts through the cylinder obtained by the FDTD simulation for different incident frequencies. Central electron density is set ton0 = 5 × 1017 m−3, and electron collision frequency is set to ν = 500 × 106 Hz.
(a)
Variablesx and y = 0
(b)
Variablesy and x = 0Figure 18
Magnitude of the electric field, in volts/meter, inside the plasma for linear cuts through the cylinder obtained by the FDTD simulation for different central electron densitiesn0. Incident wave frequency is set to fin = 10 GHz, and electron collision frequency is set to ν = 500 × 106 Hz.
(a)
Variablesx and y = 0
(b)
Variablesy and x = 0Figure 19
Magnitude of the electric field, in volts/meter, inside the plasma for linear cuts through the cylinder obtained by the FDTD simulation for different collision frequenciesν. Incident wave frequency is set to fin = 10 GHz, and central electron density is set to n0 = 5 × 1017 m−3.
(a)
Variablesx and y = 0
(b)
Variablesy and x = 0Figure 20
Comparison between the scattering amplitudes obtained by the eigenfunction method and the FDTD simulation for different variable parameters of the inhomogeneous plasma.
(a)
Variable incident frequencyfin
(b)
Variable collision frequencyν
(c)
Variable central electron densityn0One qualitative analysis is straightforward: as expected from inspecting (1), variations in the plasma frequency ωp (which depends on the plasma density) change the behavior of electromagnetic propagation inversely to variations in both wave frequency ω and electron collision frequency ν. This effect can be seen, for example, by comparing the results for fin = 5 GHz in Figure 14 with the results for n0 = 5 × 1018 m−3 in Figure 15.Additionally, another qualitative analysis of the results shows a phenomenon of wave path deflection when the propagation is through inhomogeneous plasmas. This is due to the spatial variation in the electron density, which in turn causes a spatial variation in the refraction index of the plasma medium. Continuous spatial variations in refraction indexes, in turn, are well-known to cause ray deflection.In broad terms, two different behaviors can be observed from the presented results: (1) electromagnetic waves penetrating the plasma and propagating while being conditioned by the plasma, that is, suffering dispersion, attenuation, and deflection, when appropriate to each situation’s characteristics, and (2) electromagnetic waves being reflected from the plasma and exhibiting very low penetration (or, for the inhomogeneous cases, very low penetration after a certain point in the inhomogeneous cylinder). These two different types of behavior are related to the real part of the plasma’s dielectric permittivity, with penetration possible forℜεr>0 and reflection occurring for ℜεr<0. This behavior shift is shown in Figure 21.Figure 21
Behavior shift for the plasma dielectric permittivity as a function of the local plasma densityn, the plasma collision frequency ν, and the incident wave frequency fin. Red region represents ℜεr>0, and blue region represents ℜεr<0.
(a)
n and ν
(b)
n and fin
(c)
ν and finThese results are consistent with previous qualitative results in the literature obtained by the application of the eigenfunction expansion technique, but solutions are obtained at lower computational cost and with greater resolution when applying the FDTD technique.A limitation that comes with the FDTD technique for inhomogeneous plasmas, however, is the staircasing effect on the local plasma densityn. Even with an overall spatially dependent profile, the density is taken to be constant within each grid cell, so there is an additional constraint that the spatial step must be small enough to properly approximate the desired profile up to a tolerance threshold that will depend on the application and desired error measure for the solution.For the cases presented herein, the difference between the actual plasma density function and its staircase approximation in the FDTD grid is shown in Figure22. The maximum absolute point-wise error between the approximation and the density function is less than 2.5%.Figure 22
Comparison between the actual inhomogeneous density function for the plasmas considered herein and its staircase approximation in the FDTD grid.


## 5. Conclusions

FDTD simulations are of great value in identifying different behaviors within the plasma structure and exploring the effects of the plasma parameters on the electromagnetic propagation as well as studying transient effects. In particular, for inhomogeneous cases or any case with other kinds of spatial complexity (e.g., complicated geometries for the structure), the computational cost of running analytical methods based on eigenfunction expansions is prohibitively high, and the results obtained lack good enough resolution to be precisely analyzed, two limitations that do not exist for FDTD simulations. The FDTD technique described herein allows for arbitrary time or spatially varying parameters to be incorporated in the simulation, as well as providing step-by-step transient solutions, two features that eigenfunction expansions lack.Understanding the effects of plasma inhomogeneities in electromagnetic wave propagation and scattering and being able to correctly simulate those effects are important steps in the design of plasma devices such as telecommunication antennas, especially when the inhomogeneities are time-dependent or present sharp spatial variations. The algorithm presented herein is an efficient solution that, to the authors’ knowledge, has not been previously applied to 2D plasma systems in the context of telecommunication devices.Future perspectives for this work include extending the numerical algorithm for a TEz-polarized incident wave. Due to the nature of TMz-polarized waves, the electric field was restricted to having only az-component, which caused the electron flux in the plasma to be solely in the nonsimulated axis and therefore ignored. Extending the algorithm to TEz-polarized waves creates the additional complexity of modelling the electron density variation caused by electron fluxes within the computational domain, as well as enforcing boundary conditions for the plasma currents excited within the structure.Another perspective is including ionization processes in the algorithm, which would allow the simulation of the start up and turn off of a device; so far, the plasma has been considered to be in a steady state of ionization; that is, the source responsible for ionization is considered to be active for a long time and recombination processes are ignored.With these extensions, the algorithm would be able to simulate fully self-consistent plasma systems in three spatial dimensions, thus allowing for the full simulation of an entire device like a plasma antenna or even the interaction between multiple devices operating simultaneously.



---

*Source: 3476462-2018-03-29.xml*