# Extended Linear Embedding via Green's Operators for Analyzing Wave Scattering from Anisotropic Bodies

Authors: V. Lancellotti; A. G. Tijhuis
Journal: International Journal of Antennas and Propagation
 (2014)

---

## Abstract

Linear embedding via Green’s operators (LEGO) is a domain decomposition method particularly well suited for the solution of scattering and radiation problems comprised of many objects. The latter are enclosed in simple-shaped subdomains (electromagnetic bricks) which are in turn described by means of scattering operators. In this paper we outline the extension of the LEGO approach to the case of penetrable objects with dyadic permittivity or permeability. Since a volume integral equation is only required to solve the scattering problem inside a brick and the scattering operators are inherently surface operators, the LEGO procedureper se can afford a reduction of the number of unknowns in the numerical solution with the Method of Moments and subsectional basis functions. Further substantial reduction is achieved with the eigencurrents expansion method (EEM) which employs the eigenvectors of the scattering operator as local entire-domain basis functions over a brick’s surface. Through a few selected numerical examples we discuss the validation and the efficiency of the LEGO-EEM technique applied to clusters of anisotropic bodies.


---

## Body

## 1. Introduction

It is straightforward to incorporate anisotropic materials into macroscopic Maxwell’s equations by allowing for constitutive relationships with dyadic permittivity or permeability [1]. Nevertheless solving an electromagnetic (EM) scattering or radiation problem which involves objects with dyadic constitutive parameters demands the application of some numerical method.For instance, the direct solution of Maxwell’s equations in differential form with the finite-element method (FEM) [2] is tempting, because the latter procedure leads to highly sparse matrices. On the other hand, radiation boundary conditions are naturally taken into account by means of volume integral equations (VIEs) [3–6], but their numerical treatment with the Method of Moments (MoM) produces dense matrices (e.g., [7, 8]). In either case, when the structure of concern is electrically large, one could apply an iterative technique to solve the resulting algebraic linear system. However, the convergence may be slow, unless the matrix is well conditioned; poor conditioning seems to especially plague problems that involve objects with large contrast with respect to the background medium [9].Solution strategies based on domain decomposition methods (DDMs) may help overcome the foresaid drawbacks, because rather than tackling a given EM problem as a whole, DDMs handle the local strong interactions and the distant weaker ones in a large composite structure separately. This is accomplished by enclosing “small” parts of a large structure into subdomains that are treated independently of one another. The multiple scattering that occurs among the various subdomains is then described by means of translation or transfer operators. In this way, the usage of VIEs is restricted to the solution of the EM problem within the subdomains and hence, the resulting “smaller” matrices are likely less prone to ill conditioning. Finally, a surface integral equation (SIE) for the unknown fields or equivalent currents on the subdomains’ boundaries is stated and solved numerically. Among the methods that follow this paradigm we mention the equivalence principle algorithm (EPA) [10] and the tangential equivalence principle algorithm (T-EPA) [11], the generalized surface integral equation (GSIE) method [12], and the linear embedding via Green’s operator (LEGO) technique [13, 14].The domain decomposition in LEGO is effected by enclosing the objects into 3D simple-shaped domains (Figure1) which constitute the building EM “bricks” of the model. Firstly, each brick is described through a scattering operator of equivalent surface current densities defined over a brick’s surface. Secondly, the multiple scattering occurring between any two bricks is captured by means of transfer operators. Since the mathematical framework for the separation is provided by the Surface Equivalence Principle (SEP) [1, 15], LEGO yields a system of as many coupled SIEs for the unknown surface current densities as the number of bricks in the model [14, 16].Figure 1
A 3D EM brick𝒟
k immersed in a background medium (①) encloses an anisotropic object (③) which is surrounded by a host medium (②); the EM behavior of 𝒟
k is captured by the scattering operator S
k
k.Now, handling anisotropic media with LEGO is straightforward, in that there is no limitation as to the constitutive parameters of the objects inside a brick. Needless to say, the calculation of the relevant scattering operator requires the solution of the EM scattering problem inside the brick by means of a VIE, as anticipated. And yet, once the scattering operators of all the bricks in the model have been computed, the usual equations of LEGO can be stated and solved as before [14, 16, 17].The important point is that a VIE has to be solved only for one object at a time, if many different objects form the cluster, or just once, in case the objects are all equal to each other. Furthermore, in the context of the numerical solution with the MoM, the (surface) basis functions on a brick’s boundary are usually far fewer than the (volume) basis functions within the object. Therefore, it is expected that the LEGO procedure can reduce the computational burden as compared to the workload that would be required to solve a single VIE for all the objects simultaneously. Additionally, substantial reduction of the degrees of freedom is achieved thanks to the eigencurrents expansion method (EEM) in which the eigenvectors of the scattering operator are employed as local entire-domain basis functions over a brick’s surface [14, 18, 19].The rest of this paper is organized as follows. In Section2 we outline the derivation of the scattering operator of a brick that embeds an anisotropic object, and we give detailed expressions for the special case of dyadic permittivity and scalar permeability. Besides, we review the LEGO equations involving the total inverse scattering operator of a composite structure. In Section 3 we describe the numerical solution of the LEGO equations with the MoM and the order reduction through the EEM. Lastly, in Section 4 we present an example of validation, a study of the convergence properties of the EEM, and a numerical example that involves a complex multiscale structure comprised of a magnetized plasma. Efficiency in terms of matrix sizes and computational times is briefly addressed as well.A time dependence for fields and sources in the formexp
⁡
(
j
ω
t
) is assumed and suppressed throughout.


## 2. Formulation with LEGO

### 2.1. Generalities

Simply put, the solution of a scattering problem which involves an aggregate ofN
D anisotropic bodies consists of three main steps [14, 16, 17]: (i) definition of the N
D EM bricks 𝒟
k, k
=
1
,
…
,
N
D and calculation of the scattering operators S
k
k thereof (Figure 1); (ii) calculation of the N
T
≤
N
T
,
max
⁡
=
N
D
×
(
N
D
-
1
) transfer operators T
n
k, T
k
n, n
≠
k, and n
=
1
,
…
,
N
D, which enter the total inverse scattering operator S
-
1 of the structure; and (iii) numerical inversion of resulting surface integral equation.In this section, we elaborate on steps (i) and (ii), whereas we will deal with step (iii) in Section3. At this stage, we make no hypothesis with regard to the position, the shape, and the constitutive parameters of the N
D bodies. Nevertheless, in the formulation to be developed hereinafter we consider N
D bricks of identical shape. The latter assumption—which should not constitute an issue especially for arbitrary displacements of objects and bricks—greatly facilitates the numerical solution, as will be clear in Section 4.


### 2.2.  EM Bricks with Anisotropic Media

