Finite Element Analysis of Electromagnetic Waves in Two-Dimensional Transformed Bianisotropic Media

Finite Element Analysis of Electromagnetic Waves in Two-Dimensional Transformed Bianisotropic Media

Yan Liu yanliu@xidian.edu.cn Boris Gralak Sebastien Guenneau School of Aerospace Science and Technology, Xidian University, Xi’an 710071, China Aix-Marseille Université, CNRS, Ecole Centrale Marseille, Institut Fresnel
13013 Marseille, France
Abstract

We analyse wave propagation in two-dimensional bianisotropic media with the Finite Element Method (FEM). We start from the Maxwell-Tellegen’s equations in bianisotropic media, and derive some system of coupled Partial Difference Equations (PDEs) for longitudinal electric and magnetic field components. Perfectly Matched Layers (PMLs) are discussed to model such unbounded media. We implement these PDEs and PMLs in a finite element software. We apply transformation optics in order to design some bianisotropic media with interesting functionalities, such as cloaks, concentrators and rotators. We propose a design of metamaterial with concentric layers made of homogeneous media with isotropic permittivity, permeability and magneto-electric parameters that mimic the required effective anisotropic tensors of a bianisotropic cloak in the long wavelength limit (homogenization approach). Our numerical results show that well-known metamaterials can be transposed to bianisotropic media.

keywords:
Finite Element Method, Bianisotropic media, Transformation Optics, Cloak, Concentrator, Rotator
Msc:
[2010] 00-01, 99-00
journal: Journal of LaTeX Templates

1 Introduction

In the last decade, metamaterials have attracted much attention due to their extraordinary properties, such as negative refraction, ultra refraction, anomalous dispersion and so on. Metamaterials are artificial materials engineered to gain desired properties that cannot be found in nature and usually consist of periodically arranged materials, which affect electromagnetic waves in an unconventional manner. More precisely, they exhibit new and unusual electromagnetic properties at the macro-scale, due to their structural features smaller than the operational wavelength of the electromagnetic wave. Metamaterials include the negative index materials (NIMs), single negative metamaterials, bi-isotropic and bianisotropic media, chiral media and so on.

As an application, the physicist John Pendry proposed that NIMs enable a perfect lens Pendry00 (), which allows a sub-wavelength imaging. On the other hand, J. Pendry et al in 2006 Pendry06 (); Leonhardt06 () proposed that a geometric transformation of space could distort light trajectories around a bounded region, which is thus made invisible. This powerful mathematical technique is called transformation optics, and it presents great potential thanks to the fabrication of metamaterials with specially designed properties Rahm08 (); Rahm2008 (); Chen07 (). Transformation optics has been successfully applied to design functionalities, such as invisibility cloaks, rotators, concentrators etc.

In NIMs, or single negative metamaterials (e.g. composites with or ), we have supposed that the metamaterials have independent electric and magnetic responses described by and . However, magnetoelectric coupling does occur in many other electromagnetic metamaterials, wherein the electric and magnetic fields are induced reciprocally. Anisotropic (resp. isotropic) media with such kind of property are referred as bianisotropic (resp. bi-isotropic) Marques02 (); Rill09 ().

In yanl13 (); yanl14 (), we have theoretically proved that the Pendry-Ramakrishna lens theorem is applicable to complementary bianisotropic media; on the other hand, as a powerful mathematical technique - transformation optics allows for a transformation of space between two different coordinate systems, while the Maxwell-Tellegen’s equations for bianisotropic media are also proved to be form invariant under space transformation. Although it is easy to understand the principle used to design novel devices such as cloaks, concentrators, rotators and so on, a numerical model is required to illustrate their electromagnetic response. Hence, in this paper, we would like to discuss in details the implementation of transformed bianisotropic media in a finite element model, and further show numerical results for an invisibility cloak, concentrator and rotator, the former being even achieved with a simplified multilayered design.

