1. Introduction

The precursor work proposed by Pendry has established that it is possible, by using a coordinate transformation, to exclude all fields from a certain region without perturbing the vicinity [1]. The process consists in surrounding the object we want to make invisible by an anisotropic material whose permittivity and permeability are directly obtained from the so-called coordinate transformation.

Since then, invisible cloaking has been the matter of many studies spanning different approaches; by analytically solving Maxwell equations, it has been possible to theoretically analyze the cloaking phenomenon and interpret the underlying physics [2,3]. Furthermore, the first practical realization of a cloak has confirmed the model to be experimentally feasible [4]. Naturally, numerical simulations have an important role to play here through the useful information they are able to provide. In this sense, full-wave finite-element simulations have been successfully employed by Cummers et al. to confirm the effectiveness of the transformation [5]; in addition, they used simplified permittivity and permeability to show that the arising reduced cloaking materials exhibit good performance. Nevertheless, these calculations are performed with commercial software, such as Comsol Multiphysics, which does not allow the versatility needed in some situations. In this paper, the simulation of a cylindrical cloak is carried out with the Transmission Line Modeling (TLM) method [6]. To the best knowledge of the authors, it is the first time that a time domain method is employed for a full-wave simulation of a cloaking structure.

TLM is a low-frequency numerical method which has been extensively used for the modeling of wave propagation problems, mainly of electromagnetic nature, but also for problems in acoustics or particle diffusion [6-8]. The method is not only a numerical model for solving certain phenomena, but also a conceptual approach that does not only consider the analytical equations governing the phenomenon, but deals directly with the original phenomenon by means of an equivalent transmission line circuit. However, the simulation of an invisible cloaking with TLM is not straightforward since it faces at least two difficulties.

The first one lies in the fact that the cloak requires the use of metamaterials, i.e., materials with relative constitutive parameters below unity and, in particular regions, values very close to zero. So et al. have proposed a TLM node for modeling metamaterials [9], but the required modifications in the original TLM procedure to implement the later are significant. We report here a new condensed TLM node for two-dimensional (2D) problems having the merit of being almost identical to classical TLM approaches, simplification which constitutes another contribution of this paper.

The second difficulty comes from the anisotropic nature of the coating material since, in Cartesian coordinates, the permittivity and permeability tensors exhibit off-diagonal elements. The modeling of anisotropic materials in TLM was presented in [10,11], however the proposed implementation in both cases is very complicated, and, anyway would not permit the simulation of metamaterials easily. On the other hand, Huang et al. have proposed an approach that allows substituting the anisotropic medium by a concentric layered structure of alternating homogeneous isotropic materials [12]. The angular (ε_θ, μ_θ) and radial (ε_r, μ_r) components of the electromagnetic (EM) parameters can indeed be considered as the effective permittivity or permeability of a composite made of a series of parallel layers whenever the layers are thin enough compared with the wavelength. This statement is based on the effective medium theory and permits to match the angular components with the upper Wiener bound, while the radial components correspond to the lower Wiener bound [13]. We will choose this technique to carry out the cloaking simulation.

Finally Xi, et al., have recently studied the effect of different coordinate transformations [14]. They pointed out that the performance of the non-ideal cloak may be improved by using a nonlinear transformation instead of the well-established linear one.

Thus, many developments have been proposed lately, but the connection between them is not so clear. Therefore, beyond the TLM simulation of a cloaking itself, the aim of this article is to offer a comprehensive study in order to provide some clarification.

2. TLM topics

2.1 TLM Fundamentals