Sketched in Figure1 is a quite general example of EM brick immersed in an infinite background homogeneous medium (labelled ①). The brick comprises a simple-shaped domain filled with an isotropic homogeneous host medium (tagged ②) which in turn wraps the anisotropic body (tagged ③). For the sake of completeness, we point out that in [14] we considered the case of PEC or penetrable isotropic objects and media ① and ② with the same constitutive parameters, whereas in [16] we introduced the notion of scattering operator S
¯
J
k of a material interface (i.e., a brick’s boundary, ∂
𝒟
k) which allows handling the instance of different media ① and ② efficiently.To keep the exposition lucid, with the aid of Figure2 we first outline the general procedure for deriving the scattering operator in the case of background and host medium having the same EM properties. Then, we specialize the equations and give the explicit form of the operators involved for the instance of medium ③ with dyadic permittivity ɛ
¯
3 and scalar permeability μ
3
=
μ
2
=
μ
1. (The case of medium ③ with dyadic permeability μ
¯
3 and scalar permittivity ɛ
3
=
ɛ
2
=
ɛ
1 is similar and ensues in a straightforward manner.) Moreover, the scattering operator of a material interface along with the cascading procedure described in [16] can be employed to obtain S
k
k for the general situation of Figure 1, though we will not examine this configuration in what follows.Derivation of the scattering operator of𝒟
k in the case when host (②) and background (①) medium have the same EM properties. (a) The equivalent incident currents q
k
i on ∂
𝒟
k produce incident fields F
o
i in the region occupied by the object; (b) the object “reacts” with polarization and/or magnetization currents q
o; (c) the currents q
o produce scattered fields F
t
k
s on ∂
𝒟
k; (d) the scattered fields are replaced with equivalent scattered currents q
k
s.
(a)
(b)
(c)
(d)With these preliminaries we are now ready to elaborate on the derivation ofS
k
k.The starting point is a solitary brick𝒟
k illuminated by incident fields E
k
i, H
k
i, radiated by someexternal sources. An application of the SEP [1] allows replacing the effect of the external sources with equivalent electric (J
k
i
=
n
^
k
×
H
k
i) and magnetic (M
k
i
=
E
k
i
×
n
^
k) surface current densities flowing on ∂
𝒟
k [14]. The latter are collected in the abstract 2
×
1 vector q
k
i for ease of manipulation; namely,
(1)
q
k
i
=
[
η
1
J
k
i
-
M
k
i
η
1
]
,
with η
1
=
(
μ
1
/
ɛ
1
)
1
/
2 being the intrinsic impedance in medium ①. The superscript “i” serves to remind us that the currents in question reproduce the incident fields towards the inside of the brick (see Figure 2(a)). We notice in passing that J
k
i and M
k
i can be thought of as being positioned on either side of ∂
𝒟
k, because the brick’s boundary is not a material interface by hypothesis (cf. [16]).To proceed we observe that the object inside𝒟
k is illuminated by the EM field generated by q
k
i. As suggested in Figure 2(a), such field can be related to its sources symbolically by means of a linear operator; that is,
(2)
F
o
i
=
P
o
k
q
k
i
,
where F
o
i is the abstract 2
×
1 vector
(3)
F
o
i
=
[
E
o
i
η
1
η
1
H
o
i
]
,
and thepropagator
P
o
k is a 2
×
2 matrix of dyadicsurface integral operators whose kernel is the dyadic Green’s function of media ① and ②.Now, by invoking the Volume Equivalence Principle (VEP) [1] we replace the object with volume current densities [20]. In the general case of both electrically and magnetically penetrable objects, said currents are related to the total flux densities D
o and B
o (e.g., [7]). The latter can be obtained by solving a set of VIEs; namely,(4)
X
o
o
q
o
=
F
o
i
,
where
(5)
q
o
=
[
D
o
(
ɛ
1
η
1
)
η
1
B
o
μ
1
]
,
and X
o
o is a 2
×
2 matrix of volume integral operators over the region 𝒱
o occupied by medium ③.Next, we observe that the volume “currents”q
o generate scattered fields on the boundary ∂
𝒟
k, as shown in Figure 2(c); in symbols, this reads
(6)
F
t
k
s
=
P
k
o
q
o
,
where F
t
k
s is the abstract 2
×
1 vector
(7)
F
t
k
s
=
[
0
η
1
H
t
k
s
]
,
and the propagator P
k
o is a 2
×
2 abstract matrix of dyadicvolume integral operators involving dyadic Green’s function of media ① and ②; besides, the subscript “t” stands for “tangential.” The superscript “s” (short for “scattered”) signifies that the fields of concern here propagate away from the object and outwards the brick.Lastly, as is graphically explained in Figure2(d), a second application of the SEP over ∂
𝒟
k allows substituting the scattered fields F
t
k
s with equivalent surface electric (J
k
s) and magnetic (M
k
s) current densities [14]. In symbols,
(8)
P
k
k
s
q
k
s
=
F
t
k
s
,
with q
k
s being an abstract 2
×
1 vector defined in a similar fashion to (1). The propagator P
k
k
s is exactly the same as in [14, Table 1], since it depends only on the shape of ∂
𝒟
k and not on the media inside 𝒟
k.By formally solving (2), (4), (6), and (8) for q
k
s as a function of q
k
i we find (this expression differs from [14, Equation (
11
)] for a minus sign, which is due to different definitions adopted for the operator X
o
o)
(9)
q
k
s
=
(
P
k
k
s
)
-
1
P
k
o
(
X
o
o
)
-
1
P
o
k
q
k
i
=
S
k
k
q
k
i
,
which provides us with a suitable recipe for computing the scattering operator S
k
k.In the notable instance of an anisotropic dielectric object withμ
3
=
μ
2
=
μ
1, the operators P
o
k, X
o
o, and P
k
o reduce to 1
×
2, 1
×
1, and 2
×
1 abstract matrices whose explicit expressions read
(10)
P
o
k
{
∘
}
=
[
-
j
k
1
∫
∂
𝒟
k
‍
d
2
r
′
G
¯
1
(
R
)
·
{
∘
}
-
∫
∂
𝒟
k
‍
d
2
r
′
∇
′
G
1
(
R
)
×
I
¯
·
{
∘
}
]
,
r
∈
𝒱
o
,
(11)
X
o
o
{
∘
}
=
ɛ
1
ɛ
¯
3
-
1
(
r
)
·
{
∘
}
-
k
1
2
∫
𝒱
o
‍
d
3
r
′
G
¯
1
(
R
)
·
α
¯
3
(
r
′
)
·
{
∘
}
,
r
∈
𝒱
o
,
(12)
P
k
o
{
∘
}
=
[
0
¯
·
{
∘
}
-
j
k
1
I
¯
s
·
∫
𝒱
o
‍
d
3
r
′
∇
′
G
1
(
R
)
×
α
¯
3
(
r
′
)
·
{
∘
}
]
,
r
∈
∂
𝒟
k
,
with k
1 being the wavenumber in medium ①, R
=
|
r
-
r
′
| being the distance between source (r) and observation (r
′) points, and(13)
G
¯
1
(
R
)
=
I
¯
G
1
(
R
)
-
∇
∇
′
G
1
(
R
)
k
1
2
,
G
1
(
R
)
=
exp
⁡
(
-
j
k
1
R
)
4
π
R
,
α
¯
3
(
r
′
)
=
I
¯
-
[
ɛ
¯
3
(
r
′
)
ɛ
1
]
-
1
,
and I
¯
s
=
I
¯
-
n
^
k
n
^
k denotes the surface unit dyadic on ∂
𝒟
k. The null vector in (7) and the null dyadic appearing in P
k
o in (12) are a consequence of the definition of the propagator P
k
k
s which relies only on the magnetic field H
t
k
s to yield the corresponding equivalent surface current densities on ∂
𝒟
k (cf. [14], Table 1).We conclude by observing that the propagatorsP
o
k, P
k
o and the operator X
o
o are dimensionless, as can be ascertained by inspection of (10), (12), and (11). In contrast, current and field vectors q
k
s
,
i, q
o, F
t
k
s, and F
o
i are all normalized so as to carry the physical dimension of W
1
/
2
/
m.


### 2.3.  Total Inverse Scattering Operator and LEGO Equations

When two or more bricks are considered, (9) must be modified to account for the multiple scattering that occurs between the bricks. This is readily accomplished upon recognizing that the equivalent scattered currents q
n
s on ∂
𝒟
n produce fields that reach the boundary of 𝒟
k alongside the incident fields due to the independent sources of the EM problem. By invoking the SEP, the extra incident fields can be turned into equivalent surface current densities q
k
(
n
)
i on ∂
𝒟
k. In symbols,
(14)
P
k
k
i
q
k
(
n
)
i
=
P
k
n
q
n
s
,
n
≠
k
,
where the propagators P
k
k
i, P
k
n are the same as in [14, Table 1], as they depend only on the shape of ∂
𝒟
k
,
n. The formal inversion of (14) yields
(15)
q
k
(
n
)
i
=
(
P
k
k
i
)
-
1
P
k
n
q
n
s
=
T
k
n
q
n
s
,
n
≠
k
,
where the dimensionlesstransfer operator T
k
n maps scattered currents on ∂
𝒟
n onto incident currents on ∂
𝒟
k [14].In the light of the previous discussion and thanks to (15), the extension of (9) to the case of N
D interacting bricks reads [14]
(16)
q
k
s
=
S
k
k
(
q
k
i
+
∑
n
≠
k
n
=
1
N
D
‍
T
k
n
q
n
s
)
,
k
=
1
,
…
,
N
D
,
which represents a system of N
D coupled surface integral equations. Upon grouping all the scattered and incident currents q
k
s
,
i column wise into q
s and q
i, (16) can be expressed in a compact form as
(17)
S
-
1
q
s
=
q
i
,
with
(18)
q
s
=
[
q
1
s
⋮
q
N
D
s
]
,
q
i
=
[
q
1
i
⋮
q
N
D
i
]
,
(19)
S
-
1
=
[
S
11
-
1
-
T
12
⋯
-
T
1
N
D
-
T
21
S
22
-
1
⋯
-
T
2
N
D
⋮
⋮
⋱
⋮
-
T
N
D
1
-
T
N
D
2
⋯
S
N
D
N
D
-
1
]
where S
-
1 is the total inverse scattering operator of the structure. Equation (17) along with (19) constitutes the formulation of the EM problem.Notice that (19) is a formal result not only because ordinarily S
k
k and its inverse cannot be obtained in closed form, but also because S
k
k
-
1 may not existper se when S
k
k happens to be singular. In spite of these complications, the numerical calculation of q
s through (19) is viable and stable, as long as we carefully handle the possible null or small eigenvalues of S
k
k [14, 18], as we show in Section 3.2.


## 2.1. Generalities

Simply put, the solution of a scattering problem which involves an aggregate ofN
D anisotropic bodies consists of three main steps [14, 16, 17]: (i) definition of the N
D EM bricks 𝒟
k, k
=
1
,
…
,
N
D and calculation of the scattering operators S
k
k thereof (Figure 1); (ii) calculation of the N
T
≤
N
T
,
max
⁡
=
N
D
×
(
N
D
-
1
) transfer operators T
n
k, T
k
n, n
≠
k, and n
=
1
,
…
,
N
D, which enter the total inverse scattering operator S
-
1 of the structure; and (iii) numerical inversion of resulting surface integral equation.In this section, we elaborate on steps (i) and (ii), whereas we will deal with step (iii) in Section3. At this stage, we make no hypothesis with regard to the position, the shape, and the constitutive parameters of the N
D bodies. Nevertheless, in the formulation to be developed hereinafter we consider N
D bricks of identical shape. The latter assumption—which should not constitute an issue especially for arbitrary displacements of objects and bricks—greatly facilitates the numerical solution, as will be clear in Section 4.


## 2.2.  EM Bricks with Anisotropic Media

