Integral Equations in Electromagnetics

Integral Equations in Electromagnetics Massachusetts Institute of Technology 6.635 lecture notes Most integral equations do not have a closed form solution. However, they can often be discretized and solved on a digital computer. Proof of the existence of the solution to an integral equation by discretization was first presented by Fredholm in 1903. In general, integral equations can be divided into two families: 1. When the unknown is in the integral only, the integral equation is called of the first kind. 2. When the unknown is both inside and outside the integral, the integral equation is called of the second kind. For electromagnetic applications, we can have both scalar and vector integral equations. 1 Scalar integral equations Let us consider the situation of Fig. 1: two regions are defined in space, region 2 bounded by the closed surface S and region 1 being all the remaining space, bounded by S and S ∞ , in which sources are located. #1 J n̂ S #2 S∞ Figure 1: Definition of the geometry of the problem. It is already known that we can write: (∇2 + k12 ) Φ1 (r̄) = J(r̄) , (1a) (∇2 + k22 ) Φ2 (r̄) = 0 , (1b) 1 2 Section 1. Scalar integral equations and similarly for the Green’s functions (∇2 + k12 ) g1 (r̄, r̄0 ) = −δ(r̄ − r̄ 0 ) 2 (∇ + k22 ) g2 (r̄, r̄0 ) 0 = −δ(r̄ − r̄ ) in region 1, (2a) in region 2. (2b) Upon performing (Eq. (1a)g1 (r̄, r̄0 )−Eq. (2b)Φ1 (r̄, r̄0 )), we get: [g1 (r̄, r̄0 ) ∇2 Φ1 (r̄) − Φ1 (r̄) ∇2 g1 (r̄, r̄0 )] = J(r̄)g1 (r̄, r̄0 ) + δ(r̄ − r̄ 0 )Φ1 (r̄) . (3) Upon integrating Eq. (3) over the entire volume and using the identity ∇ · (g∇Φ − Φ∇g) = − Φ∇2 g, we get: Z Z 0 0 dvJ(r̄) g(r̄, r̄ 0 ) + Φ1 (r̄0 ) . (4) dv ∇ · [g1 (r̄, r̄ )Φ1 (r̄) − Φ1 (r̄)∇g1 (r̄, r̄ )] = g∇2 Φ V V By Gauss theorem, we reduce the left-hand side integral to a surface integral. Also Z dvJ(r̄)g1 (r̄, r̄0 ) = −Φinc (r̄0 ) . (5) V We therefore obtain Z ds n̂ · [g1 (r̄, r̄0 )∇Φ1 (r̄) − Φ1 (r̄)∇g1 (r̄, r̄0 )] = −Φinc (r̄) + Φ1 (r̄0 ) , − r̄0 ∈ V . (6) S+S∞ By invoking the radiation condition, the integral over S∞ vanishes, leaving (and exchanging r̄ and r̄0 so that primed coordinates correspond to sources and unprimed ones to observation): Φ1 (r̄) = Φinc (r̄) − Z ds0 n̂ · [g1 (r̄, r̄0 )∇0 Φ1 (r̄0 ) − Φ1 (r̄0 )∇0 g1 (r̄, r̄0 )] r̄ ∈ V1 . (7) S For r̄ ∈ V2 , the wave equation has no source and therefore the integration of the delta function yields a zero value. Performing the same steps for this second case, we get the generic relation:  Φ (r̄) r̄ ∈ V 1 1 ds0 n̂ · [g1 (r̄, r̄0 )∇0 Φ1 (r̄0 ) − Φ1 (r̄0 )∇0 g1 (r̄, r̄0 )] = Φinc (r̄) − 0 S r̄ ∈ V2 Z (8) This is directly evocative of Huygens’ principle: • For r̄ ∈ V1 : the total field Φ1 (r̄) is the sum of the incident field plus the field due to the surface currents on the surface S. • For r̄ ∈ V2 : the surface source on S produces a field that exactly opposes Φinc , yielding the extinction theorem. Applying the same reasoning to region 2, we write (where there is no incident field): Z ds0 n̂ · [g2 (r̄, r̄0 )∇0 Φ2 (r̄0 ) − Φ2 (r̄0 )∇0 g2 (r̄, r̄0 )] = S  0 r̄ ∈ V1 Φ2 (r̄) r̄ ∈ V2 (9) 3 Note that the sign reversal is due to the definition of the normal vector n̂ which has to point outward from the surface. Here, since we use the same n̂ as before, we have to take it as being negative. Eqs. (8) and (9) have four independent unknowns: Φ1 , Φ2 ; n̂ · ∇Φ1 , n̂ · ∇Φ2 , (10) which can be related by the boundary conditions. Here, g(r̄, r̄ 0 ) and n̂ · ∇g(r̄, r̄ 0 ) are the kernel of the integral equation. 2 Vector integral equation For the sake of completeness, we shall write the vector wave equation as well, although we will not use it here directly. Considering the same situation as above, we know that the fields have to satisfy: ¯ , ∇ × ∇ × Ē1 (r̄) − ω 2 ²1 µ1 Ē1 (r̄) =iωµ1 J(r̄) (11a) ∇ × ∇ × Ē2 (r̄) − ω 2 ²2 µ2 Ē2 (r̄) =0 , (11b) and the Green’s functions: ∇ × ∇ × G1 (r̄, r̄0 ) − ω 2 ²1 µ1 G1 (r̄, r̄0 ) = Iδ(r̄ − r̄ 0 ) , (12a) ∇ × ∇ × G2 (r̄, r̄0 ) − ω 2 ²2 µ2 G2 (r̄, r̄0 ) = Iδ(r̄ − r̄ 0 ) . (12b) R By the same technique as before ( V dv [Eq. (11a) · G1 (r̄, r̄0 ) − Eq. (12a) · Ē1 (r̄)]dv) ¸ · Z Z 0 0 ¯ · G1 (r̄, r̄0 ) − Ē1 (r̄) . dv J(r̄) dv ∇ × ∇ × Ē1 (r̄) · G1 (r̄, r̄ ) − Ē1 (r̄) · ∇ × ∇ × G1 (r̄, r̄ ) = iωµ1 V V (13) The left-hand side can be transformed into a surface integral (left as exercise) and the right-hand side written in terms of incident field, yielding: ½ ¾ Z 0 0 0 0 ds n̂· [∇× Ē1 (r̄)]×G1 (r̄, r̄ )+ Ē1 (r̄)×∇×G1 (r̄, r̄ ) , r̄0 ∈ V1 . (14) Ē1 (r̄ ) = Ēinc (r̄ )+ S+S∞ By reciprocity of the Green’s functions, we can transform n̂ · [∇ × Ē1 (r̄)] × G1 (r̄, r̄0 ) = n̂ × [∇ × Ē1 (r̄)] · G1 (r̄, r̄0 ) = iωµ1 G1 (r̄, r̄0 ) · n̂ × H̄1 (r̄) . (15) In addition, for an unbounded homogeneous medium (∇ × G(r̄, r̄ 0 ) is reciprocal): n̂ · Ē1 (r̄) × ∇ × G1 (r̄, r̄0 ) = n̂ × Ē1 (r̄) · ∇ × G1 (r̄, r̄0 ) = −[∇ × G1 (r̄, r̄0 )] · n̂ × Ē1 (r̄) . (16) 4 Section 3. Problem with the internal resonance For Green’s functions that satisfy the radiation condition, Eq. (14) becomes: ½ ¾ Z 0 0 0 0 ds n̂ · iωµ1 G1 (r̄, r̄ ) · n̂ × H̄1 (r̄) − [∇ × G1 (r̄, r̄ )] · n̂ × Ē1 (r̄) . Ē1 (r̄ ) = Ēinc (r̄ ) + (17) S We perform the same steps for the other region to eventually obtain (and again interchanging r̄ and r̄0 ):  ¾  ½ Ē1 (r̄) r̄ ∈ V1 Ēinc (r̄) + ds0 iωµ1 G1 (r̄, r̄0 ) · n̂ × H̄1 (r̄0 ) − [∇0 × G1 (r̄, r̄0 )] · n̂ × Ē1 (r̄0 ) = 0 S r̄ ∈ V2 Z (18a) − Z ¾ ½ 0 0 0 0 0 ds iωµ2 G2 (r̄, r̄ ) · n̂ × H̄2 (r̄ ) − [∇ × G2 (r̄, r̄ )] · n̂ × Ē2 (r̄ ) = 0 S  0 r̄ ∈ V1 Ē2 (r̄) r̄ ∈ V2 (18b) Together with the boundary conditions n̂ × H̄1 (r̄) = n̂ × H̄2 (r̄) , (19a) n̂ × Ē1 (r̄) = n̂ × Ē2 (r̄) , (19b) this system can be solved for n̂ × Ē1 (r̄) and n̂ × H̄1 (r̄). Note also that Eqs. (18) can be written in terms of electric and magnetic currents, and magnetic Green’s functions: J¯eq (r̄0 ) = n̂ × H̄1 (r̄0 ) ; −M̄eq (r̄0 ) = n̂ × Ē1 (r̄0 ) , (20a) so that for example:  ½ ¾  Ē1 (r̄) r̄ ∈ V1 ¯ 0 ) + Gm (r̄, r̄0 ) · M̄ (r̄0 ) = Ēinc (r̄) + ds0 iωµ1 Ge1 (r̄, r̄0 ) · J(r̄ 1 0 S r̄ ∈ V2 Z 3 (21) Problem with the internal resonance A question arises: is the integral equation equivalent to Maxwell’s equations? Or asked differently, if we solve the integral equation and Maxwell’s equations, do we get the same solution? The answer is actually “no”, that is they are not always equivalent to each-other. The problem comes from spurious solutions at frequencies corresponding to the eigenfrequencies of the cavity enclosed by the surface S. This problem is generally referred to as the “internal resonance of the integral equation”. However, this lack of complete equivalence between the physical problem and its defining integral equation is rather minor and infrequent phenomenon, and is therefore often tolerated in practice. 5 4 Scattering by a rough surface Let us consider the 2D problem (for which we shall use a scalar integral equation) depicted in Fig. 2. z ²0 f (x) PSfrag replacements x ²1 Figure 2: Rough surface S separating two media. The integral equation in scalar form is given by: Φinc (r̄) + 4.1 Z ds0 n̂ · [Φ(r̄ 0 )∇g(r̄, r̄ 0 ) − g(r̄, r̄ 0 )∇Φ(r̄0 )] = S Dirichlet boundary conditions: EFIE  Φ(r̄) 0 r̄ ∈ V0 r̄ ∈ V1 (22) For a TE wave (Ē = E ŷ) and PEC surface, the boundary condition is Φ(r̄) = 0 , for r̄ ∈ S . (23) The integral equation becomes, for r̄ ∈ S , r̄ 0 ∈ S: Φinc (r̄) − Z ds0 g(r̄, r̄0 ) n̂ · ∇Φ(r̄ 0 ) = S  Φ(r̄) (24) 0 This equation, in which Φ represent the electric field, is referred to as the electric field integral equation (EFIE). Note that as r̄ gets closer to the surface, Φ(r̄) → 0 (from the boundary condition) so that we do not need to distinguish between approaching the surface from one side or the other. In fact, we can unify the equations and write: Z Φinc (r̄) − ds0 g(r̄, r̄0 ) n̂ · ∇Φ(r̄ 0 ) = 0 , r̄ ∈ S , r̄ 0 ∈ S . (25) S g(r̄, r̄ 0 ) In addition, has an integrable singularity as r̄ → r̄ 0 . Let us consider the surface depicted in Fig. 2, with z = f (x): s µ ¶2 p df 2 2 ds = dx + dz = dx 1 + , (26) dx such that the integral equation becomes: s µ ¶2 Z df dx0 1 + Φinc (r̄) = g(x, f (x), x0 , f (x0 )) (n̂ · ∇Φ(r̄ 0 )) , dx ∆x at z 0 = f (x0 ) , (27) 6 4.2 Neumann boundary conditions: MFIE where we can limit ∆x to [−L/2, L/2]. By letting s µ ¶2 df (n̂ · ∇Φ(r̄ 0 ))|z 0 =f (x0 ) = U (x0 ) , 1+ dx (28a) Φinc (x, f (x)) = b(x) , (28b) K(x, x0 ) = g(x, f (x), x0 , f (x0 )) , (28c) we can rewrite the integral equation as Z L/2 dx0 K(x, x0 ) U (x0 ) = b(x) , (29) −L/2 which has to be solved numerically. Before doing this, we shall write the integral equation with Neumann boundary conditions. 4.2 Neumann boundary conditions: MFIE This corresponds to a TM case with H̄ = H ŷ with a PEC surface. In this case, the boundary condition is n̂ · ∇Φ(r̄) = 0 , (30) and the integral equation becomes  Φ(r̄) r̄ ∈ V 0 Φinc (r̄) + ds0 Φ(r̄) n̂ · ∇g(r̄, r̄ 0 ) = 0 S r̄ ∈ V1 Z (31) In this case, it makes a difference if we approach the surface from the top or from the bottom. In fact Z Φinc (r̄+ ) + ds0 Φ(r̄0 ) n̂ · ∇g(r̄+ , r̄0 ) = Φ(r̄+ ) , (32a) S Z Φinc (r̄− ) + ds0 Φ(r̄0 ) n̂ · ∇g(r̄− , r̄0 ) = 0 . (32b) S These two equations seems inconsistent with one another. We shall show that in fact, they are actually consistent with each-other, due to the singularity of the Green’s functions. We shall examine what happens when we let r̄ approach the surface. Fig. 3 is an illustration of the situation at the immediate vicinity of the surface. The integral part of the equation can be written as: Z Z Z 0 0 0 0 0 0 ds0 Φ(r̄0 ) n̂ · ∇g(r̄, r̄ 0 ) , (33) ds Φ(r̄ ) n̂ · ∇g(r̄, r̄ ) + ds Φ(r̄ ) n̂ · ∇g(r̄, r̄ ) = P V S piece S where P V denotes the principal value and ‘piece’ refers to the integration over the local domain shown in Fig. 3. For this integral, we use the local coordinates: ds0 = dX 0 , n̂ = Ẑ 0 , etc (34) 7 r̄ Z0 PSfrag replacements −a r̄ = (0, Z) a X0 Figure 3: Zoom on the rough surface. The integral becomes: Z Z 0 0 0 ds Φ(r̄ ) n̂ · ∇g(r̄, r̄ ) = lim lim Φ(r̄) a→0 |Z|→0 piece a dX 0 −a ∂ g(r̄, r̄0 ) . ∂Z 0 (35) Over the small piece, we have: |r̄ − r̄0 | = p X 02 + (Z − Z 0 )2 , (36a) 1 i (1) |r̄ − g = H0 (|r̄ − r̄ 0 |) → − ln( 2π 2 µ 4¶ Z ∂g . = ∂Z 0 Z 0 =0 2π(X 02 + Z 2 ) Thus, the integral becomes Z Z = lim lim Φ(r̄) piece a r̄0 | ), 1 Z 02 a→0 |Z|→0 2π X + Z 2 −a · ¸ 0 a Φ(r̄) −1 X = lim lim tan a→0 |Z|→0 2π Z −a   1 Φ(r̄) for Z > 0 = 2 − 1 Φ(r̄) for Z < 0 (36b) (36c) dX 0 (37) 2 The two parts of the integral then become: Z Φinc (r̄) + P V ds0 Φ(r̄0 ) n̂ · ∇g(r̄, r̄ 0 ) + S Z ds0 Φ(r̄0 ) n̂ · ∇g(r̄, r̄ 0 ) − Φinc (r̄) + P V S 1 Φ(r̄) = Φ(r̄) , 2 1 Φ(r̄) = 0 , 2 (38a) (38b) (38c) which is cast into one equation: Φinc (r̄) + P V Z 1 ds0 Φ(r̄0 ) n̂ · ∇g(r̄, r̄ 0 ) = Φ(r̄) . 2 S (39) Eq. (39) is called the magnetic field integral equation (since Φ represents the magnetic field) and is an integral equation of the second kind. 8 Section 5. Solving the EFIE In the same way as before, we can define the kernel as 0 K(x, x ) = s 1+ µ df dx ¶2 n̂ · ∇g(r̄, r̄ 0 )|z=f (x),z 0 =f (x0 ) , (40) and write the MFIE as Φinc (r̄) + P V 5 Z ∞ dx0 Φ(x0 ) K(x, x0 ) = −∞ 1 Φ(x) . 2 (41) Solving the EFIE Upon using the notation introduced in Section 4.1, the problem is reduced to solving the following integral equation: Z L/2 dx0 K(x, x0 ) U (x0 ) = b(x) . (42) −L/2 Let us subdivide the integration domain into N small elements, each of length ∆ = L/N , and centered at xm (m ∈ [1, N ]). Thus, constraining the observation at these discrete locations, the integral equation becomes Z L/2 dx0 K(xm , x0 ) U (x0 ) = b(xm ) . (43) −L/2 Next, if we suppose that U (x) is constant in each interval, we replace the integral by a sum over all segments, excluding the singular term: Z X ∆x K(xm , xn ) U (xn ) + U (xm ) dx0 K(xm , x0 ) = b(xm ) . (44) m n=1 n6=m Note that we have to single out the singularity, i.e. the mth interval because K(xm , x0 ) is singular at x0 = xm . This part is known as the self-patch contribution. For a 2D problem, we have: K(xm , x0 ) = i (1) p H (k (xm − x0 )2 + (f (xm ) − f (x0 ))2 ) . 4 0 (45) Upon approximating f (x0 ) − f (xm ) ' f 0 (xm )(x0 − xm ), we write: K(xm , x0 ) = (1) p i (1) H0 (k|x0 − xm | 1 + f 0 (xm )2 ) . 4 (46) For small argument: H0 (α) ≈ 1 + i π2 ln(αγ/2) where γ is the Euler constant (γ ' 1.78). Therefore: µ · ¶¸ p γ i 1 0 0 02 k(x − xm ) 1 + f (xm ) . K(xm , x ) = 1 + i ln (47) 4 π 2 9 Therefore, the integral becomes: Z 0 0 dx K(xm , x ) = 2 m Z xm +∆x/2 0 0 dx K(xm , x ) = 2 xm Z ∆x/2 dx0 K(xm , x0 + xm ) 0 · µ ¶¸ Z i ∆x/2 1 γ 0p 0 02 ' kx 1 + f (xm ) dx 1 + i ln 2 0 π 2 ½ µ ¶¾ p 2 γk i∆x 02 1 + i ln ∆x 1 + f (xm ) , = 4 π 4e (48) where ln(e) = 1. We can therefore cast the integral equation into a matrix equation of the form: N X Amn Un = bm , (49) n=1 where Un = U (xn ) is the unknown, bm = b(xm ) ,   ∆x K(xm , xn ) ½ µ ¶¾ Amn = i∆x p γk 2 02   4 1 + i π ln 4e ∆x 1 + f (xm ) The system can now easily be solved numerically. (50a) (50b) for m 6= n for m = n (50c)