We know that, Partial Differential Equations (PDEs) are generally used to model many physical phenomena such as fluid dynamics and electromagnetism, etc. In complex media, solutions to the governing equations can be difficult to derive in closed-form by traditional analytical routes - Fourier or Laplace transform methods, power series expansion and so on, hence one often has to resort to numerical approximation of the solutions. As a particular class of numerical techniques, FEM is efficient to solve problems in heterogeneous anisotropic media Courant43 (); Zlamal68 (), and it has been widely applied in engineering design and analysis in mechanics starting from the late seventies, while researchers in photonics started to use it in the nineties. Notably, the COMSOL Multiphysics package is a finite element analysis software much used nowadays to solve various physics and engineering applications in the metamaterials’ community, as it allows entering coupled systems of PDEs, such as in opto-mechanics, thermal elasticity etc. Here, we make use of it to set up a coupled PDE system for a bianisotropic structure. To do this, we start from the Maxwell-Tellegen’s equations in bianisotropic media and derive two coupled PDEs with the longitudinal electric and magnetic fields as the unknowns; then we discuss open boundary conditions known in the photonics literature as perfectly matched layers (PMLs). To illustrate the usefulness of our numerical algorithm, we apply transformation optics to bianisotropic media in order to design some functional devices, i.e. invisibility cloak, concentrator and rotator. Their electromagnetic (EM) properties studied thanks to our PDEs model implemented in COMSOL Multiphysics prove that these well-known devices work equally well in bianisotropic media.

2 A coupled PDE sytem for bianisotropic media

We first recall the time-harmonic Maxwell’s equations (assuming a time dependence with the time variable and the angular wave frequency)

(1)

where and are the electric and magnetic field respectively, while and are the electric displacement and magnetic flux density.

Let us assume the constitutive relations in a bianisotropic medium described by:

(2)

where is the permittivity, is the permeability, and are the magnetoelectric coupling parameters. All these material parameters can be treated as rank-2 tensors. We now substitute (2) into (1) and obtain the Maxwell-Tellegen’s equations

(3)

Let us now consider an orthogonal coordinate system , and assume that the material parameters and the electromagnetic field are invariant in one direction, say, the direction, that is, , , , , are independent of .

Let us define , , and further assume following the terminology used in the book by F. Zolla et al. Zolla12 (), that , , and are -anisotropic tensors i.e. such that

(4)

We then apply (4) to (3), and have

(5)
(6)
(7)
(8)
(9)
(10)

We rewrite (5) and (6) in a matrix form

(11)

Similarly, (8) and (9) become

(12)

If we define

(13)

and

(14)

then (11) and (12) turn out to be

(15)
(16)

Furthermore, from (15), we have

(17)

Substituting it into (16), we obtain

(18)

which can be recast as

(19)

with

(20)

Let us now plug (19) into (10)

(21)

In the same way, from (16), we have

(22)

which is then applied to (15)

(23)

where

(24)

with

(25)

We plug (24) into (7)

(26)

Finally, we rewrite (21) and (26) by multiplying by on both sides of the equations, and we obtain

(27)

which are two coupled equations with the longitudinal fields and as unknowns.

3 Implementation of the PDE system in Comsol Multiphysics

The standard form of PDEs in COMSOL Multiphysics is described by the following equation reminiscent of a governing equation in elasticity theory:

(28)

where is a rank-4 tensor, is a rank -tensor, and are vectors and is a forcing term.

We rewrite the formulae in (27) as

(29)

Comparing with (28), if we consider the coupled longitudinal electric and magnetic fields, the unknown in (28) is

(30)

and the coefficients are

(31)

with entries

(32)

Simply stated, when the tensors of in (13) are represented by diagonal matrices in a Cartesian coordinate system

(33)

then the entries in (32) become

(34)

Note that each component of the coefficients is a fraction with a denominator . In order to avoid any infinite entries of while solving the PDEs, it should satisfy

(35)

Remark: This condition is essential and generally is met since the magnetoelectric coupling parameters for bianisotropic media are usually small.

4 Problem setup with open boundary conditions

In the last section, we have discussed the general PDEs in bianisotropic media, which allows us to introduce the PDEs model in COMSOL Multiphysics; however, in order to mimic the scattering and propagation properties of bianisotropic media in an open space, we need to consider the boundary conditions. We usually introduce a transformation of an infinite domain into a finite one, which should enforce that the wavelength is contracted to infinitely small values as it approaches the outer boundary of the transformed domain within a perfectly matched layer (PML). At this outer boundary, Dirichlet data could be set, thereby enforcing a vanishing field Berenger94 (); Teixeira98 ().