Sketched in Figure1 is a quite general example of EM brick immersed in an infinite background homogeneous medium (labelled ①). The brick comprises a simple-shaped domain filled with an isotropic homogeneous host medium (tagged ②) which in turn wraps the anisotropic body (tagged ③). For the sake of completeness, we point out that in [14] we considered the case of PEC or penetrable isotropic objects and media ① and ② with the same constitutive parameters, whereas in [16] we introduced the notion of scattering operator S
¯
J
k of a material interface (i.e., a brick’s boundary, ∂
𝒟
k) which allows handling the instance of different media ① and ② efficiently.To keep the exposition lucid, with the aid of Figure2 we first outline the general procedure for deriving the scattering operator in the case of background and host medium having the same EM properties. Then, we specialize the equations and give the explicit form of the operators involved for the instance of medium ③ with dyadic permittivity ɛ
¯
3 and scalar permeability μ
3
=
μ
2
=
μ
1. (The case of medium ③ with dyadic permeability μ
¯
3 and scalar permittivity ɛ
3
=
ɛ
2
=
ɛ
1 is similar and ensues in a straightforward manner.) Moreover, the scattering operator of a material interface along with the cascading procedure described in [16] can be employed to obtain S
k
k for the general situation of Figure 1, though we will not examine this configuration in what follows.Derivation of the scattering operator of𝒟
k in the case when host (②) and background (①) medium have the same EM properties. (a) The equivalent incident currents q
k
i on ∂
𝒟
k produce incident fields F
o
i in the region occupied by the object; (b) the object “reacts” with polarization and/or magnetization currents q
o; (c) the currents q
o produce scattered fields F
t
k
s on ∂
𝒟
k; (d) the scattered fields are replaced with equivalent scattered currents q
k
s.
(a)
(b)
(c)
(d)With these preliminaries we are now ready to elaborate on the derivation ofS
k
k.The starting point is a solitary brick𝒟
k illuminated by incident fields E
k
i, H
k
i, radiated by someexternal sources. An application of the SEP [1] allows replacing the effect of the external sources with equivalent electric (J
k
i
=
n
^
k
×
H
k
i) and magnetic (M
k
i
=
E
k
i
×
n
^
k) surface current densities flowing on ∂
𝒟
k [14]. The latter are collected in the abstract 2
×
1 vector q
k
i for ease of manipulation; namely,
(1)
q
k
i
=
[
η
1
J
k
i
-
M
k
i
η
1
]
,
with η
1
=
(
μ
1
/
ɛ
1
)
1
/
2 being the intrinsic impedance in medium ①. The superscript “i” serves to remind us that the currents in question reproduce the incident fields towards the inside of the brick (see Figure 2(a)). We notice in passing that J
k
i and M
k
i can be thought of as being positioned on either side of ∂
𝒟
k, because the brick’s boundary is not a material interface by hypothesis (cf. [16]).To proceed we observe that the object inside𝒟
k is illuminated by the EM field generated by q
k
i. As suggested in Figure 2(a), such field can be related to its sources symbolically by means of a linear operator; that is,
(2)
F
o
i
=
P
o
k
q
k
i
,
where F
o
i is the abstract 2
×
1 vector
(3)
F
o
i
=
[
E
o
i
η
1
η
1
H
o
i
]
,
and thepropagator
P
o
k is a 2
×
2 matrix of dyadicsurface integral operators whose kernel is the dyadic Green’s function of media ① and ②.Now, by invoking the Volume Equivalence Principle (VEP) [1] we replace the object with volume current densities [20]. In the general case of both electrically and magnetically penetrable objects, said currents are related to the total flux densities D
o and B
o (e.g., [7]). The latter can be obtained by solving a set of VIEs; namely,(4)
X
o
o
q
o
=
F
o
i
,
where
(5)
q
o
=
[
D
o
(
ɛ
1
η
1
)
η
1
B
o
μ
1
]
,
and X
o
o is a 2
×
2 matrix of volume integral operators over the region 𝒱
o occupied by medium ③.Next, we observe that the volume “currents”q
o generate scattered fields on the boundary ∂
𝒟
k, as shown in Figure 2(c); in symbols, this reads
(6)
F
t
k
s
=
P
k
o
q
o
,
where F
t
k
s is the abstract 2
×
1 vector
(7)
F
t
k
s
=
[
0
η
1
H
t
k
s
]
,
and the propagator P
k
o is a 2
×
2 abstract matrix of dyadicvolume integral operators involving dyadic Green’s function of media ① and ②; besides, the subscript “t” stands for “tangential.” The superscript “s” (short for “scattered”) signifies that the fields of concern here propagate away from the object and outwards the brick.Lastly, as is graphically explained in Figure2(d), a second application of the SEP over ∂
𝒟
k allows substituting the scattered fields F
t
k
s with equivalent surface electric (J
k
s) and magnetic (M
k
s) current densities [14]. In symbols,
(8)
P
k
k
s
q
k
s
=
F
t
k
s
,
with q
k
s being an abstract 2
×
1 vector defined in a similar fashion to (1). The propagator P
k
k
s is exactly the same as in [14, Table 1], since it depends only on the shape of ∂
𝒟
k and not on the media inside 𝒟
k.By formally solving (2), (4), (6), and (8) for q
k
s as a function of q
k
i we find (this expression differs from [14, Equation (
11
)] for a minus sign, which is due to different definitions adopted for the operator X
o
o)
(9)
q
k
s
=
(
P
k
k
s
)
-
1
P
k
o
(
X
o
o
)
-
1
P
o
k
q
k
i
=
S
k
k
q
k
i
,
which provides us with a suitable recipe for computing the scattering operator S
k
k.In the notable instance of an anisotropic dielectric object withμ
3
=
μ
2
=
μ
1, the operators P
o
k, X
o
o, and P
k
o reduce to 1
×
2, 1
×
1, and 2
×
1 abstract matrices whose explicit expressions read
(10)
P
o
k
{
∘
}
=
[
-
j
k
1
∫
∂
𝒟
k
‍
d
2
r
′
G
¯
1
(
R
)
·
{
∘
}
-
∫
∂
𝒟
k
‍
d
2
r
′
∇
′
G
1
(
R
)
×
I
¯
·
{
∘
}
]
,
r
∈
𝒱
o
,
(11)
X
o
o
{
∘
}
=
ɛ
1
ɛ
¯
3
-
1
(
r
)
·
{
∘
}
-
k
1
2
∫
𝒱
o
‍
d
3
r
′
G
¯
1
(
R
)
·
α
¯
3
(
r
′
)
·
{
∘
}
,
r
∈
𝒱
o
,
(12)
P
k
o
{
∘
}
=
[
0
¯
·
{
∘
}
-
j
k
1
I
¯
s
·
∫
𝒱
o
‍
d
3
r
′
∇
′
G
1
(
R
)
×
α
¯
3
(
r
′
)
·
{
∘
}
]
,
r
∈
∂
𝒟
k
,
with k
1 being the wavenumber in medium ①, R
=
|
r
-
r
′
| being the distance between source (r) and observation (r
′) points, and(13)
G
¯
1
(
R
)
=
I
¯
G
1
(
R
)
-
∇
∇
′
G
1
(
R
)
k
1
2
,
G
1
(
R
)
=
exp
⁡
(
-
j
k
1
R
)
4
π
R
,
α
¯
3
(
r
′
)
=
I
¯
-
[
ɛ
¯
3
(
r
′
)
ɛ
1
]
-
1
,
and I
¯
s
=
I
¯
-
n
^
k
n
^
k denotes the surface unit dyadic on ∂
𝒟
k. The null vector in (7) and the null dyadic appearing in P
k
o in (12) are a consequence of the definition of the propagator P
k
k
s which relies only on the magnetic field H
t
k
s to yield the corresponding equivalent surface current densities on ∂
𝒟
k (cf. [14], Table 1).We conclude by observing that the propagatorsP
o
k, P
k
o and the operator X
o
o are dimensionless, as can be ascertained by inspection of (10), (12), and (11). In contrast, current and field vectors q
k
s
,
i, q
o, F
t
k
s, and F
o
i are all normalized so as to carry the physical dimension of W
1
/
2
/
m.


## 2.3.  Total Inverse Scattering Operator and LEGO Equations

When two or more bricks are considered, (9) must be modified to account for the multiple scattering that occurs between the bricks. This is readily accomplished upon recognizing that the equivalent scattered currents q
n
s on ∂
𝒟
n produce fields that reach the boundary of 𝒟
k alongside the incident fields due to the independent sources of the EM problem. By invoking the SEP, the extra incident fields can be turned into equivalent surface current densities q
k
(
n
)
i on ∂
𝒟
k. In symbols,
(14)
P
k
k
i
q
k
(
n
)
i
=
P
k
n
q
n
s
,
n
≠
k
,
where the propagators P
k
k
i, P
k
n are the same as in [14, Table 1], as they depend only on the shape of ∂
𝒟
k
,
n. The formal inversion of (14) yields
(15)
q
k
(
n
)
i
=
(
P
k
k
i
)
-
1
P
k
n
q
n
s
=
T
k
n
q
n
s
,
n
≠
k
,
where the dimensionlesstransfer operator T
k
n maps scattered currents on ∂
𝒟
n onto incident currents on ∂
𝒟
k [14].In the light of the previous discussion and thanks to (15), the extension of (9) to the case of N
D interacting bricks reads [14]
(16)
q
k
s
=
S
k
k
(
q
k
i
+
∑
n
≠
k
n
=
1
N
D
‍
T
k
n
q
n
s
)
,
k
=
1
,
…
,
N
D
,
which represents a system of N
D coupled surface integral equations. Upon grouping all the scattered and incident currents q
k
s
,
i column wise into q
s and q
i, (16) can be expressed in a compact form as
(17)
S
-
1
q
s
=
q
i
,
with
(18)
q
s
=
[
q
1
s
⋮
q
N
D
s
]
,
q
i
=
[
q
1
i
⋮
q
N
D
i
]
,
(19)
S
-
1
=
[
S
11
-
1
-
T
12
⋯
-
T
1
N
D
-
T
21
S
22
-
1
⋯
-
T
2
N
D
⋮
⋮
⋱
⋮
-
T
N
D
1
-
T
N
D
2
⋯
S
N
D
N
D
-
1
]
where S
-
1 is the total inverse scattering operator of the structure. Equation (17) along with (19) constitutes the formulation of the EM problem.Notice that (19) is a formal result not only because ordinarily S
k
k and its inverse cannot be obtained in closed form, but also because S
k
k
-
1 may not existper se when S
k
k happens to be singular. In spite of these complications, the numerical calculation of q
s through (19) is viable and stable, as long as we carefully handle the possible null or small eigenvalues of S
k
k [14, 18], as we show in Section 3.2.