Since its formulation by Johns in 1971 [15], TLM has been extensively used to solve a wide range of electromagnetic problems. TLM is based on the analogy existing between Maxwell and transmission line equations to substitute the actual problem by an analogous one. In this manner, a certain medium and its electromagnetic properties are modeled by substituting the field space by a network of transmission lines formed by interconnecting unitary circuits or cells, called TLM nodes, which renders the problem discrete in space and time. Voltage pulses propagate from node to node in the mesh so constituted in substitution of the actual electric and magnetic fields in the original medium. As an example, the simplest node, constituted by four lines, is shown in Fig. 1. With such a simple node, it is only possible to model a 2D homogeneous material; more lines being required for inhomogeneous materials and/or 3-dimensional (3D) media [16]. Time is discretized through the iterative occurrence of two events at each time step. First, voltage pulses are incident upon the center of the nodes where they are scattered in all directions. All the information of this step is contained in the scattering matrix S which is previously determined to define Maxwell equations together with charge and energy conservation conditions. Secondly, the reflected voltage pulses generated in this way are transmitted to the adjacent nodes and become the incident pulses for the next time step. In this time-marching process, the value of the voltage pulses is known everywhere in the mesh at every time step, allowing the determination of the electric and magnetic fields or any related quantity. It is worth noting that most TLM nodes reported in the literature are of the condensed type in the sense that all the components of the EM field are available at the same point and at the same time, conferring to TLM a certain advantage over other methods, like Finite Differences in the Time Domain (FDTD), for example, which defines each component at slightly different points and even at a different time. This feature is especially interesting in modeling non-homogeneous media in which the imposition of precise transitions between different homogeneous materials is required [7].

 

Fig. 1. (Color online) Representation of the 2-dimensional TLM shunt node without stubs.

Download Full Size | PDF

2.2 Modeling of metamaterials with TLM

In its original form, a two-dimensional (2D) or three-dimensional (3D) TLM node is not able to model a metamaterial with anomalous permittivity or permeability values without including some modifications. In the case of metamaterials with relative permittivity or permeability below unity, the classical node can be used, provided the time step is reduced enough so as to adjust the increasing speed of light. This solution is satisfactory when modeling an isolated metamaterial but becomes problematic if a metamaterial with very low positive values of ε_r or μ_r coexists with usual materials such as vacuum. The first attempt to model metamaterials was reported by So, et al., [9] who modified the classical 2D shunt node by adding 5 extra capacitive and inductive stubs to define negative permeability and permittivity values. However, the node represents a significant modification of the impulse scattering process when compared to the classical simple version. In this paper, we will show that in fact this successful but drastic modification is not strictly necessary, it will suffice to design a simple 2D condensed node with only three new ports, a capacitive stub and two inductive stubs, to model materials with negative parameters or, what is more interesting for us, to model ε_r or μ_r values below unity without requiring a decrease in the time step value.

A 2D shunt node allows the simulation of TE modes with the electric field transverse to the propagation plane (y-oriented in this work), but the shunt node usually found in the literature represents a very particular case since it has identical dimensions along the three Cartesian directions, with identical admittance for all the link lines and it does not include magnetic losses. This particular version of 2D nodes is explained by the fact that most of the TLM effort has been concentrated on the design of general 3D condensed nodes. It has been only in later years that general condensed and lossy condensed nodes in 2D problems have been reported in the literature [17, 18], independently defining permittivity, permeability, together with electric and magnetic losses in a node with arbitrary length along each Cartesian direction. In this situation, we propose using the technique described in [18, 19] to design a 2D condensed node with a minimum number of inductive and capacitive stubs appropriately connected so as to model not only usual materials but also metamaterials with negative permittivity and/or permeability values.

In the classical node, the characteristic impedance Z ₀, or admittance Y ₀=1/Z ₀ of the main node lines are chosen to model the free space permittivity, ε ₀, and permeability, μ ₀. On the one hand, variations on ε_y with respect to the vacuum value are allowed by adding capacitance to the node. This can be achieved by connecting an open extra port, i.e., a capacitive stub, with admittance Y_yY ₀. On the other hand, variations on μ_x and μ_z require adding inductance which can be achieved by connecting two short-circuited or inductive stubs, with normalized impedance Z_xZ ₀ and Z_zZ ₀, respectively. The usual node finally obtained is shown in Fig. 2(a) and the expressions to determine Y_y, Z_x, and Z_z are:

where Δx, Δy, and Δz represent the length of the node along the x-, y-, and z-directions respectively.

The node proposed in this paper to model negative permittivity or permeability is depicted in Fig. 2(b). It differs only in the substitution of inductive stubs by capacitive stubs and vice versa. The scheme resembles that reported in [9], but the five extra stubs are substituted by only three new ports located at the node center, i.e., it is a condensed node, with similar connections to those of the classical 2D shunt node. The inductive stub, line 5, introduces in the node an inductance

 

