A Matrix Basis Formulation For The Green's Functions Of Maxwell's Equations And The Elastic Wave Equations In Layered Media
aa r X i v : . [ m a t h - ph ] A ug A MATRIX BASIS FORMULATION FOR THE GREEN’SFUNCTIONS OF MAXWELL’S EQUATIONS AND THE ELASTICWAVE EQUATIONS IN LAYERED MEDIA
WENZHONG ZHANG ∗ , BO WANG † , AND
WEI CAI ‡ Abstract.
A matrix basis formulation is introduced to represent the 3 ×
1. Introduction.
The layered media dyadic Green’s functions (LMDG) for theMaxwell’s equations and the elastic wave equations are studied and used in the in-tegral equation solvers for wave fields [12, 9, 11]. These Green’s functions are 3 × × E z - H z formulation [3, 7], etc. The Sommerfeld potential and the transverse potential reducethe number of entries to be calculated to 5 while the E z - H z approach uses merely2 scalar terms, based on the TE/TM mode decomposition. For the elastic waveequation, the dyadic Green’s function for the half-space problem was discussed in [6].The purpose of this paper is to present a new matrix representation of the 3 × ∗ Department of Mathematics, Southern Methodist University, Dallas, TX 75275. † LCSM(MOE), School of Mathematics and Statistics, Hunan Normal University, Changsha, Hu-nan, 410081, P. R. China. Department of Mathematics, Southern Methodist University, Dallas,TX 75275. This author acknowledges the financial support provided by NSFC (grant 11771137),the Construct Program of the Key Discipline in Hunan Province and a Scientific Research Fund ofHunan Provincial Education Department (No. 16B154). ‡ Corresponding author, Department of Mathematics, Southern Methodist University, Dallas, TX75275( [email protected] ). 1 reen’s functions can be naturally decomposed into independent TE and TM com-ponents within this formulation, leading to the 2-term E z - H z result, but also witha clearer explanation about the interface treatment. Meanwhile, the elastic waveGreen’s function is decomposed into S-wave components and P-wave components bymatrix rows. Third, the rotational symmetry allows us to apply fast solvers easily,e.g. the fast multipole method in layered media [13]. We also develop a vector basisformulation which is simplified from the matrix version, used for the LMDG of themix-phase elastic wave equations where the source originates in fluid medium.The rest of this paper is organized as follows. In Section 2 we establish the theoriesof the matrix basis, propose the formulation, and provide guidelines about how thetheories and the formulation are applied to the LMDG of the Maxwell’s equations andthe elastic wave equations. In Section 3 the details of the Maxwell’s Green’s functionsin layered media are explained, including a 5-term matrix-based general formulationand the concise 2-term formulation. In Section 4 we discuss the layered elastic waveGreen’s functions. A 5-term formulation separating the S-wave and the P-wave bydefinition is proposed. Then, the result is generalized to the vector case for sourcesin fluid media. The mixed interface conditions between different medium phase typesare enumerated in details.
2. The matrix basis formulation.
In this section we set up the matrix basisused for the Green’s functions of the Maxwell’s equations or the elastic wave equation,and develop basic theories of the basis coefficients.Suppose in both problems the domain has a total of L + 1 layers, indexed by0 , · · · , L from top to bottom, separated by the interfaces z = d l , 0 ≤ l ≤ L − k i , 1 ≤ i ≤ I , consist all the wave numbers inthese layers. Note that there are distinct elastic wave numbers for the S-wave and theP-wave in each solid layer. The interaction between the fixed source r ′ = ( x ′ , y ′ , z ′ )and the target r = ( x, y, z ) is studied, assuming both of them do not locate on theinterfaces.Take the 2-D Fourier transform from ( x − x ′ , y − y ′ ) to ( k x , k y ):(2.1) f ( x, y ) = 14 π Z Z R e i k x ( x − x ′ )+i k y ( y − y ′ ) b f ( k x , k y ) dk x dk y . Let ( k ρ , α ) be the polar coordinates of ( k x , k y ). Let the field F be the field extensionfrom C with certain functions of k x and k y (2.2) F = C (cid:16) i q k i − k ρ , e i √ k i − k ρ d j , e i √ k i − k ρ z , e i √ k i − k ρ z ′ : 1 ≤ i ≤ I, ≤ j ≤ L − (cid:17) , where other variables are treated as given constants. Obviously, functions in F donot depend on α . Note that k ρ ∈ F because k ρ = k i + (cid:16) i q k i − k ρ (cid:17) . Then, define the field F by the 2-term extension(2.3) F = F (i k x , i k y ) . Remark z and z ′ are indeed redundant in the definition of F for the matrix basis theory. They are included only for the convenience of statementsin the following sections. ne of our anticipations on the matrix basis is to represent the tensor Green’sfunctions with all the coefficients in F , i.e. the information of the polar angle is onlykept in the matrix basis. For this purpose, based on the observation on calculatedformulas [10], we propose the matrix basis J , · · · , J in the frequency domain asfollows. Proposition
These matrices form a basis of F × : J = , J = , J = k x k y , J = k x i k y , J = − k x − k x k y − k x k y − k y
00 0 0 , J = − i k y i k x , J = k y − i k x , J = k x k y k y − k x − k x k y
00 0 0 , J = − . (2.4)The proof is trivial.The matrix basis deserves to be divided into two groups once the products betweenthem are studied. For any subfield K ⊂ F , define vector spaces R ( K ) = span K ( J , · · · , J ) , I ( K ) = span K ( J , · · · , J ) , M ( K ) = span K ( J , · · · , J , J , · · · , J ) . (2.5) Proposition
Let K be any subfield of F containing k ρ . Then, • M ( K ) = R ( K ) ⊕ I ( K ) is the direct sum. • R ( K ) , M ( K ) are rings with matrix addition and matrix multiplication.Proof. The direct sum is obvious. For the ring property, notice the identity matrix(2.6) I = J + J ∈ R ( K ) ⊂ M ( K ) , and the product table of the matrices J , · · · , J (cid:2) J T · · · J T (cid:3) T · (cid:2) J · · · J (cid:3) = J J J − k ρ J − k ρ J − k ρ J J − k ρ J − k ρ J − k ρ J J k ρ J − k ρ J − J − J k ρ J + J k ρ J − k ρ J − k ρ J − J J − J − J J − J (2.7)which ensures the matrix multiplication is closed in either M ( K ) or R ( K ).The product table (2.7) immediately implies the product rules below. roposition Let K be any subfield of F containing k ρ . • If A ∈ R ( K ) , B ∈ R ( K ) , then A · B ∈ R ( K ) . • If A ∈ R ( K ) , B ∈ I ( K ) , then A · B ∈ I ( K ) . • If A ∈ I ( K ) , B ∈ R ( K ) , then A · B ∈ I ( K ) . • If A ∈ I ( K ) , B ∈ I ( K ) , then A · B ∈ R ( K ) . The behavior resembles the real numbers and the imaginary numbers, which is whythe letters R and I are used here. Definition
Define the linear space (2.8) R = span F { J , · · · , J } = X j =1 a j J j : a j ∈ F . The linear expansion of functions in R with basis J , · · · , J by definition is calledthe matrix basis formulation. In future sections we will claim that the matrix basis formulation can be usedto efficiently solve the Green’s functions for the Maxwell’s equations and the elasticwave equation in each layer. Since the tensors are across the layers and to be solvedtogether, theories of the block matrices are also necessary.For any subfield K ⊂ F and any p, q ∈ N , define the linear spaces of block matrices M p × q ( K ) = X j =1 K j ⊗ J j : K j ∈ K p × q , ≤ j ≤ , R p × q ( K ) = X j =1 K j ⊗ J j : K j ∈ K p × q , ≤ j ≤ , I p × q ( K ) = X j =6 K j ⊗ J j : K j ∈ K p × q , ≤ j ≤ , (2.9)where ⊗ is the Kronecker product. It is straightforward that in the above definitionsthe decompositions are unique (since in each 3 × M p × q ( K ) = R p × q ( K ) ⊕ I p × q ( K )is the direct sum. Moreover, the product rules are easily generalized to block matrices. Proposition
Let p, q, r ∈ N . Let K be any subfield of F containing k ρ . • If ¯ A ∈ R p × r ( K ) , ¯ B ∈ R r × q ( K ) , then ¯ A · ¯ B ∈ R p × q ( K ) . • If ¯ A ∈ R p × r ( K ) , ¯ B ∈ I r × q ( K ) , then ¯ A · ¯ B ∈ I p × q ( K ) . • If ¯ A ∈ I p × r ( K ) , ¯ B ∈ R r × q ( K ) , then ¯ A · ¯ B ∈ I p × q ( K ) . • If ¯ A ∈ I p × r ( K ) , ¯ B ∈ I r × q ( K ) , then ¯ A · ¯ B ∈ R p × q ( K ) . The following theorem will lead to the main result of this section.
Theorem
Suppose p, q, r ∈ N , the block matrices ¯ A ∈ R p × r ( F ) , ¯ X ∈ M r × q ( F ) and ¯ B ∈ R p × q ( F ) satisfy ¯ A · ¯ X = ¯ B . Then, there existsa “filtered” block matrix ¯ X ∈ R r × q ( F ) , i.e. each block of ¯ X has the matrix basisrepresentation, so that ¯ A · ¯ X = ¯ B . roof. We filter the solution from M r × q ( F ) to R r × q ( F ) with an intermediatestep in R r × q ( F ).First, write the direct sum decomposition ¯ X = ¯ X ⊕ ¯ X , where ¯ X ∈ R r × q ( F )and ¯ X ∈ I r × q ( F ). By Proposition 2.6 we immediately get ¯ A · ¯ X ∈ R p × q ( F ) and¯ A · ¯ X ∈ I p × q ( F ), so(2.11) ¯ A · ¯ X + ¯ A · ¯ X is the direct sum decomposition of ¯ B ∈ M p × q ( F ). Therefore ¯ A · ¯ X = ¯ B .Then, let ¯ A = X j =1 A j ⊗ J j , ¯ X = X j =1 X j ⊗ J j , ¯ B = X j =1 B j ⊗ J j where each A j ∈ F p × r , X j ∈ F r × q and B j ∈ F p × q . When treating X j as the solutionto the linear equation ¯ A · ¯ X = ¯ B , the equation is equivalent to(2.12) X u =1 5 X v =1 ( A u X v ) ⊗ ( J u J v ) = X j =1 B j ⊗ J j which, by the product table (2.7), is indeed equivalent to the linear system e A e X = e B ,where a stacked form of ¯ X is used(2.13) e A = A − k ρ A A − k ρ A − k ρ A A A − k ρ A , e X = X ... X , e B = B ... B . Let a ≤ min(5 p, r ) be the rank of e A ∈ F p × r . The diagonalization of e A implies theexistence of full-rank matrices S ∈ F p × p and T ∈ F r × r such that(2.14) e A = S · (cid:20) I a a × (5 r − a ) (5 p − a ) × a (5 p − a ) × (5 r − a ) (cid:21) · T , so the entries lower than the a -th row of S − e B are all zero, and(2.15) e X = T − · "(cid:2) I a a × (5 p − a ) (cid:3) · (cid:16) S − e B (cid:17) (5 r − a ) × q ∈ F r × q will satisfy e A e X = e B . Write e X in the stacked form(2.16) ˜ X = X ... X where each X j ∈ F r × q , the matrix(2.17) ¯ X = X j =1 X j ⊗ J j ∈ R r × q ( F )is as desired. he coming discussion in later sections on the LMDG of the Maxwell’s equationsand the elastic wave equation are generally presented in the following pattern: therestricting equations of the problem, including linear equations derived from the inter-face conditions and the radiation conditions, are re-formatted using the matrix basis J , · · · , J . The solution filtering theorem then works on the linear system of the ten-sors in the layers, so that filtered solutions are proven available using the matrix basisformulation of R . Finally, the formulation helps simplify the restricting equations,leaving only the basis coefficients to be solved.
3. Application to the Maxwell’s equations in layered media.
In thissection we first give a brief introduction about the dyadic Green’s functions of thetime-harmonic Maxwell’s equations in the free space, then we will discuss the Green’sfunctions in layered media, and the simplification using the matrix basis formulation(2.8).
Suppose in the free space the medium has constant permittivity ε andconstant permeability µ . Assuming the time dependence is harmonic, i.e. in termsof exp(i ωt ), the source-free Maxwell’s equations in the free space is simplified in the ω -domain of the Fourier transform as the equations of the electric displacement flux ~D ( r ), the electric field ~E ( r ), the magnetic flux density ~B ( r ) and the magnetic field ~H ( r ), namely, ~D = ε ~E~B = µ ~H ∇ × ~E = − i ωµ ~H ∇ × ~H = i ωε ~E ∇ · ~D = 0 ∇ · ~B = 0(3.1)and the wave number is given by(3.2) k = p ω εµ. When dealing with these equations, the Lorenz gauge condition is often introduced,which allows us to use a vector potential ~A ( r ) to represent the electric field ~E andthe magnetic field ~H as(3.3) ~E = − i ω (cid:18) I + ∇∇ k (cid:19) ~A, ~H = 1 µ ∇ × ~A, and the flux vectors ~D and ~B are replaced using ε ~E and µ ~H , respectively. The vectorpotential satisfies the Helmholtz equation(3.4) ∇ ~A + k ~A = ~ . The choice of the vector potential is not unique. Indeed, for any function φ ∈ C ( R )satisfying the Helmholtz equation ∇ φ + k φ = 0, ~A + ∇ φ can be used to replace ~A in the above identities. he dyadic Green’s functions for the free space Maxwell’s equations are definedusing a 3 × G A ( r ; r ′ ) such that the electric field dyadic Green’sfunction G E ( r ; r ′ ) and the magnetic field dyadic Green’s function G H ( r ; r ′ ) are rep-resented by(3.5) G E = − i ω (cid:18) I + ∇∇ k (cid:19) G A , G H = 1 µ ∇ × G A . The potential tensor satisfies the Helmholtz equation(3.6) ∇ G A + k G A = 1i ω δ ( r − r ′ ) I where I is the 3 × r →∞ r (cid:18) ∂∂r − i k (cid:19) G ( r ; r ′ ) = for G = G E and G = G H , here r = | r | .For the same reason, the tensor potential is not unique. A commonly used solutionto (3.6) is given by(3.8) G f A ( r ; r ′ ) = − ω e i k | r − r ′ | π | r − r ′ | I = − ω g f ( r ; r ′ ) I , where g f ( r ; r ′ ) is the free space Green’s function of the Helmholtz equation. G f A alsosatisfies the Sommerfeld radiation condition. Now suppose the space is horizontally stratified as layers 0 , · · · , L arrangedfrom top to bottom, separated by planes z = d , · · · , z = d L − where d > · · · >d L − , and each layer is homogeneous with constant permittivity ε j and constantpermeability µ j , j = 0 , · · · , L , respectively. We append the layer index to the endof the subscript of any layer-dependent variable or function to represent its valuespecified in that layer, e.g. the wave number is then(3.9) k j = q ω ε j µ j , j = 0 , · · · , L. in layer j . For simplicity, the layer index is sometimes omitted, and one may assumethe variable is a piecewise function of z . These subscript rules are applied to the restof this paper, including the later section about the elastic wave equations. The time-harmonic Maxwell’sequations in the interior of each layer has the same form as in (3.1), while the followinginterface conditions must be satisfied [9]: between the adjacent layers,(3.10) J n × ~E K = ~ , J n · ~D K = 0 , J n × ~H K = ~ , J n · ~B K = 0 . J · K represents the jump of the value at the interface, i.e. across the interface z = d ,(3.11) J f K = lim z → d + f − lim z → d − f. he dyadic Green’s functions are again given using the tensor potential G A as(3.12) G E = − i ω (cid:18) I + ∇∇ k (cid:19) G A , G H = 1 µ ∇ × G A . The tensor potential G A satisfies the Helmholtz equation(3.13) ∇ G A + k G A = 1i ω δ ( r − r ′ ) I , while the interface conditions are(3.14) J n × G E K = , J ε n · G E K = ~ T , J n × G H K = , J µ n · G H K = ~ T . In horizontally layered media n = e = (cid:2) (cid:3) T . In addition, the Green’s functionsmust satisfy the Sommerfeld radiation conditions (3.7). Take the the 2-D Fouriertransform (2.1) from ( x − x ′ , y − y ′ ) to ( k x , k y ). Suppose r = ( x, y, z ) locates in layer t , and r ′ = ( x ′ , y ′ , z ′ ) locates in layer j . The layer index is default to the target layer t when not specified. Notations of the polar coordinate pair ( k ρ , α ) are kept.We begin with the separation of the z variable from the tensor potential b G A ,which will lead to the reaction field decomposition.For the gradient alternative in the frequency domain, define the notation b ∇ by(3.15) b ∇ = i k x i k y ∂ z when a function of z follows it. The b ∇ b ∇ , b ∇ now refer to b ∇ b ∇ T and b ∇ T b ∇ , respectively.Recall that the right-hand side of the Helmholtz equation (3.13) is nontrivial ifand only if r ′ is in the same layer as r , i.e. j = t , define(3.16) b G r A ( r ; r ′ ) = b G A ( r ; r ′ ) − δ j,t b G f A ( r ; r ′ ) , where δ j,t is the Kronecker delta function. The complementary part b G r A is called thereaction field, and satisfies the homogeneous Helmholtz equation(3.17) b ∇ b G r A + k b G r A = , i.e. ∂ zz b G r A + ( k − k ρ ) b G r A = . Define(3.18) k z = q k − k ρ where the square root takes nonnegative real part. The general solutions to (3.17),when treated as an ordinary differential equation of z , is given by(3.19) b G r A = e i k z z b G r ↑ A + e − i k z z b G r ↓ A where b G ↑ A and b G ↓ A are piecewise constants with respect to z , namely,(3.20) b G r ↑ A = b G r ↑ A,t , b G r ↓ A = b G r ↓ A,t hen in layer t . Indeed, we can also write b G f A in the frequency domain in a similarform b G f A = − ωk z,j e i k z,j | z − z ′ | I = − ωk z,j e i k z,j ( z − z ′ ) { z>z ′ } I + − ωk z,j e i k z,j ( z ′ − z ) { z
2, the potential tensor b G A + ∂ z b f J + b f J + ∂ z b f J + b f J can be used to replace b G A . To eliminate thedegrees of freedom in the coefficients, define the functions b , b and b by lineartransforms of a l : b = a ,b = 1 µ ( a − ∂ z a ) ,b = 1 µ (cid:0) ∂ z a + k ρ a − k ρ ∂ z a (cid:1) , (3.40)so that b G E and b G H in (3.12) can be represented by b G E = − i ωk (cid:18) k b J + µk ρ b J + µ∂ z b J + µb J + (cid:18) k k ρ b + µk ρ ∂ z b (cid:19) J (cid:19) , b G H = 1 µ (cid:18) b J + µb J + (cid:18) k ρ ∂ z b − µk ρ b (cid:19) J − ∂ z b J (cid:19) . (3.41)Each b l has the reaction component decomposition corresponding to (3.37)(3.42) b l = b l ( k ρ , z, z ′ ) = δ j,t b f l ( k ρ , z, z ′ ) + e i k z z b r ↑ l ( k ρ , z ′ ) + e − i k z z b r ↓ l ( k ρ , z ′ ) , where b f1 = a f1 ,b f2 = 1 µ (cid:0) a f2 − ∂ z a f3 (cid:1) ,b f3 = 1 µ (cid:0) ∂ z a f1 + k ρ a f4 − k ρ ∂ z a f5 (cid:1) ,b r ∗ = a r ∗ ,b r ∗ = 1 µ ( a r ∗ − τ ∗ i k z a r ∗ ) ,b r ∗ = 1 µ (cid:0) τ ∗ i k z a r ∗ + k ρ a r ∗ − k ρ τ ∗ i k z a r ∗ (cid:1) , (3.43)due to (3.37), here ∗ ∈ {↑ , ↓} , τ ↑ = 1, τ ↓ = −
1. Specifically, since we have chosen G A = G f A = − g f / (i ω ) I , it’s clear that(3.44) a f1 = a f2 = − ω b g f , a f3 = a f4 = a f5 = 0 , where b g f = i e i k z,j | z − z ′ | / (2 k z,j ). Therefore b f1 = − ω b g f ,b f2 = − ω µ j b g f ,b f3 = − ω µ j ∂ z b g f = − ∂ z ′ b f2 . (3.45) ach b l also satisfies the Helmholtz equation(3.46) ∂ zz b j + k z b j = 0piecewisely in each layer provided z ′ = z .For solving b , b and b we take a review of the interface equations and theradiation equations. One can easily verify the interface equations (3.14), which werereinterpreted in the frequency domain as in (3.28), are equivalent to the following bycomparing the matrix basis coefficients. For example, from J · b G E = − i ω (cid:0) b J + ω − ε − ∂ z b J + ( k − ρ b + ω − ε − k − ρ ∂ z b ) J (cid:1) the continuity equations J − i ωb K = 0 , J − i ω − ε − ∂ z b K = 0 , J − i ωk − ρ b − i ω − k − ρ ε − K = 0are revealed. A complete list by items is given below: J n × b G E K = ⇔ J J · b G E K = ⇔ J b K = 0 , s ε ∂ z b { = 0 , s ε ∂ z b { = 0; J ε n · b G E K = ~ ⇔ J J · ε b G E K = ⇔ J b K = 0 , J b K = 0; J n × b G H K = ⇔ J J · b G H K = ⇔ J b K = 0 , J b K = 0 , s µ ∂ z b { = 0; J µ n · b G H K = ~ ⇔ J J · µ b G H K = ⇔ J b K = 0 . (3.47)The radiation equations (3.31) and (3.34) are reduced to(3.48) b r ↓ l, = 0 , b r ↑ l,L = 0in the top and the bottom layer, respectively, i.e. waves coming from z = ±∞ areprohibited in the reaction field decomposition. For each l , the above are a total of2 L + 2 linear equations of b r ↑ l and b r ↓ l from L + 1 layers. These linear equations aresolvable, from the knowledge of the acoustic wave equation in layered media: • − i ωb is exactly the reflection/transmission coefficient in the frequency do-main of the Green’s function of the Helmholtz equation in layered media,with piecewise constant material parameters 1 /ε . Thus we can solve b in thefrequency domain like solving the known scalar layered Helmholtz problem[14]. • Similarly, − i ωµ j b is exactly the one with piecewise constant parameters 1 /µ . • The linear system regarding b r ∗ has exactly the same coefficients as b r ∗ for theunknowns, so it’s solvable since b r ∗ are uniquely determined by the physicalproblem. Moreover,(3.49) − ∂ z ′ b = − ∂ z ′ δ j,t b f2 − e i k z z ∂ z ′ b f ↑ − e − i k z z ∂ z ′ b f ↓ satisfies every equation that b should satisfy, so by uniqueness,(3.50) b = − ∂ z ′ b , i.e. b r ∗ = − ∂ z ′ b r ∗ . The b and b functions are corresponding to the TE mode component and the TMmode component in the E z - H z formulation [3, 7], respectively. emark b G E and b G H we don’t need the intermediate, un-determined tensor potential b G A anymore. This paper derived the above formulationvia b G A because some previous work did need its formulation, such as in the integralequation applications in [5]. Remark
Remark b G r ∗ A satisfyingthe above interface equations and radiation equations for certain values of k ρ , witheach b G f A replaced by 0. It is corresponding to a pole in the frequency domain [15].In such situation we can still derive the simplified formulation using terms b , b and b , but b plays an independent role and is not anymore tied with b . Here wetake a quick review on the transverse potential and the Sommerfeld potential for-mulations and show how to reach them from the matrix basis formulation. Bothformulations restrict certain 5 entries of the 3 × b G A to be nonzero, whichuniquely determines the tensor potential. Here we claim the potential tensors in theseformulations have the matrix basis representation, and can be derived using b and b . Due to the uniqueness of b and b , it suffices to explicitly construct them.The transverse potential takes the form(3.51) b G t A = × ×× × × , where each × marks a nonzero entry. We claim b G t A = a J + a J + a J . By (3.40),(3.52) b = a , b = 1 µ a , b = − ∂ z ′ b = 1 µ ( ∂ z a − k ρ ∂ z a ) . Since a l , b l satisfy the Helmholtz equation (3.39) and (3.46), respectively, we have(3.53) a = b , a = µb , a = b − µ∂ z ∂ z ′ b k ρ k z . The Sommerfeld potential takes the form(3.54) G S A = × ×× × × and we claim G S A = a J + a J + a J . Again by (3.40),(3.55) b = a , b = 1 µ a , b = − ∂ z ′ b = 1 µ ( ∂ z a + k ρ a ) , so(3.56) a = b , a = µb , a = − µ∂ z ′ b + ∂ z b k ρ . emark b G t A , although the coefficient a has a k ρ factor in the denominator, there’s no singularity in the integrand of b G t A at k ρ = 0since they can be cancelled out with the entries of J . The same does happen tothe Sommerfeld potential b G S A , but it’s not explicitly shown in the expression of a J .Numerically we should take some care if the values as k ρ →
4. Application to the elastic wave equation in layered media.
In thissection we will apply the matrix based formulation to the dyadic Green’s functionof the elastic wave equation in layered media. The interface conditions for variouscontacting media phase types will be discussed. In the end we also give a briefdiscussion on the case when source is inside zero-viscosity fluid layer, with a simplifiedvector basis formulation. In these problems, we suppose the displacement of the mediais small, so that linearization of the elastic wave equations is in general applicable.
Suppose the homogeneous and isotropic material occupies the entirespace, with density ρ and Lam´e constants λ, µ . Define(4.1) γ = λ + 2 µ. If the material is solid, both λ, µ >
0. If the material is (compressible) liquid or gas,we assume the viscosity is negligible (also for the rest of this paper; otherwise weshould treat it in the same way as a solid layer), then µ = 0, and γ = λ . Note thatthe shear modulus µ is different from the dynamic viscosity in fluid. Suppose the timedependence is harmonic, i.e. in terms of exp(i ωt ).With external force b ( r ), the elastic wave equation in solid is a partial differentialequation of the displacement u ( r )(4.2) − ω ρ u = ∇ · T + b assuming | u | ≪
1, where the 3 × T is defined as(4.3) T ij = λδ i,j X l =1 ∂u l ∂x l + µ (cid:18) ∂u i ∂x j + ∂u j ∂x i (cid:19) , where ( x , x , x ) is used as an alternative notation for ( x, y, z ), and u = ( u , u , u )[8]. The equivalent form not using the stress tensor T is(4.4) ( λ + µ ) ∇∇ · u + µ ∇ u + ω ρ u = − b . Given source location r ′ , the dyadic Green’s function G ( r ; r ′ ) is a 3 × λ + µ ) ∇∇ · G + µ ∇ G + ω ρ G = − δ ( r − r ′ ) I . The solution is known as(4.6) G = G f ( r ; r ′ ) = 1 µ (cid:18) I + ∇∇ k s (cid:19) g s ( r ; r ′ ) − γ ∇∇ k c g c ( r ; r ′ ) , where(4.7) k s = s ω ρµ , k c = s ω ργ re the wave numbers of the S-wave and the P-wave, respectively, and(4.8) g s ( r ; r ′ ) = e i k s | r − r ′ | π | r − r ′ | , g c ( r ; r ′ ) = e i k c | r − r ′ | π | r − r ′ | are the free space Green’s functions of the Helmholtz equation with wave numbers k s and k c , respectively [8].In fluid, i.e. in liquid or gas, assuming the external force b is conservative (forthe rest of this paper as well), there’s only the P-wave propagating in the media, andthe acoustic wave equation with respect to the displacement u is λ ∇∇ · u + ω ρ u = − b , (4.9) ∇ × u = , (4.10)where the first equation repeats (4.4), and the second one is introduced because ofzero viscosity. Also consider the linearized Navier–Stokes equation(4.11) − ω ρ u = ρ (cid:18) D v Dt (cid:19) ∧ = −∇ p + b where v is the velocity, D/Dt is the material derivative and ∧ represents the Fouriertransform of time. For the pressure p we get(4.12) p = − λ ∇ · u + p where p is a constant. For the sake of convenience let p = 0. Take the divergenceof the first equation of (4.9),(4.13) ∇ p + ω ρλ p = 1 λ ∇ · b . The dyadic Green’s function of the P-wave propagation is often proposed in terms ofthe pressure p , namely(4.14) ∇ g p + ω ρλ g p = − δ ( r − r ′ )with solution g p = g f p ( r ; r ′ ) = g c ( r ; r ′ ). For the displacement u , the correspondingGreen’s function g f u satisfies(4.15) g f u = 1 ω ρ ∇ g f p which straightforwardly follows (4.9). In layered media, againwe suppose the space is horizontally stratified as layers 0 , · · · , L arranged from top tobottom, separated by planes z = d , · · · , z = d L − where d > · · · > d L − , and themedium in each layer is homogeneous with density ρ l and Lam´e constants λ l and µ l , l = 0 , · · · , L , respectively. Define γ = λ + 2 µ , and the wave numbers(4.16) k s = s ω ρµ , k c = s ω ργ are the same as in the free space. The notation for variables in different layers followsthe convention introduced at the beginning of Subsection 3.2. In fluid media, k s isnot defined because µ = 0 and S-wave does not propagate. .2.1. The interface conditions of the elastic wave equation. Before en-tering the discussion on the Green’s functions, we take a detailed study on the interfaceconditions. Recall the equation (4.2) in terms of the stress tensor T , if the externalforce b and the displacement u are not singular near the interface, then by divergencetheorem on a flat cylinder crossing the boundary with infinitesimal thickness, we getthe interface condition(4.17) J n · T K = ~ . With n = e in our problem, it suffices to consider the entries T = µ (cid:18) ∂u ∂x + ∂u ∂x (cid:19) , T = µ (cid:18) ∂u ∂x + ∂u ∂x (cid:19) , T = γ ∂u ∂x + λ (cid:18) ∂u ∂x + ∂u ∂x (cid:19) (4.18)for the continuity equations. The identity (4.17) also implies some additional regularconditions to be listed below [8]. • Across the solid-solid interface, ∂u /∂x must be regular, so J u K = 0. Thenfrom the continuity of T we further get J u K = 0, and similarly J u K = 0. • Across the solid-liquid/gas interface, we also have J u K = 0. Since in the fluidside µ = 0 so that T = T = 0, no additional condition is required. • Across the fluid-fluid interface, the conditions on T and T are no morenecessary, therefore only J T K = 0 and J u K = 0 are required. Note that inthis case T = λ ∇ · u = p , and u = − ( ω ρ ) − ∂p/∂x . • In the half-space problem where the media has an interface against the vac-uum where no acoustic wave propagates and λ = µ = 0, the zero-tractionconditions J T l K = 0 are required, l = 1 , ,
3, while the displacement on theboundary is set free.In summary, the interface equations by case are listed below: • Across the solid-solid interface, J T l K = 0 , J u l K = 0 , l = 1 , , . (4.19) • Across the solid-fluid interface,(4.20) J u K = 0 , J T l K = 0 , l = 1 , , , where from the fluid side T = T = 0 automatically holds. • Across the fluid-fluid interface,(4.21) J u K = 0 , J T K = 0 , or equivalently, in terms of the pressure p ,(4.22) s ρ ∂p∂z { = 0 , J p K = 0 . • On the solid-vacuum interface,(4.23) J T l K = 0 , l = 1 , , . On the fluid-vacuum interface,(4.24) J T K = 0 . To flexibly handle various interface circumstances, we separate these equations intofour groups: (a) J T K = 0; (b) J u K = 0; (c) J T K = J T K = 0 ; (d) J u K = J u K = 0.For each of the above cases we simply join them in need.When the source is in solid medium, the Green’s function in terms of the dis-placement G is a 3 × We begin with the case when all the layers are solid. Then the case with existing fluidor vacuum layers will follow. Suppose the source r ′ = ( x ′ , y ′ , z ′ ) locates in layer j , andthe target r = ( x, y, z ) is in layer t . Take the 2-D Fourier transform ( x − x ′ , y − y ′ )to ( k x , k y ) as defined in (2.1). Suppose ( k ρ , α ) are the polar coordinates of ( k x , k y ).Like in (3.16) for the Maxwell’s equations, define the reaction field Green’s function(4.25) G r = G − δ j,t G f , then G r satisfies the homogeneous elastic wave equation(4.26) ( λ + µ ) ∇∇ · G r + µ ∇ G r + ω ρ G r = within each layer, which, in the frequency domain, is an ordinary differential equation(4.27) ( µ J + γ J ) ∂ zz b G r +( λ + µ )( J + J ) ∂ z b G r + (cid:0) µk sz ( J + J ) + ( λ + µ ) J (cid:1) b G r = of z , here(4.28) k sz = q k s − k ρ , k cz = q k c − k ρ . The above differential equation (4.27) has the general solution b G r = X ∗∈{↑ , ↓} (cid:16) ( − τ ∗ i k sz J + J ) e τ ∗ i k sz z + ( τ ∗ i k cz J + J ) e τ ∗ i k cz z (cid:17) X r ∗ (4.29)where τ ↑ = +1 , τ ↓ = − X r ∗ = X r ∗ t is piecewise constant (with respect to the variable z ) in each layer. The free space part b G f given by (4.6) takes the following form inthe frequency domain: b G f = 1 ω ρ j ( k s,j J + k ρ J + i k sz,j τ J + i k sz,j τ J + J ) b g s − ω ρ j ( − k cz,j J + i k cz,j τ J + i k cz,j τ J + J ) b g c (4.30)for z = z ′ , where τ is the sign function of z − z ′ ,(4.31) b g s = i e i k sz,j | z − z ′ | k sz,j , b g c = i e i k cz,j | z − z ′ | k cz,j . b G f actually has the same form as in (4.29): b G f = X ∗∈{↑ , ↓} (cid:16) ( − τ ∗ i k sz,j J + J ) e τ ∗ i k sz,j z + ( τ ∗ i k cz,j J + J ) e τ ∗ i k cz,j z (cid:17) X f ∗ , (4.32) here X f ↑ = 1 { z>z ′ } ( x f ↑ J + · · · + x f ↑ J ) , X f ↓ = 1 { z 5. Conclusion. In this paper, a matrix basis formulation is proposed for han-dling the dyadic Green’s functions of the Maxwell’s equations and of the elastic waveequation in layered media. The formulation is then used to further simplify the rep-resentation and derivation of the Green’s functions in both cases. As a corollary, avector basis formulation is introduced to handle the dyadic Green’s function of theelastic wave equation with source in fluid media. REFERENCES[1] A. Sommerfeld, Partial Differential Equations in Physics, Academic Press , New York, 1964.[2] A. Erteza and B. K. Park, Nonuniqueness of resolution of Hertz vector in presence of a boundary,and the horizontal problem, IEEE Trans. Antennas Propag. , vol. 17, no. 3, pp. 376–378,1969.[3] J. A. Kong, Electromagnetic field due to dipole antennas over stratified anisotropic media, GeoPhys. vol. 38, pp. 985–996, 1972.[4] K. A. Michalski, On the scalar potential of a point charge associated with a time-harmonicdipole in a layered medium, IEEE Trans. Antennas Propag. , vol. 35, no. 11, pp. 1299–1301, 1987.[5] K. A. Michalski and D. Zheng, Electromagnetic scattering and radiation by surfaces of arbitraryshape in layered media, Part I: Theory, IEEE Trans. Antennas Propag. , vol. 38, pp. 335–344, 1990.[6] T. J. Gosling and J. R. Willis, A line-integral representation for the stresses due to an arbitrarydislocation in an isotropic half-space, J. Mech. Phys. Solids , vol. 42, no. 8, pp. 1199–1221,1994.[7] W. C. Chew, J. L. Xiong and M. A. Saville, A Matrix-Friendly Formulation of Layered MediumGreen’s Function, IEEE Antennas Wirel. Propag. Lett. , vol. 5, pp. 490–494, 2006.[8] W. C. Chew, M. S. Tong and B. Hu, Integral Equation Methods for Electromagnetic and ElasticWaves, Synthesis Lectures on Computational Electromagnetics , vol. 3, no. 1, 2008.[9] W. Cai, Computational Methods for Electromagnetic Phenomena: Electrostatic in Solvation,Scattering, and Electron Transport, Cambridge University Press , 2013.[10] M. H. Cho and W. Cai, Efficient and Accurate Computation of Electric Field Dyadic GreensFunction in Layered Media, J. Sci. Comput. vol. 71, pp. 1319–1350, 2017.[11] D. Chen, M. H. Cho and W. Cai, Accurate and efficient Nystrom volume integral equationmethod for electromagnetic scattering of 3-D metamaterials in layered media, SIAM J.Sci. Comput. vol. 40, no. 1, pp. B259–B282, 2018.[12] W. C. Chew, Waves and Fields in Inhomogenous Media, Wiley-IEEE, Press (February 2, 1999).[13] B. Wang, W. Zhang and W. Cai, Fast Multipole Method for 3-D Helmholtz Equation in LayeredMedia, SIAM J. Sci. Comput. vol. 41, no. 6, pp. A3954–A3981, 2019.2714] B. Wang, D. Chen, B. Zhang, W. Zhang, M. H. Cho and W. Cai, Taylor expansion based fastmultipole method for 3-D Helmholtz equations in layered media, J. Comp. Phys. vol. 401,109008, 2020.[15] W. Zhang, B. Wang and W. Cai, Exponential Convergence for Multipole and Local Expansionsand Their Translations for Sources in Layered Media: Two-Dimensional Acoustic Wave,