## 3. Numerical Solution with MoM and EEM

The numerical solution of (17) with the Method of Moments (MoM) in Galerkin’s form and the eigencurrents expansion method (EEM) was extensively presented in [14, 18]. For ease of reference we briefly recall the key points of the procedure herein below.


### 3.1. Reduction to a Weak Form with MoM

We employ the MoM with subsectional basis functions and a symmetric inner product repeatedly in order to (a) turn (2) and (6) into a weak form; and (b) formally invert (4), (8), and (14). In particular, we model the bricks’ boundaries and the volume of the anisotropic objects with 3D triangular-faceted tessellations and tetrahedral meshes, respectively. Then, we associate 2
N
F (N
O) surface (volume) divergence-conforming linear vector functions [21, 22] with the surface (volume) mesh on (inside) each brick in order to expand the surface currents J
k
s
,
i and M
k
s
,
i (the flux density D
o or B
o).As a result, we can write the algebraic counterparts ofS
k
k and T
k
n as
(20)
[
S
k
k
]
=
[
P
k
k
s
]
-
1
[
P
k
o
]
[
X
o
o
]
-
1
[
P
k
o
]
,
(21)
[
T
k
n
]
=
[
P
k
k
i
]
-
1
[
P
k
n
]
,
n
≠
k
,
where [
S
k
k
] and [
T
k
n
] have size 2
N
F
×
2
N
F and the rank of [
X
o
o
] is N
O.Upon using (20) and (21) as building blocks, the algebraic counterpart of (17) can be written formally as
(22)
[
S
]
-
1
[
q
s
]
=
[
q
i
]
,
where the 2
N
F
N
D
×
1 vectors [
q
s
,
i
] store the MoM expansion coefficients of q
k
s
,
i. As anticipated, though, explicitly building and then numerically inverting the linear system (22) are far from being trivial for two reasons.First and foremost, the inverse scattering operators[
S
k
k
]
-
1 may not be defined in some cases. Specifically, if the volume basis functions N
O in the tetrahedral mesh are fewer than the number 2
N
F of surface basis functions on ∂
𝒟
k, then, as is clear from (20), the scattering operator issingular, with its rank being min
⁡
{
2
N
F
,
N
O
}
=
N
O. At any rate, even though 2
N
F
≤
N
O, [
S
k
k
] may still be nearly singular due to the presence of a few very small eigenvalues. On the whole, inverting [
S
k
k
] is not only inconvenient, but also a potentially ill-conditioned problem.Secondly, difficulties posed by the inversion of[
S
k
k
] aside, the size of the total inverse scattering operator, may be quite large, if one considers a problem modelled with a moderate to large number of EM bricks. As a consequence, inversion with direct solvers such as LU factorization [23] may not be feasible. Then, iterative solvers seem an obvious alternative, but the convergence rate is expected to be poor, since [
S
]
-
1 is ill conditioned, in the light of the previous observations.To circumvent the hurdles above we apply the EEM [14].


### 3.2. Compression and Inversion with Eigencurrents

In the EEM we employ the eigenvectors[
V
k
k
] of [
S
k
k
] as local entire-domain basis functions to expand the scattered currents [
q
k
s
] on ∂
𝒟
k. Since these eigenvectors in tandem with the low-level basis functions define surface currents over ∂
𝒟
k, we refer to the columns of [
V
k
k
] as to the eigencurrents of 𝒟
k.The eigencurrents can be separated into two subsets:coupled anduncoupled. The former—associated with the larger eigenvalues of [
S
k
k
]—produce substantial radiation and thus contribute to the multiple scattering that occurs among the N
D bricks. Conversely, the uncoupled eigencurrents—which correspond with the smaller or null eigenvalues of [
S
k
k
]—yield scattered fields that are practically confined around a brick and hence serve to describe the scattered currents on a brick. We denote by N
C
≤
N
C
,
max
⁡
=
min
⁡
{
2
N
F
,
N
O
} the number of coupled eigencurrents.The notion of coupled-uncoupled eigencurrents allows us to partition the matrix[
S
]
-
1 (expressed in the basis of the eigencurrents) into four subblocks [14, Equation (
31
)]; namely,
(23)
[
S
^
]
-
1
=
[
[
S
~
C
C
]
-
1
[
S
~
C
U
]
-
1
[
S
~
U
C
]
-
1
[
S
~
U
U
]
-
1
]
,
where the subscripts “C” (“U”) stand for coupled (uncoupled).In the spectral reordered representation (23) of the total inverse scattering operator the subblock [
S
~
U
U
]
-
1—relevant to the interaction between any two uncoupled eigencurrents from either two distinct bricks or the same brick—can be approximated with its diagonal, that is, the reciprocal of the eigenvalues associated with all the uncoupled eigencurrents, because the latter interact only with themselves. Likewise, the off-diagonal subblocks [
S
~
C
U
]
-
1, [
S
~
U
C
]
-
1 can be zeroed altogether, as they account for the negligible interaction between a pair of coupled and uncoupled eigencurrents. Lastly, the diagonal subblock [
S
~
C
C
]
-
1—relative to the interaction of any two coupled eigencurrents from either two distinct bricks or the same brick—has to be retained in full.In the end, the inverse of[
S
^
]
-
1 is approximately given by
(24)
[
S
^
]
≈
[
[
S
~
C
C
]
[
0
]
[
0
]
[
Λ
U
U
]
]
,
where [
Λ
U
U
] is a diagonal matrix that contains the (
2
N
F
-
N
C
)
×
N
D eigenvalues germane to the uncoupled eigencurrents of the problem. Since [
Λ
U
U
] is known the moment the eigenvalues of [
S
k
k
] have been determined, evidently, the occurrence of small or null eigenvalues does not pose any problem at all in (24). Furthermore, since the number N
C of coupled eigencurrents contributed by 𝒟
k is usually far smaller than the number 2
N
F of low-level basis functions over ∂
𝒟
k, the calculation of [
S
~
C
C
] requires inverting a relatively small matrix—which can be carried out with direct solvers.In conclusion, the EEM facilitates the inversion of[
S
]
-
1, as anticipated. Once the expansion coefficients [
q
~
k
s
] in the eigencurrents basis have been computed, [
q
k
s
] can be retrieved through [
q
k
s
]
=
[
V
k
k
]
[
q
~
k
s
].


## 3.1. Reduction to a Weak Form with MoM

We employ the MoM with subsectional basis functions and a symmetric inner product repeatedly in order to (a) turn (2) and (6) into a weak form; and (b) formally invert (4), (8), and (14). In particular, we model the bricks’ boundaries and the volume of the anisotropic objects with 3D triangular-faceted tessellations and tetrahedral meshes, respectively. Then, we associate 2
N
F (N
O) surface (volume) divergence-conforming linear vector functions [21, 22] with the surface (volume) mesh on (inside) each brick in order to expand the surface currents J
k
s
,
i and M
k
s
,
i (the flux density D
o or B
o).As a result, we can write the algebraic counterparts ofS
k
k and T
k
n as
(20)
[
S
k
k
]
=
[
P
k
k
s
]
-
1
[
P
k
o
]
[
X
o
o
]
-
1
[
P
k
o
]
,
(21)
[
T
k
n
]
=
[
P
k
k
i
]
-
1
[
P
k
n
]
,
n
≠
k
,
where [
S
k
k
] and [
T
k
n
] have size 2
N
F
×
2
N
F and the rank of [
X
o
o
] is N
O.Upon using (20) and (21) as building blocks, the algebraic counterpart of (17) can be written formally as
(22)
[
S
]
-
1
[
q
s
]
=
[
q
i
]
,
where the 2
N
F
N
D
×
1 vectors [
q
s
,
i
] store the MoM expansion coefficients of q
k
s
,
i. As anticipated, though, explicitly building and then numerically inverting the linear system (22) are far from being trivial for two reasons.First and foremost, the inverse scattering operators[
S
k
k
]
-
1 may not be defined in some cases. Specifically, if the volume basis functions N
O in the tetrahedral mesh are fewer than the number 2
N
F of surface basis functions on ∂
𝒟
k, then, as is clear from (20), the scattering operator issingular, with its rank being min
⁡
{
2
N
F
,
N
O
}
=
N
O. At any rate, even though 2
N
F
≤
N
O, [
S
k
k
] may still be nearly singular due to the presence of a few very small eigenvalues. On the whole, inverting [
S
k
k
] is not only inconvenient, but also a potentially ill-conditioned problem.Secondly, difficulties posed by the inversion of[
S
k
k
] aside, the size of the total inverse scattering operator, may be quite large, if one considers a problem modelled with a moderate to large number of EM bricks. As a consequence, inversion with direct solvers such as LU factorization [23] may not be feasible. Then, iterative solvers seem an obvious alternative, but the convergence rate is expected to be poor, since [
S
]
-
1 is ill conditioned, in the light of the previous observations.To circumvent the hurdles above we apply the EEM [14].


