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

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.


Introduction
A plasma dipole is an antenna with a radiating structure based on a plasma element instead of a metallic conductor [1][2][3][4][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][9][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][14][15][16][17][18].Some aspects of plasma antennas have already been treated with FDTD techniques [19][20][21][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 Figure 1.The method described herein is more efficient than equivalent eigenfunction expansions and can also provide field solutions with greater resolution and accuracy.
The remaining of this work is organized as follows.In Section 2, 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.

Theory
The electromagnetic wave propagation in a cold, collisional, isotropic plasma with no background magnetic field can be characterized by the relative permittivity with 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 Applying an inverse Fourier transform and some algebraic manipulations allows (3) to be written, in the time domain, as The relevant equations that govern the evolution of the electromagnetic fields within the plasma are then given, in the time domain, by

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 N x × N y spatial cells for a simulation run for N t 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  Region 2  3 International Journal of Antennas and Propagation 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 Electromagnetic fields are spatially discretized over a finite Cartesian grid as per Figure 2, 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
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.
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 where The update equation for the electric field requires knowledge of the plasma current at time step n + 1/2, which is obtained by a time average such that the update equation for the electric field within the plasma is given by where These update equations naturally take into account any spatial variations in both the electron density and the (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 space

Incident Wave Generation.
The incident wave under consideration for the FDTD technique is a plane wave propagating towards the right (positive x-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.

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.x which can be readily implemented concurrently to the leapfrogging algorithm at all grid points by means of 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 Section 3.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.
With the equivalent electric and magnetic currents at the fictitious boundary, the electromagnetic potentials for a twodimensional problem can be written in terms of Green's functions as 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 Additionally, by considering an observation point in the far field, the asymptotic expression for the Hankel function can be used, resulting in, for the z-component of the electric field,  where the angle ψ is given by cos ψ = ρˆ• ρˆ0.The numerical scattering amplitude is then given by 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.

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 N t = 1200, and number of spatial cells N = N x × N y = 63001.Figure 5 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.
Figures 6-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,  International Journal of Antennas and Propagation 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 12 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.
These 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 (Figures 6(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 finitedifference algorithms.

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  where n 0 is the plasma density at the center of the plasma, ρ is the distance from the center, and r is the cylinder radius.
Figure 13 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.Figures 14-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.
One 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 f in = 5 GHz in Figure 14 with the results for n 0 = 5 × 10 18 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   16 International Journal of Antennas and Propagation 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.These 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 density n.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 Figure 22.The maximum absolute point-wise error between the approximation and the density function is less than 2.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 timedependent 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 a z-component, which caused  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   19 International Journal of Antennas and Propagation dimensions, thus allowing for the full simulation of an entire device like a plasma antenna or even the interaction between multiple devices operating simultaneously.

Figure 1 :
Figure 1: Plane wave incidence on the plasma cylinder.

1 Figure 2 :
Figure2: 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).
Region separation for the TFSF technique.Region 1 consists of total fields while region 2 consists only of scattered fields Region 1 b) Detailed view of the field components adjacent to the TFSF boundary.A computational cell at the corner of the TFSF boundary is shown

Figure 3 :
Figure 3: Total-field scattered-field technique used to generate the incident fields and calculate the scattered far fields.

3. 4 .Figure 5 :
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 to f in = 10 GHz, electron density is set to n = 5 × 10 17 m −3 , and electron collision frequency is set to ν = 500 × 10 6 Hz.

Figure 6 : 3 40 3 40 8 (
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 to n 0 = 5 × 10 17 m −3 , and electron collision frequency is set to ν = 500 × 10 6 Hz.

Figure 7 :
Figure 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 to f in = 10 GHz, and electron collision frequency is set to ν = 500 × 10 6 Hz.

Figure 8 :
Figure 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 to f in = 10 GHz, and electron density is set to n 0 = 5 × 10 17 m −3 .

Figure 9 :
Figure9: 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 to n 0 = 5 × 10 17 m −3 , and electron collision frequency is set to ν = 500 × 10 6 Hz.

Figure 10 :
Figure 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 to f in = 10 GHz, and electron collision frequency is set to ν = 500 × 10 6 Hz.

Figure 11 :Figure 12 :
Figure 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 to f in = 10 GHz, and electron density is set to n 0 = 5 × 10 17 m −3 .

Figure 13 :
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 to f in = 10 GHz, central electron density is set to n 0 = 5 × 10 17 m −3 , and electron collision frequency is set to ν = 500 × 10 6 Hz.

Figure 14 : 3 40
Figure 14: Magnitude of the electric field, in volts/meter, within the inhomogeneous plasma cylinder obtained by the FDTD simulation for different incident frequencies f in .Central electron density is set to n 0 = 5 × 10 17 m −3 , and electron collision frequency is set to ν = 500 × 10 6 Hz.

Figure 15 :
Figure 15: Magnitude of the electric field, in volts/meter, within the inhomogeneous plasma cylinder obtained by the FDTD simulation for different central electron densities n 0 .Incident wave frequency is set to f in = 10 GHz, and electron collision frequency is set to ν = 500 × 10 6 Hz.

Figure 16 :−
Figure 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 f in = 10 GHz, and central electron density is set to n 0 = 5 × 10 17 m −3 .

Figure 17 :−
Figure17: 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 to n 0 = 5 × 10 17 m −3 , and electron collision frequency is set to ν = 500 × 10 6 Hz.

Figure 18 :
Figure18: 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 densities n 0 .Incident wave frequency is set to f in = 10 GHz, and electron collision frequency is set to ν = 500 × 10 6 Hz.

Figure 19 :
Figure 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 f in = 10 GHz, and central electron density is set to n 0 = 5 × 10 17 m −3 .
ν and f in

Figure 21 :
Figure 21: Behavior shift for the plasma dielectric permittivity as a function of the local plasma density n, the plasma collision frequency ν, and the incident wave frequency f in .Red region represents ℜ ε r > 0, and blue region represents ℜ ε r < 0.

Figure 22 :
Figure 22:  Comparison between the actual inhomogeneous density function for the plasmas considered herein and its staircase approximation in the FDTD grid.