The PML is shown to be equivalent to an analytic continuation of Maxwell’s equations to complex variables’s spatial domain Teixeira97 (). The form of the electromagnetic field inside a PML region can be obtained from the solutions in real space from a mapping to a complex space through this analytic continuation. Physically speaking, an equivalent artificial medium is achieved from this map in the PML region that is described by complex, anisotropic and inhomogeneous parameters, even if the original ones are real, isotropic and homogeneous. The analytic continuation means to alter the eigenfunctions of Maxwell’s equations in such a way that the propagating modes are mapped continuously to exponentially decaying modes, which allows for reflectionless absorption of electromagnetic waves. In other words, an absorbing medium dissipating the outgoing wave can be achieved through proper complex map: an artificial medium situated in a region can be added to a finite medium in order to mimic an open space within which the field is damped to a negligible value Agha08 (). Importantly, the impedance of the equivalent medium is the same as that of the initial medium since all the parameters undergo the same transform, what ensures non-reflecting features on the interface between the medium and the PMLs.

F. L. Teixeira et.al Teixeira97 (); Teixeira98 () derive the Maxwellian PML’s for arbitrary bianisotropic and dispersive media. In the Cartesian coordinate system, if an analytic continuation is defined by the transformation function

(36)

with the complex stretching variables Chew94 () and stands for , , . In the complex space, the nabla operator in Maxwell’s equations changes as

(37)

it reads as

(38)

with

(39)

since and commute for , is a diagonal tensor, and .

Based on this, we then can expand the Maxwell’s equations in the new complex space, and derive that the PMLs with bianisotropic constitutive parameters are

(40)
Figure 1: Computational domain enclosed by perfectly matched layers (PMLs) with functionality in the , and - direction Rahm08 ().

Take a full wave simulation for two dimensional structures in COMSOL Multiphysics as an illustration, Fig. 1 shows the computational domain terminated by PMLs in Cartesian coordinates. The PMLs are designed to decrease the distribution of the waves along the , and - directions, respectively. Specifically, assuming the parameter of the matrix (), then we choose

(41)

with . Then the transformed parameters in the PMLs are

(42)

where

(43)

can be derived by introducing the definition of , and in (41) along different directions. Similar concept of PMLs can be extended to the anisotropic/bianisotropic matrix and other coordinate system Teixeira97 (); Teixeira1998 ().

5 FEM implementation

To check both our PDEs model and PMLs for bianisotropic media, some conceptual devices designed from transformation optics are numerically studied in the following subsections, these structures are all achieved by bianisotropic media.

5.1 Invisibility cloak

Electromagnetic (EM) metamaterials such as invisibility cloaks can be designed through the blowup of a point Leonhardt06 (); Pendry06 (), using transformation optics which is a versatile mathematical tool enabling a deeper analytical insight into the scattering properties of EM fields in metamaterials. In this subsection, a cloak consisting of a bianisotropic medium is investigated, and the geometric transformation proposed by J. Pendry Pendry06 () is proved to be applicable to design such kind of cloak.

We start with a free space with permittivity , permeability and magnetoelectric coupling parameters , in the Cartesian coordinate system, and first map it onto polar coordinates as defined by

(44)

a disk with results as shown in Fig. 2(a). Let us then introduce a geometric transformation defined by (45) which maps the field within this disk onto an annulus , i.e. the original point (red) in the original disk (Fig. 2(a)) is blowup to a disk with in coordinates as shown in panel (b) of Fig. 2. Based on this coordinate transformation, the schematic diagram of the ray trajectories in the cloak is shown in panel (c): the rays are distorted and guided around the cloak.

Figure 2: Geometric transformation for invisibility cloak: (a) Disk (grey) with radius and the original point (red) is at ; (b) The original point is mapped to a disk with , i.e the disk with is compressed to an annulus , the mapping function is defined by (45); (c) Ray trajectories in the cloak.

The mapping function is defined as Pendry06 ()

(45)

which provides the Jacobian matrix

(46)

Finally, we go back to the Cartesian coordinates , which are radially contracted, and the compound Jacobian matrix is

(47)

On the other hand, according to the coordinates transformation (44), we have

(48)

with

(49)

Then, (47) becomes

(50)

The transformation matrix is correspondingly

(51)

More explicitly, we apply the matrix into (51) which we denote by

(52)

with

(53)

Hence, the parameters in the annulus are

(54)

while the parameters of the outer region are unchanged, and that of the disk can have any value without affecting the electromagnetic scattering.