## 3.2. Compression and Inversion with Eigencurrents

In the EEM we employ the eigenvectors[
V
k
k
] of [
S
k
k
] as local entire-domain basis functions to expand the scattered currents [
q
k
s
] on ∂
𝒟
k. Since these eigenvectors in tandem with the low-level basis functions define surface currents over ∂
𝒟
k, we refer to the columns of [
V
k
k
] as to the eigencurrents of 𝒟
k.The eigencurrents can be separated into two subsets:coupled anduncoupled. The former—associated with the larger eigenvalues of [
S
k
k
]—produce substantial radiation and thus contribute to the multiple scattering that occurs among the N
D bricks. Conversely, the uncoupled eigencurrents—which correspond with the smaller or null eigenvalues of [
S
k
k
]—yield scattered fields that are practically confined around a brick and hence serve to describe the scattered currents on a brick. We denote by N
C
≤
N
C
,
max
⁡
=
min
⁡
{
2
N
F
,
N
O
} the number of coupled eigencurrents.The notion of coupled-uncoupled eigencurrents allows us to partition the matrix[
S
]
-
1 (expressed in the basis of the eigencurrents) into four subblocks [14, Equation (
31
)]; namely,
(23)
[
S
^
]
-
1
=
[
[
S
~
C
C
]
-
1
[
S
~
C
U
]
-
1
[
S
~
U
C
]
-
1
[
S
~
U
U
]
-
1
]
,
where the subscripts “C” (“U”) stand for coupled (uncoupled).In the spectral reordered representation (23) of the total inverse scattering operator the subblock [
S
~
U
U
]
-
1—relevant to the interaction between any two uncoupled eigencurrents from either two distinct bricks or the same brick—can be approximated with its diagonal, that is, the reciprocal of the eigenvalues associated with all the uncoupled eigencurrents, because the latter interact only with themselves. Likewise, the off-diagonal subblocks [
S
~
C
U
]
-
1, [
S
~
U
C
]
-
1 can be zeroed altogether, as they account for the negligible interaction between a pair of coupled and uncoupled eigencurrents. Lastly, the diagonal subblock [
S
~
C
C
]
-
1—relative to the interaction of any two coupled eigencurrents from either two distinct bricks or the same brick—has to be retained in full.In the end, the inverse of[
S
^
]
-
1 is approximately given by
(24)
[
S
^
]
≈
[
[
S
~
C
C
]
[
0
]
[
0
]
[
Λ
U
U
]
]
,
where [
Λ
U
U
] is a diagonal matrix that contains the (
2
N
F
-
N
C
)
×
N
D eigenvalues germane to the uncoupled eigencurrents of the problem. Since [
Λ
U
U
] is known the moment the eigenvalues of [
S
k
k
] have been determined, evidently, the occurrence of small or null eigenvalues does not pose any problem at all in (24). Furthermore, since the number N
C of coupled eigencurrents contributed by 𝒟
k is usually far smaller than the number 2
N
F of low-level basis functions over ∂
𝒟
k, the calculation of [
S
~
C
C
] requires inverting a relatively small matrix—which can be carried out with direct solvers.In conclusion, the EEM facilitates the inversion of[
S
]
-
1, as anticipated. Once the expansion coefficients [
q
~
k
s
] in the eigencurrents basis have been computed, [
q
k
s
] can be retrieved through [
q
k
s
]
=
[
V
k
k
]
[
q
~
k
s
].


## 4. Results

### 4.1. Example of Validation

An extended numerical code has been developed for the calculation of[
S
k
k
] and [
T
n
k
] and the subsequent solution of (22) with the EEM in the instance of N
D
≥
1 objects with either dyadic permittivity or dyadic permeability.Among many other tests, as a sanity check we have examined the problem of a solitary dielectric sphere withscalar permittivity and embedded it in a cubic brick 𝒟
1. We have computed the scattering operator [
S
11
] with the newly developed code and the old one [14], which employs surface integral equations (e.g., the Poggio-Miller-Chang-Harrinton-Wu-Tai formulation (PMCHWT) [24]) at the interface between media ② and ③. With the two codes we have obtained the equivalent scattered currents q
1
s on ∂
𝒟
1 and the scattered fields in the Fraunhofer region when the sphere is illuminated with a plane wave. The number of surface basis functions over ∂
𝒟
1 is 2
N
F = 1152; the number of basis functions for the solution of the scattering problem inside 𝒟
1 is N
O
=
2
×
552 for the PMCHWT approach and N
O
=
1454 for the approach based on (4) and (11). (Geometrical and physical data are given in the caption of Figure 3.)Example of validation in near field: total equivalent current densitiesM
1
,
tot
i, J
1
,
tot
i, over a brick that embeds a dielectric sphere (not shown) in free space; (a), (c) results obtained with LEGO and PMCHWT; (b), (d) results obtained with LEGO and VIE (4). (See insets in Figure 4 for the surface and volumetric meshes of the sphere.)Data: brick’s edge 10 cm, sphere’s radius 4 cm, ɛ
1
=
ɛ
2
=
ɛ
0, ɛ
¯
3
=
2
ɛ
0
I
¯, E
i
=
x
^
exp
⁡
(
-
j
2
π
f
z
/
c
0
) V/m, and f
=
1 GHz.
(a)
|
M
1
,
tot
| LEGO/SIE
(b)
|
M
1
,
tot
| LEGO/VIE
(c)
|
J
1
,
tot
| LEGO/SIE
(d)
|
J
1
,
tot
| LEGO/VIEThe total equivalent surface current densitiesM
1
,
tot
=
M
1
i
-
M
1
s and J
1
,
tot
=
J
1
i
-
J
1
s yielded by the two codes are plotted in Figure 3; as can be seen, the current distributions are in excellent agreement. This notable result is a consequence of the algebraic scattering operator [
S
11
] being theoretically independent of the actual formulation of the EM problem inside 𝒟
1. Besides, if the scattered currents [
q
1
s
] are correct, then we may expect also the far-field quantities (which from [
q
1
s
] are derived) to be in excellent agreement. The latter is confirmed by the plot of the bistatic radar cross-section (RCS) given in Figure 4.Figure 4
Example of validation in far field: bistatic RCS of a dielectric sphere embedded in a cubic brick; (−) solution obtained through LEGO and PMCHWT [14], (∘) solution obtained through LEGO and VIE (14). Insets: surface and volumetric meshes of the sphere. (See Figure 3 for the surface of mesh of the cubic brick.)Data: brick’s edge 10 cm, sphere’s radius 4 cm, ɛ
1
=
ɛ
2
=
ɛ
0, ɛ
¯
3
=
2
ɛ
0
I
¯, E
i
=
x
^
exp
⁡
(
-
j
2
π
f
z
/
c
0
) V/m, and f
= 1 GHz.As a final remark, we notice that the solution with VIE will usually require more basis functionsN
O than what is necessary to solve the SIE (when the latter makes sense). This is obviously due to the need for distributing basis functions inside a volume rather than a surface. In fact, our numerical tests show that in order to attain the same level of accuracy of LEGO/SIE and LEGO/VIE, the size of the triangular facets over the sphere’s surface (see insets of Figure 4) should be comparable. If that criterion is met, then the requirement for having a uniform mesh in 𝒱
o somehow constrains the number of tetrahedra and hence basis functions in the object.


### 4.2.  Convergence with Number of Coupled Eigencurrents