Fig. 2. (Color online) Representation of (a) the classical 2D condensed shunt node with stubs, and (b) the modified 2D condensed shunt node allowing the simulation of metamaterials.

Download Full Size | PDF

which is equivalent to a frequency-dependent negative capacitance

Adding C_eq to the capacitance of lines 1 to 4 allows a reduction of the total ε_y below the vacuum value or even to achieve negative permittivity. Concretely, the equation defining ε_y becomes

Equivalently, the negative inductance introduced at frequency ω by the series connected capacitive ports, lines 6 and 7, yields

which again allows modeling an anomalous relative value of the magnetic permeability below unity or zero. It is worth noting that Eqs. (4) differ only by the factor -Δt ² ω ²/ 4 compared to Eqs. (1), but connectivity of the lines in Figs. 2 are identical, i.e., the scattering matrix, is unaltered with respect to the classical node. Relative permittivity or permeability less than 1 are now easily reachable for a given frequency and, in particular, the zero value which becomes a natural value with this approach without requiring a zero value of Δt.

Summarizing, the set of Eqs. (4) constitutes the only modification of the original 2D shunt node, the scattering matrix S remaining absolutely unchanged. Concretely, if the EM parameters are greater than 1 in a certain region of the space, Eqs. (1) have to be chosen while Eqs. (4) have to be used if the EM relative parameters take values below unity.

3. Cloaking structures modeling

Let consider a plane wave incident upon a Perfectly Electric Conducting (PEC) infinite cylinder of radius R ₁. The incident wave is a 2 GHz TE_y polarized plane wave with electric field parallel to the axis of the cylinder. The incident wave propagation vector, k, is x-oriented while the magnetic field is z-oriented.

The PEC cylinder is surrounded by a cloaking material with outer radius R ₂. According to the coordinate transformation, the anisotropic relative permittivity and permeability of the cloaking material in cylindrical coordinates have the following radius dependence:

Moreover, since the polarization is TE, the only constitutive parameters of interest are ε_y, μ_r, and μ_θ.

It has been shown that the anisotropic cloaking shell can be mimicked by concentric alternating layers made up of two homogeneous isotropic sub-layers A and B [12], chosen in this study to have the same thickness. The permittivity of each one is thus given by Eq. (5a), while the permeability μ_A, and μ_B can be obtained from Wiener formula:

which may be combined with Eqs. (5b, 5c) to give

The 2D TLM shunt node presented in section 2 is used to model this problem and obtain the electric field distribution. A TLM mesh formed with 2800 nodes along the x- and z-directions is used for the modeling. The node size is Δx=Δz=5.10^-4 m and the maximum allowable time step for these lengths, Δt=1.1785ps, is used. Typically, the boundary conditions used in cloaking studies are chosen to be absorbent along the direction of propagation and Perfectly Magnetic Conducting (PMC) for vertical directions, i.e., along the z-direction in this work, which is equivalent to artificially simulate a periodic structure of cylinders along the z-direction. However, the cloaking we are numerically considering is not ideal, and so a certain scattering is expected at each cylinder. In this manner, if we want to characterize the cloaking performance of the coating for a given cylinder, the existence of neighbor cylinders in the periodic structure would artificially magnify the non-ideality of the results. In other words, we are interested in the study of the capability of a single coating structure to make a certain region invisible, and deviations of this ideal goal are better identified if only one coating element are considered. In this sense, the PMC boundaries are not a satisfactory choice but, on the other hand, they cannot be directly substituted by absorbent boundaries because the plane wave condition corresponding to excitation would be broken. The problem is worked out by employing the Huygens surface technique which consists of dividing the mesh into an inner total field (incident and scattered) region and an outer scattered field region. The source conditions for the plane wave excitation are enforced on the interface separating these regions, which presents two benefits: first, absorbing boundary conditions can be directly applied to the limits of the mesh, since excitation has been removed at these points; and second, an illustrative map of the scattered field is directly obtained at the outer region. Finally, the far field is obtained in all cases using the transformation proposed in [20]. The technique considers the tangential EM fields on an imaginary surface surrounding the cylinder in the scattered field region as equivalent electric and magnetic sources. A numerical integration together with a frequency-dependent amplitude transformation which allows using a 3D scheme on a 2D problem provides the far field response.

 