Figure 3: (a) Schematic diagram of bianisotropic cloak with a L-shaped obstacle inside the cloaking region, and a point source (-polarized with the electric field along with Hz) locates outside the shell. Plots of Re: (b) a point source in a pure matrix with parameters , and ; (c) same point source in the matrix with a presence of L-shaped obstacle, the parameters of obstacle are , and ; (d) same point source that radiates in the presence of an L-shaped obstacle surrounded by an invisibility cloak.

Fig. 3(a) shows our designed bianisotropic cloak, the cloaking region is the innermost circle of radius um wherein a L-shaped obstacle locates, the grey annulus represents the distorted space with radius , where um. Firstly, the EM distribution of a point source (-polarized with the electric field along ) with frequency Hz radiating in a matrix with parameters is shown in Fig. 3(b). are the permittivity, permeability of the vacuum, respectively; while with the velocity of light in vacuum should be promised to ensure convergence of the numerical algorithm as indicated in (35). If there is a L-shaped obstacle with , , in the matrix, it leads to an interacting with the source, the plot of Re is indicated in panel (c). However, when the L-shaped obstacle is surrounded by the designed cloak, then the obstacle seems to be invisible for the outer observer, as shown in panel (d).

5.2 Concentrator

A cylindrical concentrator is designed to focus the incident waves with wave vectors perpendicular to the cylinder axis, enhancing the electromagnetic energy density in a given area Rahm08 ().

We start with a map from Cartesian coordinates to cylindrical ones, a cylindrical lens with is shown in Fig. 4(a), then a radial transformation is introduced to distort the space as

(55)
Figure 4: (a) Schematic diagram of a concentrator: a space is compressed into a circle with radius (the innermost in red) with the expense of an expansion of space (grey color) between and ; the dotted circle corresponds to a virtual interface used in the geometric transformation, see (55). (b) Ray trajectories in the concentrator.

As shown in Fig. 4(a), a space is compressed into a cylindrical region with radius (the innermost circle) with the expense of an expansion of space(grey color) between and , where an intermediate circle is located at . Fig. 4(b) shows the ray trajectories in such a concentrator. The associated Jacobian matrix for the transformation in (55) is

(56)

Finally we go back to Cartesian coordinates. The compound Jacobian matrix is then

(57)

Furthermore, the transformation matrix is

(58)

Let us expand this formula in the distorted cylindrical coordinates , and denote as in (52), finally we have

(59)

with

(60)

Apply matrix to (54), we can obtain the formulae of these parameters in the annulus as well as the innermost circle, while the outside space of the concentrator is bianisotropic matrix with .

Fig. 5(a) is the schematic diagram of our concentrator, the radii of the inner and outer circles are um, um, while the intermediate circle has a radius um. The parameters of each region are defined in (54) along with transfer matrix in (59), here the magnetoelectric coupling parameters are chosen as the same values as those of the invisibility cloak. An -polarized (the electric field is perpendicular to the - plane) plane wave with frequency Hz is incident from above, the plot of Re is depicted in panel (b): the waves are compressed in the inner circle as expected.

Figure 5: (a) Schematic diagram of a bianisotropic concentrator. (b) Plot of Re under an -polarized incidence with frequency Hz, the parameters in the region , and the annulus are given by (54) with transformation matrix in (59).

5.3 Rotator

Fig. 6 shows a rotator allowing a rotation of the incident wave, the distribution of the field inside a disk appears as it comes from a different angle comparing with that of outside , for example, a angle rotation is achieved here in the inner disk.

Figure 6: (a) Schematic diagram of the rotator: incident waves are rotated in the annulus and reach the disk of radius with a degree rotation. (b) Ray trajectories in rotator.

Be different from the geometric transformation introduced in the concentrator, an angular transformation is introduced in this case instead of the radial transformation. The mapping is performed in the region : the domain in is transformed by rotating a fixed angle from the region in virtual space, while the rotation angle is continually changed from the fixed angle to zero in region . In other words, the transformation algorithm can be expressed as

(61)

where is the rotation angle for the inner disk, and it is reduced to zero as the radius approaches to . is an arbitrary continuous function of .

Based on this transformation, we can derive the associated Jacobian matrix as

(62)

again, the compound Jacobian matrix can be obtained by , as well as the transformation matrix . In annulus , the notations in (52) become

(63)

with . In regions