To test the efficiency of the EEM [14] in compressing the algebraic system (22) we have considered the plane-wave scattering from a linear arrangement of N
D
=
4 anisotropic spheres. The EM problem has been modelled with N
D
=
4 cubic bricks each containing one sphere, as shown in Figure 5. The number of surface basis functions over ∂
𝒟
1 is 2
N
F = 1152, whereas the number of volume basis functions within the sphere is N
O
=
1454. (Geometrical data and constitutive parameters are given in the caption of Figure 5.)Figure 5
For convergence study, LEGO model withN
D
=
4 bricks of a linear arrangement of four dielectric spheres.Data: brick’s edge 10 cm, sphere’s radius 4 cm, ɛ
1
=
ɛ
2
=
ɛ
0, and ɛ
¯
3
=
(
4
x
^
x
^
+
2
y
^
y
^
+
2
z
^
z
^
)
ɛ
0.Notice that in the LEGO setup for this problem, only[
S
11
] needs to be computed, as the four bricks are identical. Furthermore, thanks to the limited translational symmetry of the structure in Figure 5 only N
T
=
6 transfer operators out of N
T
,
max
⁡
=
12 are independent and must be computed. This would not be the case, if the bricks had different shapes.We have obtained a reference solution by solving (22) with no reduction at all, that is, with N
C
,
max
⁡
=
2
N
F coupled eigencurrents. Subsequently, we have applied the EEM to its fullest, with N
C
∈
{
40,80
} coupled eigencurrents, and have solved the problem by using the compressed scattering operator (24). The RCSs computed with N
C
∈
{
40,80,2
N
F
} are plotted in Figure 6 for comparison.Figure 6
For convergence study, bistatic RCS of a linear arrangement of four dielectric spheres (see Figure5) for increasing number N
C of coupled eigencurrents. Inset: closeup of the RCS.Data: E
i
=
z
^
exp
⁡
(
-
j
2
π
f
x
/
c
0
) V/m, f
= 1 GHz.It can be seen thatN
C
=
40 coupled eigencurrents (
⊳
) are already sufficient for the solution to capture the trend and most details of the RCS versus the elevation angle θ; in particular, position of nulls and local maxima are predicted, even though the levels thereof depart from the reference values. However, with N
C
=
80 eigencurrents the RCS (∘) is practically indistinguishable from the reference curve (−). From a numerical viewpoint this means an accurate solution can be achieved by inverting a matrix of rank 80
×
4
=
320 as compared to the rank 1152
×
4
=
4608 of the original uncompressed system (22).We elaborate a bit further on this numerical experiment as regards the assumptions that constitute the groundwork of the reduction scheme based on the notion of coupled-uncoupled eigencurrents. In Section3.2 we posit that the coupled (uncoupled) eigencurrents correspond with the larger (smaller) eigenvalues λ
p
(
k
), p
=
1
,
…
,
2
N
F, of the (algebraic) scattering operator [
S
k
k
]. As an example, for the problem of Figure 5 the spectrum of [
S
11
] is plotted in Figure 7. It can be ascertained that the eigenvalues decrease exponentially to zero (i.e., the threshold of numerical noise in double-precision floating point arithmetics). In fact, |
λ
40
(
1
)
| and |
λ
80
(
1
)
| differ by approximately one order of magnitude. In the light of the conspicuous, fast decay of the spectrum, it can be understood why the eigencurrents associated with the eigenvalues λ
p
(
1
), p
>
80 contribute naught to the fields scattered by 𝒟
1 and hence can rightly be termed uncoupled.Figure 7
For convergence study, eigenvalues of[
S
k
k
], k
=
1
,
…
,
4, as a function of their index p
=
1
,
…
,
2
N
F for the bricks of the model in Figure 5 at f
= 1 GHz.In [18] we demonstrated that forisotropic objects the spectrum of [
S
k
k
] is mostly affected by the constitutive parameters of medium ③ (if media ① and ② are the same), the size and the position of the object with respect to ∂
𝒟
k; a weaker dependence on the working frequency was observed. For the time being, we mention that it seems reasonable to expect similar trends also in the case of dyadic constitutive parameters of medium ③, although a thorough and extensive numerical campaign is required to support this supposition. (Also see Section 4.3.)


### 4.3.  Example of Application: Arrays of Microplasma Rods

We now apply LEGO and EEM to solve a relatively large and multiscale problem, namely, the scattering from arrays of microplasma rods; these structures appear to have found application as reconfigurable metamaterials (e.g., [25]).We consider four arrays comprised of2
×
20, 4
×
20, 6
×
20, and 8
×
20
z-aligned cylindrical rods arranged in a Cartesian lattice in the x
O
y plane; shown in Figure 8 are the tetrahedral model of a rod and the 3D surface mesh of the cuboidal embedding brick. The lattice constants in both x and y directions are the same as a brick’s width. (See the caption of Figure 8 for applicable geometrical and mesh data.) Modelling the small cylindrical rods requires a reasonably fine mesh capable of accurately reproducing the small geometrical features. By contrast, the surface mesh on a brick’s boundary need not be as fine because the brick has got a simple shape. Thereby, we anticipate that the bare introduction of LEGO bricks can already afford reduction of unknowns.Modelling arrays of microplasma rods with LEGO: (a) volumetric mesh of a rod; (b) surface mesh of the brick that embeds the rod.Data: rod’s diameter 1 mm, rod’s height 10 mm, brick’s width 2.5 mm, brick’s height 11.5 mm, N
O
=
1876, and 2
N
F
=
1056.
(a)
(b)The plasma is a collisional magnetized gas of electrons with densityn
e
=
1
0
19 m
-
3, z-aligned static confining magnetic intensity B
z
=
1mT, and collision frequency ν
e
=
ω
p
e
/
2, where ω
p
e denotes the electron (angular) plasma frequency. With these values and for frequencies f
=
{
10,20,30
} GHz, the plasma dyadic permittivity in Cartesian coordinates reads
(25)
ɛ
¯
3
=
ɛ
0
(
S
j
D
0
-
j
D
S
0
0
0
P
)
,
where the Stix parameters S, D, and P [26] are given explicitly in Table 1. Finally, as the plasma rods are immersed in free space, then ɛ
1
=
ɛ
2
=
ɛ
0.Table 1
Plasma parameters for the problem of Figures8, 9, and 10.
f (GHz)
S
D
P
10
-
1.6735
-
j
3.7954
(
25
-
j
70
)
·
1
0
-
4
-
1.6735
-
j
3.7954
20
-
0.3402
-
j
0.9513
(
-
6
-
j
18
)
·
1
0
-
4
-
0.3402
-
j
0.9513
30
0.2681
-
j
0.3463
(
-
4.33
-
j
5.28
)
·
1
0
-
4
0.2681
-
j
0.3463By applying the EEM to (22) with N
C
=
50 coupled eigencurrents we have computed the equivalent scattered currents [
q
s
] and afterwards the forward and backward RCS for the four arrays of plasma rods. The results are plotted in Figure 9 as a function of the number of parallel rows along the x-axis for an x-directed z-polarized incident plane wave. (See the caption of Figure 9 for relevant data.) The minimum of reflection appears to occur for the array of 6
×
20 rods at f
=
10GHz, which is well below the electron plasma frequency ω
p
e
/
(
2
π
)
≈
28GHz.Figure 9
Microplasma rods: forward and backward RCS (along thex-direction) as a function of the number of rows (each is comprised of 20 plasma rods) for three values of frequency. Remark: the markers are joined with dotted lines for visualization’s sake.Data: E
i
=
z
^
exp
⁡
(
-
j
2
π
f
x
/
c
0
) V/m, f
∈
{
10,20,30
} GHz, and d
=
2.5mm (lattice constant and bricks’ width).The spectrum of[
S
11
] for the brick in Figure 8 at the three frequencies of interest is plotted in Figure 10; the vertical dashed line (at p
=
N
C
=
50) separates the regions of eigenvalues corresponding with the coupled (left) and the uncoupled (right) eigencurrents. Based on the convergence properties examined in Section 4.2 and the decay pattern which is not dissimilar from the one shown in Figure 7, we can assume that N
C
=
50 coupled eigencurrents are sufficient to determine accurate far-field quantities.Figure 10
Microplasma rods: eigenvalues of[
S
11
], as a function of their index p
=
1
,
…
,
2
N
F relative to the brick of the model in Figure 8 for three frequency values. Inset: closeup of the spectrum showing trend as frequency is increased.Besides, the plot of Figure10 seems to confirm that the spectrum of a scattering operator is only weakly dependent on frequency [18] even in the case of embedded anisotropic objects. In particular, the first eigenvalues are seen to become smaller in magnitude as the frequency increases. This behavior can be explained intuitively by observing that on average the electrical distance between a rod and its surrounding brickincreases as the frequency does so. Now, since the scattered field decays away from the object that has generated it and in the basis of the eigencurrents the scattered field (for a single brick) is given by [
q
~
p
s
]
=
λ
p
(
1
)
[
q
~
p
i
], then it follows that |
λ
p
(
1
)
| must be smaller if the distance from the brick’s surface and the included object is larger, everything else unchanged.


### 4.4.  Efficiency

We conclude this section by discussing the time requirements of LEGO-EEM with the aid of the problems of Figures5 and 8. For the sake of completeness we mention that the simulations for the spheres (Figure 5) have been run on a PC equipped with an Intel Core i7-3770 CPU 3.40-GHz processor and 8-GB RAM, whereas the simulations for the microplasma rods (Figure 8) have been run on a PC equipped with an Intel Core Quad CPU 2.66-GHz processor and 4-GB RAM.In the extension of LEGO to anisotropic bodies the calculation of the scattering operator (20) appears to be relatively time consuming, owing to the need for computing volume integrals in tetrahedra and (formally) inverting [
X
o
o
], which is usually larger than that in case of SIEs. Therefore, the time T
[
S
k
k
] taken to obtain [
S
k
k
] appears to be roughly proportional to N
O
2.For the problems at hand,T
[
S
11
] is about 49 min for the anisotropic sphere case and about 59 min for the plasma rod case. Nevertheless, the calculation of [
S
k
k
] (and its spectral decomposition) has to be performed only once, if the structure can be modelled by means of N
D identical bricks, as we have shown. This is especially useful if one plans on reutilizing the bricks over again in other simulations, as long as the frequency is not altered, of course.When this expedient is applied to the case of the anisotropic spheres with different numbers of coupled eigencurrents, the overall calculation timeT
LEGO actually becomes T
LEGO
-
T
[
S
11
], which drops drastically, as can be ascertained from Table 2.Table 2
Computation times for the problem of Figure5.
N
C
N
C
N
D
T
[
S
11
]
T
LEGO
T
LEGO
-
T
[
S
11
]
1152
4608
49 min
65 min
15 min†
40
160
—
54 min
≈4 min
80
320
—
53 min†
≈4 min
†Times derived indirectly, as these simulations have not been run as such.Likewise, we have to compute[
S
k
k
] only once (k
=
1) for the four arrays of microplasma rods. The calculation times T
LEGO
-
T
[
S
11
] are listed in Table 3. Interestingly, T
LEGO
-
T
[
S
11
] grows only linearly with the number of bricks N
D, as shown by the plot in Figure 11. This result is a consequence of the limited translational symmetry of the arrays; in practice, the transfer operators [
T
k
n
] and [
T
n
k
] come in clusters of identical specimens and hence it is sufficient to compute only one element of each cluster [14, 17]. Then, as long as the number of clusters of transfer operators grows linearly with N
D, so does the time T
LEGO
-
T
[
S
11
]. It is worthwhile noticing, though, that if one examines different EM problems involving afixed number N
D of objects with increasingly larger size (i.e., ever more unknowns N
O), the complexity—being dominated by the calculation of [
S
11
]—will rather grow quadratically.Table 3
Computation times for the problem of Figures8 and 9.
N
O
N
D
2
N
F
N
D
N
C
N
D
T
LEGO
-
T
[
S
11
]
37520
21120
1000
17 min
75040
42240
2000
34 min
112560
63360
3000
52 min
150080
84480
4000
73 minFigure 11
Microplasma rods: computational times minus the time taken to compute[
S
11
] for the problems of Figures 8 and 9 at f
=
10GHz.Finally, the matrix sizes relative to the problem of arrays of microplasma rods are listed in Table3: N
O
N
D is the size of algebraic system that one would need to invert if the problem was formulated directly with a single VIE; 2
N
F
N
D is the size of the LEGO system (22); N
C
N
D is the size of [
S
~
C
C
]
-
1 in (23). It is apparent that the LEGO approach is already competitive over the baseline VIE, whereas further and substantial reduction is afforded by the EEM.