Fig. 3. (a). (Color online) Far field pattern for unprotected cylinder and two different cloakings. (b) Total and scattered amplitude field plot at the vicinity of the coated cylinder for μ and ε sampled at the inner boundary of each layer.

Download Full Size | PDF

For the numerical simulation, the cloaking shell is made up of 10 pairs of layers A-B (in total, 20 sub-layers) and R ₁=0.1 m while R ₂=0.2 m, so that the thin layer condition is satisfied [12]. In the first modeling, ε_y, μ_r, and μ_θ are calculated at the center of the isotropic layers (i.e., on the interface separating layer A and layer B), the corresponding far field pattern is depicted in Fig. 3(a). For comparison, this field has been normalized by the maximum value of the field when no cloaking structure is present. This reference far field pattern, which is checked to perfectly coincide with the theoretical one [21], can be observed in the same figure. A significant reduction of the scattering is noticed in almost all directions at the notable exception of the forward scattering. It is known that the forward scattering of a PEC cylinder is more intense for TE mode compared to the TM mode case [21], so we are dealing with the worst case polarization. By judiciously adjusting the simulation parameters, we will show below that this scattering reduction can be enhanced.

Instead of calculating the parameters ε_y, μ_r, and μ_θ at the center of the isotropic layers, we now calculate them at the boundary between each layer. However, this process gives an infinite value of μ_θ for the first layer (i.e. the layer directly adjacent to the inner PEC cylinder) and, therefore, the electromagnetic parameters will not be calculated exactly at the inner boundary of the layers but at a distance 10^-4 δ of the later, where δ is the thickness of one layer. The far field pattern is displayed in Fig. 3(a) and, compared to the previous case, showing an even stronger scattering reduction. This reduction now occurs for all the direction and the forward scattering has been in particular reduced for 4.3 dB, while the backward scattering reduces for 15.1 dB. This confirms that the cloaking performance is highly sensitive to the parameters’ value of the inner layer; indeed, it is the most affected layer by this change. The simulated electric field is plotted in Fig. 3(b). It is straightaway clear that, once in the cloaking shell, the wave front is deviated and conducted around the cloaked cylinder, while outside the object under consideration it emerges almost unperturbed. Moreover, the Huygens technique allows visualizing in the same picture the scattered field at the outer almost green square region, with the incident field directly filtered out, that ostensibly takes the form of concentric waves. As it has been discussed in [22], the presence of the concentric scattering is mainly due to the absence of the magnetic surface current on the inner boundary of the cloak. Indeed, for an ideal cloaking, the permittivity and permeability for r=R ₁ is enforced by Eqs. (5) to be ε_y→0, μ_r→0, and μ_θ→∞. These extreme values induce a magnetic surface current exactly located on the inner boundary of the cloak once illuminated by the incident wave. Such a current shields the cloaked region and then avoids the field to enter this area. On the other hand, the inner layer μ_θ is not infinite in our simulation, so this surface current does not exist here and, then, scattering occurs on the interface of the PEC cylinder, making the cloaking structure of Fig. 3(b) imperfect.

4. Effect of using reduced parameters on the quality of the cloak

In order to facilitate experimental realization or in order to avoid infinite values for the electromagnetic constants, reduced material parameters are often used for the cloaking shell. These reduced parameters usually consist of choosing one of the constants to be 1, the penalty being only the loss of the reflectionless property as claimed in some papers [4, 5]. Nevertheless, Yan, et al., have showed that using reduced parameters experiments more scattering than intended and is more than just nonzero reflectance [23]. To verify this statement for the proposed layered cloaking structure, the far field calculation is calculated for two types of reduced parameters and displayed in Fig. 4(a). One corresponds to μ_θ=1, while the other one to ε_y=1 (the choice μ_r is not allowed since the corresponding ε_y would not be spatially uniform), both with parameters calculated at the center of the layers. As expected, the performance of the cloaking is deteriorated; however, it is worth noting that in both cases the forward scattering is now reduced especially for the choice ε_y=1 (for 7dB). Finally, it is worth noting from Fig. 4(a) that the greater the reduction in the forward direction, the lesser the reduction in the backward direction, which seems to be associated to energy conservation.

