Finite Element Analysis of Electromagnetic Waves in Two-Dimensional Transformed Bianisotropic Media
FFinite Element Analysis of Electromagnetic Waves inTwo-Dimensional Transformed Bianisotropic Media
Yan Liu a, ∗ , Boris Gralak b , Sebastien Guenneau b a School of Aerospace Science and Technology, Xidian University, Xi’an 710071, China b Aix-Marseille Université, CNRS, Ecole Centrale Marseille, Institut Fresnel13013 Marseille, France
Abstract
We analyse wave propagation in two-dimensional bianisotropic media with theFinite Element Method (FEM). We start from the Maxwell-Tellegen’s equationsin bianisotropic media, and derive some system of coupled Partial DifferenceEquations (PDEs) for longitudinal electric and magnetic field components. Per-fectly Matched Layers (PMLs) are discussed to model such unbounded media.We implement these PDEs and PMLs in a finite element software. We applytransformation optics in order to design some bianisotropic media with inter-esting functionalities, such as cloaks, concentrators and rotators. We proposea design of metamaterial with concentric layers made of homogeneous mediawith isotropic permittivity, permeability and magneto-electric parameters thatmimic the required effective anisotropic tensors of a bianisotropic cloak in thelong wavelength limit (homogenization approach). Our numerical results showthat well-known metamaterials can be transposed to bianisotropic media.
Keywords:
Finite Element Method, Bianisotropic media, TransformationOptics, Cloak, Concentrator, Rotator ∗ Corresponding author
Email address: [email protected] (Yan Liu)
Preprint submitted to Journal of L A TEX Templates December 18, 2015 a r X i v : . [ phy s i c s . op ti c s ] D ec . Introduction In the last decade, metamaterials have attracted much attention due to theirextraordinary properties, such as negative refraction, ultra refraction, anoma-lous dispersion and so on. Metamaterials are artificial materials engineered togain desired properties that cannot be found in nature and usually consist ofperiodically arranged materials, which affect electromagnetic waves in an uncon-ventional manner. More precisely, they exhibit new and unusual electromagneticproperties at the macro-scale, due to their structural features smaller than theoperational wavelength of the electromagnetic wave. Metamaterials include thenegative index materials (NIMs), single negative metamaterials, bi-isotropic andbianisotropic media, chiral media and so on.As an application, the physicist John Pendry proposed that NIMs enable aperfect lens [1], which allows a sub-wavelength imaging. On the other hand, J.Pendry et al in 2006 [2, 3] proposed that a geometric transformation of spacecould distort light trajectories around a bounded region, which is thus madeinvisible. This powerful mathematical technique is called transformation op-tics, and it presents great potential thanks to the fabrication of metamaterialswith specially designed properties [4, 5, 6]. Transformation optics has been suc-cessfully 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 andmagnetic responses described by ε and µ . However, magnetoelectric couplingdoes occur in many other electromagnetic metamaterials, wherein the electricand magnetic fields are induced reciprocally. Anisotropic (resp. isotropic) mediawith such kind of property are referred as bianisotropic (resp. bi-isotropic) [7, 8].In [9, 10], we have theoretically proved that the Pendry-Ramakrishna lenstheorem is applicable to complementary bianisotropic media; on the other hand,as a powerful mathematical technique - transformation optics allows for a trans-formation of space between two different coordinate systems, while the Maxwell-2ellegen’s equations for bianisotropic media are also proved to be form invariantunder space transformation. Although it is easy to understand the principle usedto design novel devices such as cloaks, concentrators, rotators and so on, a nu-merical model is required to illustrate their electromagnetic response. Hence, inthis paper, we would like to discuss in details the implementation of transformedbianisotropic media in a finite element model, and further show numerical re-sults for an invisibility cloak, concentrator and rotator, the former being evenachieved with a simplified multilayered design.We know that, Partial Differential Equations (PDEs) are generally used tomodel many physical phenomena such as fluid dynamics and electromagnetism,etc. In complex media, solutions to the governing equations can be difficult toderive in closed-form by traditional analytical routes - Fourier or Laplace trans-form methods, power series expansion and so on, hence one often has to resortto numerical approximation of the solutions. As a particular class of numeri-cal techniques, FEM is efficient to solve problems in heterogeneous anisotropicmedia [11, 12], and it has been widely applied in engineering design and analy-sis in mechanics starting from the late seventies, while researchers in photonicsstarted to use it in the nineties. Notably, the COMSOL Multiphysics package isa finite element analysis software much used nowadays to solve various physicsand engineering applications in the metamaterials’ community, as it allows en-tering coupled systems of PDEs, such as in opto-mechanics, thermal elasticityetc. Here, we make use of it to set up a coupled PDE system for a bianisotropicstructure. To do this, we start from the Maxwell-Tellegen’s equations in bian-isotropic media and derive two coupled PDEs with the longitudinal electric andmagnetic fields as the unknowns; then we discuss open boundary conditionsknown in the photonics literature as perfectly matched layers (PMLs). To illus-trate the usefulness of our numerical algorithm, we apply transformation opticsto bianisotropic media in order to design some functional devices, i.e. invisibilitycloak, concentrator and rotator. Their electromagnetic (EM) properties studiedthanks to our PDEs model implemented in COMSOL Multiphysics prove thatthese well-known devices work equally well in bianisotropic media.3 . A coupled PDE sytem for bianisotropic media We first recall the time-harmonic Maxwell’s equations (assuming a time de-pendence e − iωt with t the time variable and ω the angular wave frequency) ∇ × E = iω B , ∇ × H = − iω D . (1)where E and H are the electric and magnetic field respectively, while D and B are the electric displacement and magnetic flux density.Let us assume the constitutive relations in a bianisotropic medium describedby: D = ε E + iξ H , B = − iζ E + µ H . (2)where ε is the permittivity, µ is the permeability, ξ and ζ are the magnetoelectriccoupling parameters. All these material parameters can be treated as rank-2tensors. We now substitute (2) into (1) and obtain the Maxwell-Tellegen’sequations ∇ × E = ωζ E + iωµ H , ∇ × H = − iωε E + ωξ H . (3)Let us now consider an orthogonal coordinate system ( x , x , x ) , and assumethat the material parameters and the electromagnetic field are invariant in onedirection, say, the x direction, that is, ε , µ , ξ , E , H are independent of x .Let us define E = ( E , E , E ) T , H = ( H , H , H ) T , and further assumefollowing the terminology used in the book by F. Zolla et al. [13], that ε , µ , ξ and ζ are z -anisotropic tensors i.e. such that ε = ε ε ε ε
00 0 ε , µ = µ µ µ µ
00 0 µ ,ξ = ξ ξ ξ ξ
00 0 ξ , ζ = ζ ζ ζ ζ
00 0 ζ . (4)4e then apply (4) to (3), and have ∂ x E = ωζ E + ωζ E + iωµ H + iωµ H , (5) − ∂ x E = ωζ E + ωζ E + iωµ H + iωµ H , (6) ∂ x E − ∂ x E = ωζ E + iωµ H , (7) ∂ x H = − iωε E − iωε E + ωξ H + ωξ H , (8) − ∂ x H = − iωε E − iωε E + ωξ H + ωξ H , (9) ∂ x H − ∂ x H = − iωε E + ωξ H . (10)We rewrite (5) and (6) in a matrix form ∂ x E ∂ x E = ω − ζ ζ ζ − ζ E − E + iω − µ µ µ − µ H − H . (11)Similarly, (8) and (9) become ∂ x H ∂ x H = iω ε − ε − ε ε E − E + ω − ξ ξ ξ − ξ H − H . (12)If we define ε T = ε − ε − ε ε , µ T = − µ µ µ − µ ,ξ T = − ξ ξ ξ − ξ , ζ T = − ζ ζ ζ − ζ , (13)and A = ∂ x ∂ x , E = E − E , H = H − H . (14)then (11) and (12) turn out to be AE = ωζ T E + iωµ T H , (15)5 H = iωε T E + ωξ T H . (16)Furthermore, from (15), we have E = ω − ζ − T AE − iζ − T µ T H . (17)Substituting it into (16), we obtain AH = iωε T ( ω − ζ − T AE − iζ − T µ T H ) + ωξ T H = iε T ζ − T AE + ωε T ζ − T µ T H + ωξ T H . (18)which can be recast as H = B − (cid:0) AH − iε T ζ − T AE (cid:1) , (19)with B = ωξ T + ωε T ζ − T µ T . (20)Let us now plug (19) into (10) (cid:104) ∂ x ∂ x (cid:105) H = (cid:104) ∂ x ∂ x (cid:105) (cid:0) B − AH − iB − ε T ζ − T AE (cid:1) = (cid:104) ∂ x ∂ x (cid:105) (cid:26) − iω ( µ T + ζ T ε − T ξ T ) − (cid:27) AE + (cid:104) ∂ x ∂ x (cid:105) (cid:26) ω ( ξ T + ε T ζ − T µ T ) − (cid:27) AH = − iωε E + ωξ H . (21)In the same way, from (16), we have H = ω − ξ − T AH − iξ − T ε T E , (22)which is then applied to (15) AE = ωζ T E + iωµ T ( ω − ξ − T AH − iξ − T ε T E )= ωζ T E + iµ T ξ − T AH + ωµ T ξ − T ε T E , (23)where E = C − AE − iC − µ T ξ − T AH , (24)6ith C = ωζ T + ωµ T ξ − T ε T . (25)We plug (24) into (7) [ ∂ x ∂ x ] E = [ ∂ x ∂ x ] (cid:0) C − AE − iC − µ T ξ − T AH (cid:1) = [ ∂ x ∂ x ] (cid:26) ω ( ζ T + µ T ξ − T ε T ) − (cid:27) AE − [ ∂ x ∂ x ] (cid:26) iω ( ε T + ξ T µ − T ζ T ) − (cid:27) AH = ωζ E + iωµ H . (26)Finally, we rewrite (21) and (26) by multiplying by ω on both sides of theequations, and we obtain (cid:104) ∂ x ∂ x (cid:105) (cid:8) − i ( µ T + ζ T ε − T ξ T ) − (cid:9) (cid:104) ∂ x ∂ x (cid:105) T E + (cid:104) ∂ x ∂ x (cid:105) (cid:8) ( ξ T + ε T ζ − T µ T ) − (cid:9) (cid:104) ∂ x ∂ x (cid:105) T H = − iω ε E + ω ξ H , (cid:104) ∂ x ∂ x (cid:105) (cid:8) ( ζ T + µ T ξ − T ε T ) − (cid:9) (cid:104) ∂ x ∂ x (cid:105) T E + (cid:104) ∂ x ∂ x (cid:105) (cid:8) − i ( ε T + ξ T µ − T ζ T ) − (cid:9) (cid:104) ∂ x ∂ x (cid:105) T H = ω ζ E + iω µ H . (27)which are two coupled equations with the longitudinal fields E and H asunknowns.
3. Implementation of the PDE system in Comsol Multiphysics
The standard form of PDEs in COMSOL Multiphysics is described by thefollowing equation reminiscent of a governing equation in elasticity theory: ∇ · ( − c ∇ u − α u + γ ) + a u + β · ∇ u = f . (28)where c is a rank-4 tensor, α is a rank -tensor, γ and β are vectors and f is aforcing term. 7e rewrite the formulae in (27) as (cid:104) ∂ x ∂ x (cid:105) (cid:8) ( µ T + ζ T ε − T ξ T ) − (cid:9) (cid:104) ∂ x ∂ x (cid:105) T E + (cid:104) ∂ x ∂ x (cid:105) (cid:8) i ( ξ T + ε T ζ − T µ T ) − (cid:9) (cid:104) ∂ x ∂ x (cid:105) T H = ω ε E + iω ξ H , (cid:104) ∂ x ∂ x (cid:105) (cid:8) i ( ζ T + µ T ξ − T ε T ) − (cid:9) (cid:104) ∂ x ∂ x (cid:105) T E + (cid:104) ∂ x ∂ x (cid:105) (cid:8) ( ε T + ξ T µ − T ζ T ) − (cid:9) (cid:104) ∂ x ∂ x (cid:105) T H = iω ζ E − ω µ H . (29)Comparing with (28), if we consider the coupled longitudinal electric and mag-netic fields, the unknown u in (28) is u = E H . (30)and the coefficients are c = c c c c , a = ω ε iω ξ iω ζ − ω µ . (31)with entries c = ( µ T + ζ T ε − T ξ T ) − , c = i ( ξ T + ε T ζ − T µ T ) − ,c = i ( ζ T + µ T ξ − T ε T ) − , c = ( ε T + ξ T µ − T ζ T ) − . (32)Simply stated, when the tensors of ε, µ, ξ, ζ in (13) are represented by diagonalmatrices in a Cartesian coordinate system v = diag( ν , ν , ν ) , v = ε, µ, ξ, ζ . (33)8hen the entries in (32) become c = − ε µ ε − ξ ζ − ε µ ε − ξ ζ ,c = iζ µ ε − ξ ζ iζ µ ε − ξ ζ ,c = iξ µ ε − ξ ζ iξ µ ε − ξ ζ ,c = µ µ ε − ξ ζ µ µ ε − ξ ζ . (34)Note that each component of the coefficients c ij ( i, j = 1 , is a fraction with adenominator ( ε kk µ kk − ξ kk ζ kk ) ( k = 1 , . In order to avoid any infinite entriesof c while solving the PDEs, it should satisfy ε kk µ kk − ξ kk ζ kk (cid:54) = 0 , ( k = 1 , (35) Remark : This condition is essential and generally is met since the magneto-electric 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 me-dia, which allows us to introduce the PDEs model in COMSOL Multiphysics;however, in order to mimic the scattering and propagation properties of bian-isotropic 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 val-ues as it approaches the outer boundary of the transformed domain within aperfectly matched layer (PML). At this outer boundary, Dirichlet data could beset, thereby enforcing a vanishing field [14, 15].9he PML is shown to be equivalent to an analytic continuation of Maxwell’sequations to complex variables’s spatial domain [16]. The form of the electro-magnetic field inside a PML region can be obtained from the solutions in realspace from a mapping to a complex space through this analytic continuation.Physically speaking, an equivalent artificial medium is achieved from this mapin the PML region that is described by complex, anisotropic and inhomogeneousparameters, even if the original ones are real, isotropic and homogeneous. Theanalytic continuation means to alter the eigenfunctions of Maxwell’s equationsin such a way that the propagating modes are mapped continuously to exponen-tially decaying modes, which allows for reflectionless absorption of electromag-netic waves. In other words, an absorbing medium dissipating the outgoing wavecan be achieved through proper complex map: an artificial medium situated ina region can be added to a finite medium in order to mimic an open space withinwhich the field is damped to a negligible value [17]. Importantly, the impedanceof the equivalent medium is the same as that of the initial medium since all theparameters undergo the same transform, what ensures non-reflecting featureson the interface between the medium and the PMLs.F. L. Teixeira et.al [16, 15] derive the Maxwellian PML’s for arbitrary bian-isotropic and dispersive media. In the Cartesian coordinate system, if an ana-lytic continuation is defined by the transformation function ˜ u = (cid:90) u s u ( u (cid:48) ) du (cid:48) . (36)with s u ( u (cid:48) ) the complex stretching variables [18] and u stands for x , y , z . Inthe complex space, the nabla operator in Maxwell’s equations changes as ∇ → ˜ ∇ = ˆ x ∂∂ ˜ x + ˆ y ∂∂ ˜ y + ˆ z ∂∂ ˜ z = ˆ x s x ∂∂x + ˆ y s y ∂∂y + ˆ z s z ∂∂z , (37)it reads as ˜ ∇ = = S · ∇ , (38)with = S = ˆ x ˆ x ( 1 s x ) + ˆ y ˆ y ( 1 s y ) + ˆ z ˆ z ( 1 s z ) . (39)10ince s u ( u ) and ∂/∂u (cid:48) commute for u (cid:54) = u (cid:48) , = S is a diagonal tensor, and det( = S ) =( s x s y s z ) − .Based on this, we then can expand the Maxwell’s equations in the new com-plex space, and derive that the PMLs with bianisotropic constitutive parametersare v PML = (det = S ) − [ = S · v · = S ] , v = ε, µ, ξ, ζ . (40) object matrix PML y PML x PML xy PML y PML x PML xy PML xy PML xy Figure 1: Computational domain enclosed by perfectly matched layers (PMLs) with function-ality in the x , y and x - y direction [4]. Take a full wave simulation for two dimensional structures in COMSOLMultiphysics as an illustration, Fig. 1 shows the computational domain ter-minated by PMLs in Cartesian coordinates. The PMLs are designed to de-crease the distribution of the waves along the x , y and x - y directions, respec-tively. Specifically, assuming the parameter of the matrix v = diag( v xx , v yy , v zz ) ( v = ε, µ, ξ, ζ ), then we choose = S PML x : s x = a − b ∗ i, s y = 1 , s z = 1 , PML y : s x = 1 , s y = a − b ∗ i, s z = 1 , PML xy : s x = a − b ∗ i, s y = a − b ∗ i, s z = 1 . (41)with a, b ∈ Z + . Then the transformed parameters in the PMLs are v (cid:48) = diag( v xx L xx , v yy L yy , v zz L zz ) , v = ε, µ, ξ, ζ . (42)11here L xx = s y s z s x , L yy = s x s z s y , L zz = s x s y s z . (43)can be derived by introducing the definition of s x , s y and s z in (41) along differ-ent directions. Similar concept of PMLs can be extended to the anisotropic/bianisotropicmatrix and other coordinate system [16, 19].
5. FEM implementation
To check both our PDEs model and PMLs for bianisotropic media, someconceptual devices designed from transformation optics are numerically studiedin the following subsections, these structures are all achieved by bianisotropicmedia.
Electromagnetic (EM) metamaterials such as invisibility cloaks can be de-signed through the blowup of a point [3, 2], using transformation optics whichis a versatile mathematical tool enabling a deeper analytical insight into thescattering properties of EM fields in metamaterials. In this subsection, a cloakconsisting of a bianisotropic medium is investigated, and the geometric trans-formation proposed by J. Pendry [2] is proved to be applicable to design suchkind of cloak.We start with a free space with permittivity ε = ε , permeability µ = µ andmagnetoelectric coupling parameters ξ = ξ , ζ = ζ in the Cartesian coordinatesystem, and first map it onto polar coordinates ( r, θ, z ) as defined by x = r cos θ, y = r sin θ, z = z. (44)a disk with R results as shown in Fig. 2(a). Let us then introduce a geometrictransformation defined by (45) which maps the field within this disk onto anannulus R < r ≤ R , i.e. the original point (red) r = 0 in the original disk(Fig. 2(a)) is blowup to a disk with r = R in coordinates ( r (cid:48) , θ (cid:48) , z (cid:48) ) as shownin panel (b) of Fig. 2. Based on this coordinate transformation, the schematic12iagram of the ray trajectories in the cloak is shown in panel (c): the rays aredistorted and guided around the cloak. R R R T.O(a) (b) (c)
Figure 2: Geometric transformation for invisibility cloak: (a) Disk (grey) with radius r = R and the original point (red) is at r = 0 ; (b) The original point is mapped to a disk with r (cid:48) = R , i.e the disk with r = R is compressed to an annulus R < r (cid:48) ≤ R , the mappingfunction is defined by (45); (c) Ray trajectories in the cloak. The mapping function is defined as [2] r (cid:48) = R + r ( R − R ) /R , ≤ r ≤ R θ (cid:48) = θ, < θ ≤ πz (cid:48) = z . (45)which provides the Jacobian matrix J rr (cid:48) = ∂ ( r, θ, z ) ∂ ( r (cid:48) , θ (cid:48) , z (cid:48) ) = diag( R R − R , ,
1) 0 ≤ r ≤ R . (46)Finally, we go back to the Cartesian coordinates ( x (cid:48) , y (cid:48) , z (cid:48) ) , which are radiallycontracted, and the compound Jacobian matrix is J xx (cid:48) = J xr J rr (cid:48) J r (cid:48) x (cid:48) (47)On the other hand, according to the coordinates transformation (44), we have J xr = ∂ ( x, y, z ) ∂ ( r, θ, z ) = R ( θ ) diag(1 , r, J r (cid:48) x (cid:48) = ∂ ( r (cid:48) , θ (cid:48) , z (cid:48) ) ∂ ( x, y, z ) = diag(1 , r (cid:48) , R − ( θ (cid:48) ) (48)with R ( θ ) = cos( θ ) − sin( θ ) 0sin( θ ) cos( θ ) 00 0 1 (49)13hen, (47) becomes J xx (cid:48) = J xr J rr (cid:48) J r (cid:48) x (cid:48) = R ( θ ) diag(1 , r,
1) diag (cid:18) R R − R , , (cid:19) diag(1 , r (cid:48) , R − ( θ (cid:48) )= R ( θ ) diag (cid:18) R R − R , rr (cid:48) , (cid:19) R ( − θ (cid:48) ) . (50)The transformation matrix is correspondingly T − xx (cid:48) = (cid:2) J T xx (cid:48) J xx (cid:48) / det( J xx (cid:48) ) (cid:3) − = R ( θ (cid:48) ) diag (cid:18) r (cid:48) − R r (cid:48) , r (cid:48) r (cid:48) − R , R ( R − R ) r (cid:48) − R r (cid:48) (cid:19) R ( − θ (cid:48) ) . (51)More explicitly, we apply the matrix R ( θ (cid:48) ) into (51) which we denote by T − xx (cid:48) = T T T T
00 0 T , (52)with T = 1 − R r (cid:48) cos ( θ (cid:48) ) + R r (cid:48) − R sin ( θ (cid:48) ) ,T = T = R ( R − r (cid:48) ) r (cid:48) ( r (cid:48) − R ) sin( θ (cid:48) ) cos( θ (cid:48) ) ,T = 1 − R r (cid:48) sin ( θ (cid:48) ) + R r (cid:48) − R cos ( θ (cid:48) ) ,T = R ( R − R ) r (cid:48) − R r (cid:48) . (53)Hence, the parameters in the annulus R < r (cid:48) ≤ R are v = v T − xx (cid:48) , v = ε, µ, ξ, ζ . (54)while the parameters of the outer region r (cid:48) > R are unchanged, and thatof the disk r (cid:48) ≤ R can have any value without affecting the electromagneticscattering.Fig. 3(a) shows our designed bianisotropic cloak, the cloaking region is theinnermost circle of radius R = 0 . um wherein a L-shaped obstacle locates, thegrey annulus represents the distorted space with radius R < r (cid:48) ≤ R , where R = 0 . um. Firstly, the EM distribution of a point source ( s -polarized with14 point source z E R R Figure 3: (a) Schematic diagram of bianisotropic cloak with a L-shaped obstacle inside thecloaking region, and a point source ( s -polarized with the electric field along z with f =8 . × Hz) locates outside the shell. Plots of Re ( E z ) : (b) a point source in a pure matrixwith parameters ε = ε I , µ = µ I and ξ = ζ = 0 . /c I ; (c) same point source in thematrix with a presence of L-shaped obstacle, the parameters of obstacle are ε = (1+5 × i ) ε I , µ = µ I and ξ = ζ = 0 . /c I ; (d) same point source that radiates in the presence of anL-shaped obstacle surrounded by an invisibility cloak. the electric field along z ) with frequency f = 8 . × Hz radiating in a matrixwith parameters v = v I , ( v = ε, µ, ξ, ζ ) is shown in Fig. 3(b). v = ε , µ arethe permittivity, permeability of the vacuum, respectively; while | ξ ζ | (cid:54) = 1 /c with c the velocity of light in vacuum should be promised to ensure convergenceof the numerical algorithm as indicated in (35). If there is a L-shaped obstaclewith ε = (1 + 5 × i ) ε I , µ = µ I , ξ = ζ = 0 . /c I in the matrix, it leadsto an interacting with the source, the plot of Re ( E z ) is indicated in panel (c).15owever, when the L-shaped obstacle is surrounded by the designed cloak, thenthe obstacle seems to be invisible for the outer observer, as shown in panel (d). A cylindrical concentrator is designed to focus the incident waves with wavevectors perpendicular to the cylinder axis, enhancing the electromagnetic energydensity in a given area [4].We start with a map from Cartesian coordinates to cylindrical ones, a cylin-drical lens with r = R is shown in Fig. 4(a), then a radial transformation isintroduced to distort the space as r (cid:48) = R R r ≤ r ≤ R R − R R − R r − R − R R − R R R < r ≤ R θ (cid:48) = θ, ≤ θ < πz (cid:48) = z . (55) R R R (a) (b) Figure 4: (a) Schematic diagram of a concentrator: a space is compressed into a circle withradius R (the innermost in red) with the expense of an expansion of space (grey color)between R and R ; the dotted circle corresponds to a virtual interface used in the geometrictransformation, see (55). (b) Ray trajectories in the concentrator. As shown in Fig. 4(a), a space is compressed into a cylindrical region withradius R (the innermost circle) with the expense of an expansion of space(greycolor) between R and R , where an intermediate circle is located at R . Fig.16(b) shows the ray trajectories in such a concentrator. The associated Jacobianmatrix for the transformation in (55) is J rr (cid:48) = diag( R R , ,
1) 0 ≤ r ≤ R , diag( R − R R − R , , R < r ≤ R . (56)Finally we go back to Cartesian coordinates. The compound Jacobian matrixis then J xx (cid:48) = J xr J rr (cid:48) J r (cid:48) x (cid:48) = R ( θ (cid:48) )diag( R R , rr (cid:48) , R ( − θ (cid:48) ) 0 ≤ r ≤ R R ( θ (cid:48) )diag( R − R R − R , rr (cid:48) , R ( − θ (cid:48) ) R < r ≤ R (57)Furthermore, the transformation matrix is T − xx (cid:48) = J − xx (cid:48) J − T xx (cid:48) det( J xx (cid:48) )= R ( θ (cid:48) )diag(1 , , R R ) R ( − θ (cid:48) ) 0 ≤ r ≤ R R ( θ (cid:48) )diag( R − R R − R rr (cid:48) , R − R R − R r (cid:48) r , R − R R − R rr (cid:48) ) R ( − θ (cid:48) ) R < r ≤ R (58)Let us expand this formula in the distorted cylindrical coordinates ( r (cid:48) , θ (cid:48) , z (cid:48) ) ,and denote T − xx (cid:48) as in (52), finally we have ≤ r (cid:48) ≤ R : T = T = 1 , T = T = 0 , T = R R .R < r (cid:48) ≤ R : T = A cos ( θ (cid:48) ) + B sin ( θ (cid:48) ) ,T = T = ( A − B ) sin( θ (cid:48) ) cos( θ (cid:48) ) ,T = A sin ( θ (cid:48) ) + B cos ( θ (cid:48) ) , T = C . (59)with A = 1 + R − R R − R R r (cid:48) , B = r (cid:48) r (cid:48) + R − R R − R R ,C = ( R − R R − R ) (1 + R − R R − R R r (cid:48) ) . (60)Apply matrix T − xx (cid:48) to (54), we can obtain the formulae of these parametersin the annulus as well as the innermost circle, while the outside space of theconcentrator is bianisotropic matrix with v = v I .17ig. 5(a) is the schematic diagram of our concentrator, the radii of theinner and outer circles are R = 0 . um, R = 0 . um, while the intermediatecircle has a radius R = 0 . um. The parameters of each region are defined in(54) along with transfer matrix T − xx (cid:48) in (59), here the magnetoelectric couplingparameters ξ = ζ = 0 . /c are chosen as the same values as those of theinvisibility cloak. An s -polarized (the electric field is perpendicular to the x - y plane) plane wave with frequency . × Hz is incident from above, the plotof Re ( E z ) is depicted in panel (b): the waves are compressed in the inner circleas expected. (a) (b)-1 1 (um) R R R Figure 5: (a) Schematic diagram of a bianisotropic concentrator. (b) Plot of Re ( E z ) under an s -polarized incidence with frequency f = 8 . × Hz, the parameters in the region r (cid:48) ≤ R ,and the annulus R < r (cid:48) ≤ R are given by (54) with transformation matrix T − xx (cid:48) in (59). Fig. 6 shows a rotator allowing a rotation of the incident wave, the distri-bution of the field inside a disk r ≤ R appears as it comes from a differentangle comparing with that of outside r > R , for example, a π/ angle rotationis achieved here in the inner disk.Be different from the geometric transformation introduced in the concentra-tor, an angular transformation is introduced in this case instead of the radialtransformation. The mapping is performed in the region r < R : the domain in r < R is transformed by rotating a fixed angle from the region r (cid:48) = a in virtual18 R R (a) (b) Figure 6: (a) Schematic diagram of the rotator: incident waves are rotated in the annulus R < r ≤ R and reach the disk of radius R with a π/ degree rotation. (b) Ray trajectoriesin rotator. space, while the rotation angle is continually changed from the fixed angle tozero in region R < r < R . In other words, the transformation algorithm canbe expressed as r ≤ R : r (cid:48) = r, θ (cid:48) = θ + θ , z (cid:48) = z .R < r ≤ R : r (cid:48) = r, θ (cid:48) = θ + θ f ( R ) − f ( r ) f ( R ) − f ( R ) , z (cid:48) = z .r > R : r (cid:48) = r, θ (cid:48) = θ, z (cid:48) = z . (61)where θ is the rotation angle for the inner disk, and it is reduced to zero as theradius approaches to r = R . f ( r ) is an arbitrary continuous function of r .Based on this transformation, we can derive the associated Jacobian matrixas J rr (cid:48) = diag(1 , , r ≤ R and r > R θ f (cid:48) ( r ) f ( R ) − f ( R ) 1 00 0 1 R < r ≤ R (62)again, the compound Jacobian matrix can be obtained by J xx (cid:48) = J xr J rr (cid:48) J r (cid:48) x (cid:48) ,as well as the transformation matrix T − xx (cid:48) = J − xx (cid:48) J − T xx (cid:48) det( J xx (cid:48) ) . In annulus19 < r (cid:48) (= r ) ≤ R , the notations in (52) become T = 1 + r (cid:48) D sin ( θ (cid:48) ) + 2 r (cid:48) D sin( θ (cid:48) ) cos( θ (cid:48) ) ,T = T = r (cid:48) D [sin ( θ (cid:48) ) − cos ( θ (cid:48) )] − r (cid:48) D sin( θ (cid:48) ) cos( θ (cid:48) ) ,T = 1 + r (cid:48) D cos ( θ (cid:48) ) − r (cid:48) D sin( θ (cid:48) ) cos( θ (cid:48) ) , T = 1 . (63)with D = θ f (cid:48) ( r ) f ( R ) − f ( R ) . In regions r (cid:48) ≤ R and r (cid:48) ≥ R , T − xx (cid:48) = I . Again,substituting the matrix T − xx (cid:48) into (54), the expressions for those parameters canbe derived in different regions. (a) (b)-1 1 (um) R R Figure 7: (a) Schematic diagram of the bianisotropic rotator. (b) Plot of Re ( E z ) under an s -polarized incidence with frequency f = 8 . × Hz, the parameters in the region r (cid:48) ≤ R ,annulus R < r (cid:48) ≤ R are given by (54) with transformation matrix T − xx (cid:48) in (63). Similarly, the FEM analysis for the rotator is shown in Fig. 7. The radiiof the inner and outer circles are R = 0 . um, R = 0 . um, respectively. Theparameters of each region are defined in (54) along with transfer matrix T − xx (cid:48) in(63). An s -polarized incident wave with frequency f = 8 . × Hz is assumedto radiate from above. A rotation of θ = π/ of the fields in the inner circlecan be observed in panel (b) comparing with the incidence. Note that, θ canbe an arbitrary values in [ − π, π ] . 20 . A multilayered bianisotropic cloak through homogenization Let us now consider a layered medium periodic along x and invariant along x and x . A periodic cell consists of a heterogeneous bi-isotropic mediumdescribed by piecewise constant scalar valued functions ε ( x ) , µ ( x ) , ξ ( x ) , ζ ( x ) . Using two-scale expansion techniques applied to the Maxwell-Tellegen’sequations like in [20] one gets the homogenized equations ∇ × E eff = ωζ eff E eff + iωµ eff H eff , ∇ × H eff = − iωε eff E eff + ωξ eff H eff . (64)where ε eff , µ eff , ξ eff and ζ eff are anisotropic homogenized tensors such that ε eff = < ε − > − < ε >
00 0 < ε > ,µ eff = < µ − > − < µ >
00 0 < µ > ,ξ eff = < ξ − > − < ξ >
00 0 < ξ > ,ζ eff = < ζ − > − < ζ >
00 0 < ζ > . (65)where < f ( x ) > = (cid:82) d f ( x ) dx , with d the periodic cell size. Note also thatthese formulae can be deduced from the high-order homogenization approach in[21], by keeping only leading order terms.One can use these formulae to design a concentric multilayered bianisotropiccloak of inner radius R = 140 nm and outer radius R = 300 nm, consistingof 22 homogeneous layers (see table 1 for the values of the inverse relative21ermittivity ε r , permeability µ r , magnetoelectric coupling parameters ξ r and ζ r ) in a bianisotropic background with permittivity ε = ε , permeability µ = µ and magnetoelectric coupling parameters ξ = ξ = 0 . /c , ζ = ζ = 0 . /c . Table 1: Parameters of the layered cloak. ( ε = ε r ε , µ = µ r µ , ξ = ξ r ξ , ζ = ζ r ζ ) layer 1 (inner) 2 3 4 5 6thickness [nm] 10.0 3.75 7.5 7.5 7.5 7.5 ε − r = µ − r = ζ − r = ξ − r ε − r = µ − r = ζ − r = ξ − r ε − r = µ − r = ζ − r = ξ − r ε − r = µ − r = ζ − r = ξ − r nm) with aconcentric multilayered bianisotropic cloak is much reduced in Fig. 8(b) com-paring with Fig. 8(a), and the shadow zone behind the obstacle is also muchreduced with the layered cloak. The interval of cloaking frequencies from to THz is fairly broadband. 22 a) (b)(c) (d)
Figure 8: (a) Scattering of an infinite conducting obstacle (of radius nm) by a point sourceat a frequency
THz; (b) Scattering like in (a) when the obstacle is surrounded by a layeredcloak. (c) Scattering like in (a) for a point source at a frequency
THz; (d) Scattering likein (b) for a point source at a frequency
THz. Color scale corresponds to normalized valuesof Re ( E z ) .
7. Conclusion
In this paper, we have discussed the PDEs model for the analysis of the prop-agation properties in bianisotropic media with an invariance along one direction:two coupled PDEs with scalar solutions are derived for this electromagnetic sys-tem. The spatially varying coefficients within these two PDEs can be scalar ormatrix valued. Regarding the numerical model, specially designed perfectlymatched layers (PMLs) are introduced in order to account for the unbounded23omain. To do this, we consider a transformation which maps an infinite domainonto a finite one, and we further add the damping. This provides a reflection-less bianisotropic layer for all incident angles, and this can be considered as anextension of the PMLs introduced by J. Berenger nearly twenty years ago [14].We implement these PMLs in scattering problems in bianisotropic media in-cluding an invisibility cloak, a field concentrator and a field rotator designed byintroducing the proper geometric transformation. Importantly, these functionalstructures are proved to work well when they are achieved by bianisotropic me-dia. We finally, show that one can achieve near-cloaking with a layered cloakwith piecewise constant permittivity, permeability and magnetoelectric couplingparameters. The solutions of singular systems considered in this paper representimportant benchmarks for the validation of non-trivial numerical calculationsfor bianisotropic media.
8. Acknowledgments
S.Guenneau acknowledges a European funding through ERC Starting GrantNo. 279673.