## 4.1. Example of Validation

An extended numerical code has been developed for the calculation of[
S
k
k
] and [
T
n
k
] and the subsequent solution of (22) with the EEM in the instance of N
D
≥
1 objects with either dyadic permittivity or dyadic permeability.Among many other tests, as a sanity check we have examined the problem of a solitary dielectric sphere withscalar permittivity and embedded it in a cubic brick 𝒟
1. We have computed the scattering operator [
S
11
] with the newly developed code and the old one [14], which employs surface integral equations (e.g., the Poggio-Miller-Chang-Harrinton-Wu-Tai formulation (PMCHWT) [24]) at the interface between media ② and ③. With the two codes we have obtained the equivalent scattered currents q
1
s on ∂
𝒟
1 and the scattered fields in the Fraunhofer region when the sphere is illuminated with a plane wave. The number of surface basis functions over ∂
𝒟
1 is 2
N
F = 1152; the number of basis functions for the solution of the scattering problem inside 𝒟
1 is N
O
=
2
×
552 for the PMCHWT approach and N
O
=
1454 for the approach based on (4) and (11). (Geometrical and physical data are given in the caption of Figure 3.)Example of validation in near field: total equivalent current densitiesM
1
,
tot
i, J
1
,
tot
i, over a brick that embeds a dielectric sphere (not shown) in free space; (a), (c) results obtained with LEGO and PMCHWT; (b), (d) results obtained with LEGO and VIE (4). (See insets in Figure 4 for the surface and volumetric meshes of the sphere.)Data: brick’s edge 10 cm, sphere’s radius 4 cm, ɛ
1
=
ɛ
2
=
ɛ
0, ɛ
¯
3
=
2
ɛ
0
I
¯, E
i
=
x
^
exp
⁡
(
-
j
2
π
f
z
/
c
0
) V/m, and f
=
1 GHz.
(a)
|
M
1
,
tot
| LEGO/SIE
(b)
|
M
1
,
tot
| LEGO/VIE
(c)
|
J
1
,
tot
| LEGO/SIE
(d)
|
J
1
,
tot
| LEGO/VIEThe total equivalent surface current densitiesM
1
,
tot
=
M
1
i
-
M
1
s and J
1
,
tot
=
J
1
i
-
J
1
s yielded by the two codes are plotted in Figure 3; as can be seen, the current distributions are in excellent agreement. This notable result is a consequence of the algebraic scattering operator [
S
11
] being theoretically independent of the actual formulation of the EM problem inside 𝒟
1. Besides, if the scattered currents [
q
1
s
] are correct, then we may expect also the far-field quantities (which from [
q
1
s
] are derived) to be in excellent agreement. The latter is confirmed by the plot of the bistatic radar cross-section (RCS) given in Figure 4.Figure 4
Example of validation in far field: bistatic RCS of a dielectric sphere embedded in a cubic brick; (−) solution obtained through LEGO and PMCHWT [14], (∘) solution obtained through LEGO and VIE (14). Insets: surface and volumetric meshes of the sphere. (See Figure 3 for the surface of mesh of the cubic brick.)Data: brick’s edge 10 cm, sphere’s radius 4 cm, ɛ
1
=
ɛ
2
=
ɛ
0, ɛ
¯
3
=
2
ɛ
0
I
¯, E
i
=
x
^
exp
⁡
(
-
j
2
π
f
z
/
c
0
) V/m, and f
= 1 GHz.As a final remark, we notice that the solution with VIE will usually require more basis functionsN
O than what is necessary to solve the SIE (when the latter makes sense). This is obviously due to the need for distributing basis functions inside a volume rather than a surface. In fact, our numerical tests show that in order to attain the same level of accuracy of LEGO/SIE and LEGO/VIE, the size of the triangular facets over the sphere’s surface (see insets of Figure 4) should be comparable. If that criterion is met, then the requirement for having a uniform mesh in 𝒱
o somehow constrains the number of tetrahedra and hence basis functions in the object.


## 4.2.  Convergence with Number of Coupled Eigencurrents

To test the efficiency of the EEM [14] in compressing the algebraic system (22) we have considered the plane-wave scattering from a linear arrangement of N
D
=
4 anisotropic spheres. The EM problem has been modelled with N
D
=
4 cubic bricks each containing one sphere, as shown in Figure 5. The number of surface basis functions over ∂
𝒟
1 is 2
N
F = 1152, whereas the number of volume basis functions within the sphere is N
O
=
1454. (Geometrical data and constitutive parameters are given in the caption of Figure 5.)Figure 5
For convergence study, LEGO model withN
D
=
4 bricks of a linear arrangement of four dielectric spheres.Data: brick’s edge 10 cm, sphere’s radius 4 cm, ɛ
1
=
ɛ
2
=
ɛ
0, and ɛ
¯
3
=
(
4
x
^
x
^
+
2
y
^
y
^
+
2
z
^
z
^
)
ɛ
0.Notice that in the LEGO setup for this problem, only[
S
11
] needs to be computed, as the four bricks are identical. Furthermore, thanks to the limited translational symmetry of the structure in Figure 5 only N
T
=
6 transfer operators out of N
T
,
max
⁡
=
12 are independent and must be computed. This would not be the case, if the bricks had different shapes.We have obtained a reference solution by solving (22) with no reduction at all, that is, with N
C
,
max
⁡
=
2
N
F coupled eigencurrents. Subsequently, we have applied the EEM to its fullest, with N
C
∈
{
40,80
} coupled eigencurrents, and have solved the problem by using the compressed scattering operator (24). The RCSs computed with N
C
∈
{
40,80,2
N
F
} are plotted in Figure 6 for comparison.Figure 6
For convergence study, bistatic RCS of a linear arrangement of four dielectric spheres (see Figure5) for increasing number N
C of coupled eigencurrents. Inset: closeup of the RCS.Data: E
i
=
z
^
exp
⁡
(
-
j
2
π
f
x
/
c
0
) V/m, f
= 1 GHz.It can be seen thatN
C
=
40 coupled eigencurrents (
⊳
) are already sufficient for the solution to capture the trend and most details of the RCS versus the elevation angle θ; in particular, position of nulls and local maxima are predicted, even though the levels thereof depart from the reference values. However, with N
C
=
80 eigencurrents the RCS (∘) is practically indistinguishable from the reference curve (−). From a numerical viewpoint this means an accurate solution can be achieved by inverting a matrix of rank 80
×
4
=
320 as compared to the rank 1152
×
4
=
4608 of the original uncompressed system (22).We elaborate a bit further on this numerical experiment as regards the assumptions that constitute the groundwork of the reduction scheme based on the notion of coupled-uncoupled eigencurrents. In Section3.2 we posit that the coupled (uncoupled) eigencurrents correspond with the larger (smaller) eigenvalues λ
p
(
k
), p
=
1
,
…
,
2
N
F, of the (algebraic) scattering operator [
S
k
k
]. As an example, for the problem of Figure 5 the spectrum of [
S
11
] is plotted in Figure 7. It can be ascertained that the eigenvalues decrease exponentially to zero (i.e., the threshold of numerical noise in double-precision floating point arithmetics). In fact, |
λ
40
(
1
)
| and |
λ
80
(
1
)
| differ by approximately one order of magnitude. In the light of the conspicuous, fast decay of the spectrum, it can be understood why the eigencurrents associated with the eigenvalues λ
p
(
1
), p
>
80 contribute naught to the fields scattered by 𝒟
1 and hence can rightly be termed uncoupled.Figure 7
For convergence study, eigenvalues of[
S
k
k
], k
=
1
,
…
,
4, as a function of their index p
=
1
,
…
,
2
N
F for the bricks of the model in Figure 5 at f
= 1 GHz.In [18] we demonstrated that forisotropic objects the spectrum of [
S
k
k
] is mostly affected by the constitutive parameters of medium ③ (if media ① and ② are the same), the size and the position of the object with respect to ∂
𝒟
k; a weaker dependence on the working frequency was observed. For the time being, we mention that it seems reasonable to expect similar trends also in the case of dyadic constitutive parameters of medium ③, although a thorough and extensive numerical campaign is required to support this supposition. (Also see Section 4.3.)


## 4.3.  Example of Application: Arrays of Microplasma Rods