5. Effect of non-linear transformation on the quality of the cloak

Similarly, it is of interest to study the effect of a nonlinear transformation on the effectiveness of the proposed layered cloaking structure. A nonlinear transformation is characterized by a factor α that changes the electromagnetic parameters of Eqs. (5a), (5b) and (5c) ([14]. The far field pattern, for α=1 (corresponding to a linear transformation), α=0.8, and α=0.5, is depicted in Fig. 4(b) for the full parameters of Eqs. (5a), (5b) and (5c) calculated at the center of the layers. It can be observed that the choice α=0.8 globally improves the cloaking performance, particularly for angle contained between 90° and 270° in the backward direction. Concerning the α=0.5 transformation, it appears that the result is not convincing although it offers some improvement for some directions.

 

Fig. 4. (Color Online) (a) Far field pattern for reduced parameters. (b) Far field pattern for nonlinear transformation.

Download Full Size | PDF

6. Frequency dispersion

Schurig, et al., have presented experimental results demonstrating that a cloaking structure can be constructed in the microwave regime from a metamaterial consisting of split-ring resonators [4]. On the other hand, it is well-known that a dual distributed L-C network is able to reproduce the exotic parameters of a metamaterial which are no longer constant, and are, instead, explicit functions of frequency. No doubt that a cloaking structure could also be built with such a network. In this sense, TLM would be able to provide interesting information since a dual L-C network immersed in a classical TLM mesh is precisely the scheme used for the node shown in Fig. 2(b) proposed in this work. This means that the TLM solution is expected to reproduce the behavior of the cloaking material at the design frequency of 2 GHz, but a simple Fourier Transform of the time domain TLM results for a wide-band incident field would also provide realistic information on the behavior at other frequencies. Let us consider, the most efficient cloaking of Fig. 3 obtained with the electromagnetic parameters calculated at the inner boundary of each layers. By Fourier transform of the response to a Gaussian-shaped incident field, information at other frequencies can be directly obtained from a time-domain simulation. Figure 5 shows the 2D backscattering width normalized to λ, σ_2D/λ, for the cloaking device compared to that of the simple conducting cylinder. As expected, 2 GHz is the frequency for which the effectiveness of the cloaking material is optimum. However, the resonant nature of dual L-C networks causes that regions of reduced and magnified radiation alternate on both sides of the 2 GHz design frequency.

 

Fig. 5. (Color online) Two-dimensional scattering width, σ2D/λ, versus frequency of a conducting cylinder and cloaked cylinder at 2 GHZ.

Download Full Size | PDF

7. Conclusion

A TLM simulation of cloaking structures, made up of alternating isotropic layers mimicking an anisotropic shell, has been proposed in this paper. Previously, we have implemented a new elegant technique to simulate metamaterial with TLM. This new approach lets the usual TLM procedure almost unchanged; the scattering matrix is in particular exactly the same, which provides a lot of flexibility in the programming task. We have therefore tackled the cloaking modeling, seizing all the possibilities offered by a time-domain method. In this way, the Huygens surface technique has allowed using more appropriate boundary conditions avoiding the artificial simulation of an infinite quantity of periodic cylinders. Furthermore, this technique allows, in a single simulation, to obtain total field zones and scattering field zones. We have shown that it is possible to improve the cloaking effectiveness by judiciously choosing the electromagnetic parameters of the shell. In this sense, it has been established that the parameters of each composite layer should be calculated at the inner boundary of each one to obtain better cloaking performances. Finally, the TLM mesh being equivalent to an L-C network seems to be extremely well suited to predict the behavior of a cloak material which would be made up of embedded lumped capacitors and inductors. The time domain response to a Gaussian pulse shows that an operating frequency range exists, with alternating regions of reduced and amplified radiation.