We now apply LEGO and EEM to solve a relatively large and multiscale problem, namely, the scattering from arrays of microplasma rods; these structures appear to have found application as reconfigurable metamaterials (e.g., [25]).We consider four arrays comprised of2
×
20, 4
×
20, 6
×
20, and 8
×
20
z-aligned cylindrical rods arranged in a Cartesian lattice in the x
O
y plane; shown in Figure 8 are the tetrahedral model of a rod and the 3D surface mesh of the cuboidal embedding brick. The lattice constants in both x and y directions are the same as a brick’s width. (See the caption of Figure 8 for applicable geometrical and mesh data.) Modelling the small cylindrical rods requires a reasonably fine mesh capable of accurately reproducing the small geometrical features. By contrast, the surface mesh on a brick’s boundary need not be as fine because the brick has got a simple shape. Thereby, we anticipate that the bare introduction of LEGO bricks can already afford reduction of unknowns.Modelling arrays of microplasma rods with LEGO: (a) volumetric mesh of a rod; (b) surface mesh of the brick that embeds the rod.Data: rod’s diameter 1 mm, rod’s height 10 mm, brick’s width 2.5 mm, brick’s height 11.5 mm, N
O
=
1876, and 2
N
F
=
1056.
(a)
(b)The plasma is a collisional magnetized gas of electrons with densityn
e
=
1
0
19 m
-
3, z-aligned static confining magnetic intensity B
z
=
1mT, and collision frequency ν
e
=
ω
p
e
/
2, where ω
p
e denotes the electron (angular) plasma frequency. With these values and for frequencies f
=
{
10,20,30
} GHz, the plasma dyadic permittivity in Cartesian coordinates reads
(25)
ɛ
¯
3
=
ɛ
0
(
S
j
D
0
-
j
D
S
0
0
0
P
)
,
where the Stix parameters S, D, and P [26] are given explicitly in Table 1. Finally, as the plasma rods are immersed in free space, then ɛ
1
=
ɛ
2
=
ɛ
0.Table 1
Plasma parameters for the problem of Figures8, 9, and 10.
f (GHz)
S
D
P
10
-
1.6735
-
j
3.7954
(
25
-
j
70
)
·
1
0
-
4
-
1.6735
-
j
3.7954
20
-
0.3402
-
j
0.9513
(
-
6
-
j
18
)
·
1
0
-
4
-
0.3402
-
j
0.9513
30
0.2681
-
j
0.3463
(
-
4.33
-
j
5.28
)
·
1
0
-
4
0.2681
-
j
0.3463By applying the EEM to (22) with N
C
=
50 coupled eigencurrents we have computed the equivalent scattered currents [
q
s
] and afterwards the forward and backward RCS for the four arrays of plasma rods. The results are plotted in Figure 9 as a function of the number of parallel rows along the x-axis for an x-directed z-polarized incident plane wave. (See the caption of Figure 9 for relevant data.) The minimum of reflection appears to occur for the array of 6
×
20 rods at f
=
10GHz, which is well below the electron plasma frequency ω
p
e
/
(
2
π
)
≈
28GHz.Figure 9
Microplasma rods: forward and backward RCS (along thex-direction) as a function of the number of rows (each is comprised of 20 plasma rods) for three values of frequency. Remark: the markers are joined with dotted lines for visualization’s sake.Data: E
i
=
z
^
exp
⁡
(
-
j
2
π
f
x
/
c
0
) V/m, f
∈
{
10,20,30
} GHz, and d
=
2.5mm (lattice constant and bricks’ width).The spectrum of[
S
11
] for the brick in Figure 8 at the three frequencies of interest is plotted in Figure 10; the vertical dashed line (at p
=
N
C
=
50) separates the regions of eigenvalues corresponding with the coupled (left) and the uncoupled (right) eigencurrents. Based on the convergence properties examined in Section 4.2 and the decay pattern which is not dissimilar from the one shown in Figure 7, we can assume that N
C
=
50 coupled eigencurrents are sufficient to determine accurate far-field quantities.Figure 10
Microplasma rods: eigenvalues of[
S
11
], as a function of their index p
=
1
,
…
,
2
N
F relative to the brick of the model in Figure 8 for three frequency values. Inset: closeup of the spectrum showing trend as frequency is increased.Besides, the plot of Figure10 seems to confirm that the spectrum of a scattering operator is only weakly dependent on frequency [18] even in the case of embedded anisotropic objects. In particular, the first eigenvalues are seen to become smaller in magnitude as the frequency increases. This behavior can be explained intuitively by observing that on average the electrical distance between a rod and its surrounding brickincreases as the frequency does so. Now, since the scattered field decays away from the object that has generated it and in the basis of the eigencurrents the scattered field (for a single brick) is given by [
q
~
p
s
]
=
λ
p
(
1
)
[
q
~
p
i
], then it follows that |
λ
p
(
1
)
| must be smaller if the distance from the brick’s surface and the included object is larger, everything else unchanged.


## 4.4.  Efficiency

We conclude this section by discussing the time requirements of LEGO-EEM with the aid of the problems of Figures5 and 8. For the sake of completeness we mention that the simulations for the spheres (Figure 5) have been run on a PC equipped with an Intel Core i7-3770 CPU 3.40-GHz processor and 8-GB RAM, whereas the simulations for the microplasma rods (Figure 8) have been run on a PC equipped with an Intel Core Quad CPU 2.66-GHz processor and 4-GB RAM.In the extension of LEGO to anisotropic bodies the calculation of the scattering operator (20) appears to be relatively time consuming, owing to the need for computing volume integrals in tetrahedra and (formally) inverting [
X
o
o
], which is usually larger than that in case of SIEs. Therefore, the time T
[
S
k
k
] taken to obtain [
S
k
k
] appears to be roughly proportional to N
O
2.For the problems at hand,T
[
S
11
] is about 49 min for the anisotropic sphere case and about 59 min for the plasma rod case. Nevertheless, the calculation of [
S
k
k
] (and its spectral decomposition) has to be performed only once, if the structure can be modelled by means of N
D identical bricks, as we have shown. This is especially useful if one plans on reutilizing the bricks over again in other simulations, as long as the frequency is not altered, of course.When this expedient is applied to the case of the anisotropic spheres with different numbers of coupled eigencurrents, the overall calculation timeT
LEGO actually becomes T
LEGO
-
T
[
S
11
], which drops drastically, as can be ascertained from Table 2.Table 2
Computation times for the problem of Figure5.
N
C
N
C
N
D
T
[
S
11
]
T
LEGO
T
LEGO
-
T
[
S
11
]
1152
4608
49 min
65 min
15 min†
40
160
—
54 min
≈4 min
80
320
—
53 min†
≈4 min
†Times derived indirectly, as these simulations have not been run as such.Likewise, we have to compute[
S
k
k
] only once (k
=
1) for the four arrays of microplasma rods. The calculation times T
LEGO
-
T
[
S
11
] are listed in Table 3. Interestingly, T
LEGO
-
T
[
S
11
] grows only linearly with the number of bricks N
D, as shown by the plot in Figure 11. This result is a consequence of the limited translational symmetry of the arrays; in practice, the transfer operators [
T
k
n
] and [
T
n
k
] come in clusters of identical specimens and hence it is sufficient to compute only one element of each cluster [14, 17]. Then, as long as the number of clusters of transfer operators grows linearly with N
D, so does the time T
LEGO
-
T
[
S
11
]. It is worthwhile noticing, though, that if one examines different EM problems involving afixed number N
D of objects with increasingly larger size (i.e., ever more unknowns N
O), the complexity—being dominated by the calculation of [
S
11
]—will rather grow quadratically.Table 3
Computation times for the problem of Figures8 and 9.
N
O
N
D
2
N
F
N
D
N
C
N
D
T
LEGO
-
T
[
S
11
]
37520
21120
1000
17 min
75040
42240
2000
34 min
112560
63360
3000
52 min
150080
84480
4000
73 minFigure 11
Microplasma rods: computational times minus the time taken to compute[
S
11
] for the problems of Figures 8 and 9 at f
=
10GHz.Finally, the matrix sizes relative to the problem of arrays of microplasma rods are listed in Table3: N
O
N
D is the size of algebraic system that one would need to invert if the problem was formulated directly with a single VIE; 2
N
F
N
D is the size of the LEGO system (22); N
C
N
D is the size of [
S
~
C
C
]
-
1 in (23). It is apparent that the LEGO approach is already competitive over the baseline VIE, whereas further and substantial reduction is afforded by the EEM.


## 5. Conclusion

We have shown that the extension of the LEGO procedure for solving the EM scattering from clusters of anisotropic bodies basically requires modifying the way the scattering operators are computed. In fact, as long as the intrinsic modularity of LEGO is properly exploited, then incorporating new types of EM bricks into an existing code is in principle straightforward. Besides, modularity implies that when calculations are suitably organized, the scattering operators can be reutilized in subsequent runs, thus enabling, for example, the efficient analysis of increasingly larger structures.Along the same lines, the EEM can be applied as before because its implementation is independent of the specific types of EM bricks being considered. In this regard, we have demonstrated that when anisotropic objects are enclosed in EM bricks, the EEM performs as efficiently as it does in the instance of conducting or penetrable isotropic media [18].All in all, the features above make the LEGO method an attractive technique for the numerical simulations of moderately large and yet finite-size structures comprised of many separate objects.We conclude by observing that although we have presented results for homogeneous anisotropic media, the procedure described in Section2.2 applies as well to the case of spatially dependent permittivity or permeability without any modification [20].



---

*Source: 467931-2014-03-19.xml*