Electromagnetic scattering by discrete random media. II: The coherent field

Journal of Quantitative Spectroscopy & Radiative Transfer 230 (2019) 86–105

Contents lists available at ScienceDirect

Journal of Quantitative Spectroscopy & Radiative Transfer journal homepage: www.elsevier.com/locate/jqsrt

Review

Electromagnetic scattering by discrete random media. II: The coherent field Adrian Doicu a,∗, Michael I. Mishchenko b a b

Deutsches Zentrum für Luft- und Raumfahrt (DLR), Institut für Methodik der Fernerkundung (IMF), Oberpfaffenhofen 82234, Germany NASA Goddard Institute for Space Studies, 2880 Broadway, New York, NY 10025, USA

a r t i c l e

i n f o

Article history: Received 22 January 2019 Revised 14 March 2019 Accepted 18 March 2019 Available online 19 March 2019

a b s t r a c t The computation of the coherent field in the case of a plane electromagnetic wave obliquely incident on a discrete random layer with non-scattering boundaries is addressed. For dense media, the analysis is based on a special-form solution for the conditional configuration-averaged exciting field coefficients, and is restricted to the computation of the so-called zeroth-order fields without a special treatment of the boundary regions. In this setting, we calculate the coherent fields reflected and transmitted by the layer, and the coherent field inside the layer. We found that these fields are analytically equivalent to plane electromagnetic waves, and investigated the fulfillment of the boundary conditions for the electric fields at the layer interfaces. The results are then particularized to the cases of normal incidence and a semi-infinite discrete random medium. For sparsely distributed particles, we present a self-consistent derivation of the coherent field and discuss the Twersky and Foldy approximations. © 2019 Elsevier Ltd. All rights reserved.

Contents 1. 2. 3.

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Recapitulation of the results of Ref. [1]. . . . . . . . . . . . . . . . . . . . . . . . Coherent field in dense media . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1. The coherent field reflected by the layer . . . . . . . . . . . . . . . . . 3.2. The coherent field transmitted by the layer . . . . . . . . . . . . . . 3.3. The coherent field inside the layer. . . . . . . . . . . . . . . . . . . . . . 3.3.1. The method of Fikioris and Waterman . . . . . . . . . . . . 3.3.2. The method of Tsang and Kong. . . . . . . . . . . . . . . . . . 3.3.3. Sparse-medium approximation . . . . . . . . . . . . . . . . . . 3.4. Normal incidence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5. Semi-infinite discrete random medium . . . . . . . . . . . . . . . . . . 4. Coherent field in a sparse medium. . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1. Conditional configuration-averaged exciting field coefficients 4.2. Coherent field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3. The Twersky and Foldy approximations. . . . . . . . . . . . . . . . . . 5. Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Appendix A. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

∗

Corresponding author. E-mail address: [email protected] (A. Doicu).

https://doi.org/10.1016/j.jqsrt.2019.03.011 0022-4073/© 2019 Elsevier Ltd. All rights reserved.

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . .

87 87 89 90 91 94 94 96 97 99 100 100 101 102 103 103 104 104 104

A. Doicu and M.I. Mishchenko / Journal of Quantitative Spectroscopy & Radiative Transfer 230 (2019) 86–105

1. Introduction

u = u0 +



A 0i u0i +

i

In the first part of this series [1], we derived the dispersion equation and gave a prescription for computing the configurationaveraged exciting field in the case of a plane electromagnetic wave obliquely incident on a discrete random layer with non-scattering boundaries. Using these results as a starting point we continue our survey with the analysis of the coherent field inside and outside the scattering medium. To compute these fields in the case of dense media, we will use the techniques developed by Waterman and Truell [2], Fikioris and Waterman [3,4], and Tsang and Kong [5,6] for calculating the coherent field in a semi-infinite discrete random medium. For sparse media, we will use a somewhat different technique. To reveal the difference between these two approaches, we consider the following Foldy–Lax model for the scattering by N particles randomly distributed throughout a volume V [1]:

u = u0 +



A 0i ui ,

(1)

i

ui = u0i +



Ai j u j .

(2)

j=i

In Eqs. (1) and (2), u is the total field, u0 is the incident field, ui is the field exciting particle i centered at Ri , A0i = A0i (r, Ri ) is a (linear) operator describing the scattering from particle i to the observation point r, u0i is the incident field at the origin of particle i, and the operator Ai j = Ai j (Ri , R j ) describes the scattering from particle j to particle i. Note that the vectors r and Ri extend from the common origin O of the laboratory coordinate system. Taking the configuration average of Eqs. (1) and (2) under the assumptions that the positions of all the particles are equally probable within the   volume   V, and that the Lax quasi-crystalline approximation u j = u j applies, yields ij

j



 

u = u0 + n0

 

3

A 0i u i d Ri ,

(3)

i



 

ui = u0i + n0 i

 

A i j u j g( R i j ) d R j , 3

(4)

j

 

where n0 = N/V is the particle number concentration, ui is the i conditional configuration average of ui with the position of par  ticle i held fixed, ui is the conditional configuration average of ij

ui with the positions of particles i and j held fixed, and g(Ri j ) = g(Ri j ) = g(Ri , R j ) with Ri j = Ri − R j is the pair correlation function. solve the integral equation (4) for   For dense media, we first  ui , and then compute u from the integral representation (3). i For sparse media, we consider the iterated solution of Eq. (4) with g( R i j ) = 1 ,



 

ui = u0i + n0 i

3



Ai j u0i d R j + n20

3

3

Ai j A jk u0k d Rk d R j + · · · , (5)

sum   up the series (5) (assuming that it converges), and compute u from Eq. (3). For sparse media, an alternative approach is the following [7– 9]: (i) consider the iterated solution of Eq. (2) under the Twersky approximation,

ui = u0i +

 j=i

Ai j u0 j +



Ai j A jk u0k + · · · ;

(6)

j=i k=i, j

(ii) insert Eq. (6) in Eq. (1) and obtain an order-of-scattering expansion for the total field

+

  i

 i

87

A 0i A i j u0 j

j=i

A0i Ai j A jk u0k i · · · ;

(7)

j=i k=i, j

(iii) take the configuration average of the expansion (7) under the assumption of independent particle positions and get

  3 3 3 u = u0 + n0 A0i u0i d Ri + n20 A0i Ai j u0 j d R j d Ri  3 3 3 + n30 A0i Ai j A jk u0k d Rk d R j d Ri + · · · ;

 

(8)

and finally, (iv) sum up the series (8). Our paper is organized as follows. In Section 2, we summarize the results of Ref. [1]. Section 3 deals with the computation of the coherent field for a dense distribution of particles. Specifically, the coherent fields reflected and transmitted by the layer, and the coherent field inside the layer will be computed. For the latter, several solution methods based on the approaches by Waterman and Truell [2], and Tsang and Kong [6], as well as on a simplified approach relying on a so called sparse-medium approximation for the integration domain, will be examined. The results are then particularized to the cases of normal incidence and a semi-infinite discrete random medium. In Section 4 we present a self-consistent derivation of the coherent field for a sparse distribution of particles, and discuss the Twersky and Foldy approximations. The paper is concluded with a short discussion. 2. Recapitulation of the results of Ref. [1]. In this section, we summarize the results established in Ref. [1]. We consider the problem of electromagnetic scattering by a discrete random medium. More specifically, we consider a group of N identical spherical particles of radius a centered at quasi-random positions R1 , R2 ,…, RN ∈ D, where the domain D is a laterally infinite plane-parallel layer with imaginary (non-scattering) boundaries z = 0 and z = H. The wavenumbers of the non-absorbing, non-magnetic background medium and the particles are k1 and k2 = mk1 , respectively, where m is the relative refractive index. We denote by f = n0V0 the particle volume concentration, where n0 = N/V is their number concentration, V is the cumulative volume occupied by the particles, and V0 = (4/3 )π a3 is the volume of each particle. The particulate medium is illuminated by a plane electromagnetic wave with the propagation direction given by the unit vector  s = s(θ0 , ϕ0 ) and the amplitude E 0 ( s ), that is,

E0 (r ) = E 0 ( s )ejk1s·r ,

(9)

( ( E 0 ( s ) = E0θ θ s ) + E0ϕ ϕ s ),

(10)

√ where j = −1; θ and ϕ hereinafter denote the corresponding spherical-coordinate angles; and E0θ and E0ϕ are the polarized components of the amplitude vector. Under the quasi-crystalline approximation, the   conditional configuration-averaged exciting field coefficients ei (Ri ) satisfy i the integral equation

  ei

i

(Ri ) = ejk1s ·Ri e0 

+ n0

D−D2a (Ri )

 

Q ( k1 Ri j ) e j

j

(R j )g(Ri j ) d3 R j ,

(11)

where D2a (Ri ) is a complete or a truncated sphere of radius 2a centered at the origin of particle i; k1s = k1 s = k1s⊥ + k z is the wave vector of the incident field; k1z (k1s⊥ ) = 1z (k1s⊥ ) k21 − k21s⊥ ; e0 = 4π E 0 ( s ) · x ( s ) is the vector of the incident field coefficients; x( r ) = [(−j )n mmn ( r ), j(−j )n nmn ( r )]T is a vector concatenating the vector spherical harmonics mmn ( r ) and

88

A. Doicu and M.I. Mishchenko / Journal of Quantitative Spectroscopy & Radiative Transfer 230 (2019) 86–105

T (k R )T; T (k R ) is the translation manmn ( r ), Q(k1 Ri j ) = T31 1 ij 31 1 i j trix relating the vectors of radiating and regular vector spherical wave functions X3 (k1 rj ) and X1 (k1 ri ), respectively, that is, X3 (k1 r j ) = T31 (k1 Ri j )X1 (k1 ri ) for r j = ri + Ri j and ri < Rij ; and T is the particle-centered transition matrix of a spherical particle. Also, the asterisk denotes a complex-conjugate value, while T stands for “transposed”. For the definitions of the vector spherical harmonics and the vector spherical wave functions we refer to Appendix A of Ref. [1]. To solve the integral equation (11) we look for a solution in the form

  ei

( Ri ) = i



ejKb ·Ri eb ,

(12)

e iη = i

eb = n0 (Jb1a + Jb2a )Teb , b = ±;

(13)

2. two inhomogeneous systems of equations corresponding to the generalized Ewald–Oseen extinction theorem:



n0 Jb2z0 Teb = 0,

(14)

b=±

ejbKz (K⊥ )H Jb2zH Teb = 0;

(15)

b=±

Jb1a ,

The expressions for the matrices in Ref. [1]. For the η-polarized incident field

Jb1a ,

Jb2z0 ,

and

are given

in which case, the equations of the generalized Lorenz–Lorentz law and the generalized Ewald–Oseen extinction theorem read as ebη = n0 (Jb1a + Jb2a )Tebη , b = ± and e0η +

eiη = ejk1s ·Ri e0η + n0

n0 Jb2z0 Tebη = 0,

(21)

(22)



ejbKz H Jb2zH Tebη = 0,

(23)

b=±

respectively. The subsequent procedure is to express the components e1bηmn and e2bηmn of the vector ebη = [e1bηmn , e2bηmn ]T as n ( e1bηmn = 4π j η sb ) · m−mn ( sb )x1bηn ,

n+1

e2bηmn = −4π j

( η sb ) · n−mn ( sb )x2bηn ,

(24)

(25)

and to solve Eqs. (21)–(23) for the azimuth-independent coefficients x1bηn and x2bηn .

The homogeneous systems of equations of the generalized Lorenz–Lorentz law are used to derive a dispersion equation for the effective wavenumber K. In Ref. [1] it was shown that the two homogeneous systems of equations (21) are identical to the homogeneous system of equations for a semi-infinite discrete random medium at normal incidence, and that the dispersion equation is direction and polarization independent. After computing the effective wavenumber K, the components of the vectors ebη can be expressed in terms of two arbitrary constants, one for e+η and the other one for e−η . These two constants are determined from the two inhomogeneous systems of equations of the generalized Ewald–Oseen extinction theorem (22) and (23) provided that each system of equations reduces to a single scalar equation. To derive these scalar equations we used the addition theorem for vector spherical harmonics, which yields the relations (cf. Appendix B of Ref. [1])

(26)

  ( ( [nmn ( k )  n−mn ( k  )] · θ k ) = −χn n(n + 1 )Nn ( k · k )θ k) m





D−D2a (Ri )



Q ( k 1 R i j ) e j η g( R i j ) d R j . 3

j

In Eq. (17), we have e0η = 4π xη ( s ) and x( r) = meaning that the components

e10ηmn

[e10ηmn , e20ηmn ]T are computed as n ( e10ηmn = 4π j η s ) · m−mn ( s ),

n+1



(16)

the polarized components eiη of ei , defined through the relation  ei = η=θ ,ϕ E0η eiη , satisfy the integral equation

e20ηmn = −4π j

(20)

b=±

m

Jb2zH

( E0η ( r ) = η s )ejk1s·r ,

i

ejKb ·Ri ebη ,

  ( ( [mmn ( k )  m−mn ( k  )] · θ k ) = χn n(n + 1 )Mn ( k · k )θ k ),

and 3. the Snell law: K sin θT = k1 sin θ0 and ϕT = ϕ0 .

 



b=±

1. two homogeneous systems of equations corresponding to the generalized Lorenz–Lorentz law:



 

b=±

where b stands for the signs + and − , and the vectors eb are unknown and have to be determined. In terms of the effective wavenumber in the particulate medium K, the wave vectors Kb in  Eq.  (12) are defined by Kb = K sb = K⊥ + bKz (K⊥ )zˆ , where Kz (K⊥ ) = K 2 − K⊥2 ,  s+ =  sT =  sT (θT , ϕT ) is the direction of the transmitted wave,  s− =  sTR =  sTR (θTR , ϕT ) is the direction of the transmitted wave reflected by the upper boundary, and θTR = π − θT . The analysis is performed under the assumptions that the representation (12) is valid for all zi in D even in the critical domains 0 ≤ zi < 2a and H − 2a < zi ≤ H. The solution obtained under these assumptions is referred to as the zeroth-order solution. Substituting Eq. (12) in Eq. (11), integrating over the positions of the particles, separating the upward and downward propagating waves, and then balancing the waves with the wavenumbers k1 and K, we are led to

e0 +

As in Eq. (12), the solution of the above integral equation is sought in the form

( η s ) · n−mn ( s ).

and

e20ηmn



(17)

( r )η r ), η=θ ,ϕ xη (

of the vector e0η =

(27) for a θ -polarized incidence, and

  ( ( [mmn ( k )  m−mn ( k  )] · ϕ k ) = −χn n(n + 1 )Nn ( k · k )ϕ k ), m

(28)

  ( ( [nmn ( k )  n−mn ( k  )] · ϕ k ) = χn n(n + 1 )Mn ( k · k )ϕ k) m

(29) (18)

(19)

for a ϕ -polarized incidence. In Eqs. (26)–(29), the coefficient χ n is given by

χn

1 = 2π n ( n + 1 )



2n + 1 , 2

(30)

A. Doicu and M.I. Mishchenko / Journal of Quantitative Spectroscopy & Radiative Transfer 230 (2019) 86–105

89

and for the directions  k= k(θ , ϕ ) and  k =  k (θ  , ϕ ), the functions Mn ( k · k ) and Nn ( k · k ) are given, respectively, by

Mn ( k · k ) = πn1 (x ),

(31)

Nn ( k · k ) = xπn1 (x ) −



n(n + 1 )Pn (x ),

(32)



|m|

where x =  k · k , πn1 (x ) = Pn1 (x )/ 1 − x2 , and Pn are the normalized associated Legendre functions of degree n and order m. The scalar equations are found to be

1 = −j

  8π 2 n0  b χn n(n + 1 ) k1 k1z Kz − bk1z n b=±





× Mn ( s · sb )Tn1 x1bθ n − Nn ( s · sb )Tn2 x2bθ n ,

0=

 b=±



(33)

  b ejbKz H χn n(n + 1 ) Kz + bk1z n

α ( r − Ri ) =

× Mn ( sR ·  sb )Tn1 x1bθ n − Nn ( sR ·  sb )Tn2 x2bθ n



(34)

for a θ -polarized incidence, and

1=j

n0 

8π k1 k1z 2



b=±

b Kz − bk1z



Fig. 1. Scattering by a layer of spherical particles.

0 , r ∈ Da ( Ri ), 1, r ∈ / Da ( Ri ),

where Da (Ri ) is a sphere of radius a centered at Ri , we express this statement mathematically as





E (r ) =

χn n(n + 1 )

n

(39)

N 



α ( r − Ri )

E0 ( r ) +

i=1



× Nn ( s · sb )Tn1 x1bϕ n − Mn ( s · sb )Tn2 x2bϕ n ,

(35)

+

N 



Escti (r )

i=1

N 

[1 − α (r − Ri )]Einti (r ).

(40)

i=1

0=

 b=±



  b ejbKz H χn n(n + 1 ) Kz + bk1z  n

× Nn ( sR ·  sb )Tn1 x1bϕ n − Mn ( sR ·  sb )Tn2 x2bϕ n

Using the identity N 



(36)

for a ϕ -polarized incidence. In Eqs. (33)–(36), Kz and k1z are given by

Kz = Kz (k1s⊥ ) =



K 2 − k21s⊥ =

2

K 2 − k21 sin

θ0 ,

(37)

α ( r − Ri ) = 1 −

i=1

N  [1 − α (r − Ri )],

(41)

i=1

and the relation giving the field exciting particle i, Eexci (r ) =  E0 (r ) + Nj=i Esct j (r ) for r ∈ Da (Ri ), we obtain

E ( r ) = E0 ( r ) +

N 

Escti (r )

i=1

and

k1z = k1z ( k 1s ⊥ ) =



k21 − k21s⊥ =

2

k21 − k21 sin

θ0 ,

(38)

respectively. Regarding the behavior of the coefficients x1bη,2n , the following two results have been established through a representative numerical analysis: 1. |x1−,η2n |  |x1+,η2n |, and

2. x1+,n2 ≈ 1 for small values of the volume concentration f, e.g., f = 0.01. 3. Coherent field in dense media A rigorous method for computing the coherent field in a discrete random medium has been given by Waterman and Truell [2]. This method which takes into account whether the observation point is inside or outside the discrete random layer is summarized below. If the observation point r is outside of any particle, the total field is the sum of the incident and all scattered fields, i.e.,  E0 ( r ) + N i=1 Escti (r ); while inside particle i, the total field is the internal field Einti (r ) when excited by Eexci (r ). Defining the indicator function

+

N 

[1 − α (r − Ri )][Einti (r ) − Eexci (r )].

(42)

i=1

Taking the configuration average of Eq. (42) over the particle positions (for the equations governing the configuration averaging process, we refer to Section 4 of Ref. [1]) we find that for an external observation point r located in   the domains z ≤ −a or z ≥ H + a, the coherent field Ec (r ) = E(r ) is

Ec ( r ) = E0 ( r ) + n0





D



Escti (r ) d Ri , 3

(43)

i

while for an external observation point residing in the critical domains −a < z ≤ 0 or H ≤ z < H + a, a truncated sphere of radius a should be excluded from the integration domain D (Fig. 2). For an internal observation point r situated in the domain a ≤ z ≤ H − a, the coherent field is

Ec ( r ) = E0 ( r ) + n0  + n0

Da ( r )





 D−Da (r )



Escti (r ) d Ri 3

i







[ Einti (r ) − Eexci (r ) ]d Ri , i

i

3

(44)

where Da (r) is a complete sphere of radius a centered at r; for an internal observation point residing in the critical domains 0 ≤ z < a or H − a < z ≤ H, Da (r) is a truncated sphere of radius a (Fig. 3).

90

A. Doicu and M.I. Mishchenko / Journal of Quantitative Spectroscopy & Radiative Transfer 230 (2019) 86–105

problem slightly in the sense that the single particle probability density function p(Ri ) = 1/V is replaced by

1 a ( R i ; r ) , V − V0

p( R i ; r ) =

1 , Ri ∈ D − Da ( r ) , 0, otherwise.

a ( R i ; r ) =

Fig. 2. For an external observation point P and for −a < z ≤ 0 or H ≤ z < H + a, the domain of integration is D with a truncated sphere of radius a excluded. O is the origin of the laboratory coordinate system.

The new probability density function (50), which is similar to the hole-correction approximation, can be regarded as a conditional probability given the position of the observation point r; it is also normalized, and we have

 i



D









 

i

i

(i ) Escti (r ) = Esct (r ) = XT3 (k1 ri )T ei , i i i









 

i

i

(i ) Eexci (r ) = Eexc (r ) = XT1 (k1 ri ) ei , i i i









(i )

Einti (r ) = Einti (ri ) = i

 

i

XT1

 

(k2 ri )Tint ei i ,

(45)

i

where ri = r − Ri , ei is the solution of the integral equation (11), i and, for a spherical particle, Tint is the particle-centered “transition matrix” relating the expansion coefficients of the internal field to those of the exciting field. The first integral in Eq. (44) corresponds to the configurations in which the observation point r is external to all particles. The first term in the second integral of Eq. (44) gives the internal contribution (the point r is inside the particle), while the second term cancels out the external point contribution since the probability that r is an outside point diminishes [2]. As for the exciting field coefficients, we will consider only the zeroth-order solution for the coherent field, in which case the exclusion domain Da (r) is a complete sphere of radius a. In particular, we will use Eq. (43) to compute the coherent fields reflected and transmitted by the layer, and Eq. (44) to compute the coherent field inside the layer. To integrate over all positions of particle i we use a local coordinate system centered at the observation point. In this regard, we make the change of variable p = Ri − r = −ri , implying p = p⊥ + pz z with p⊥ = Ri⊥ − r⊥ and pz = zi − z, and use the symmetry relations (cf. Eq. (256) of Ref. [1])

M3mn (−k1 p ) = (−1 )n M3mn (k1 p ),

(−k1 p ) = (−1 )

n+1

N3mn

( k1 p ).

It should be pointed out that the coherent field component given by the first integral in Eq. (44) alters the statistics of the

 D−Da (r )



Escti (r ) d Ri



3

i



Escti (r ) d Ri . 3

i

The coherent field reflected by the layer is

ERη ( r ) = n0





D

b=± mn

Tn1 e1bηmn M3mn (k1 ri )



+ Tn2 e2bηmn N3mn (k1 ri ) ejKb ·Ri d Ri , z ≤ −a. 3

(51)

For z ≤ −a < 0 ≤ zi , we use the integral representations of the radiating vector spherical wave functions in terms of plane electromagnetic waves, that is (p = Ri − r = −ri ) [10]

M3mn (k1 p ) =

N3mn (k1 p ) =



1

mmn ( k+ )ej[k⊥ ·p⊥ +k1z (k⊥ )(zi −z )]

2π j

n



1

2

d k⊥ , k1 k1z ( k ⊥ ) (52)

jnmn ( k+ )ej[k⊥ ·p⊥ +k1z (k⊥ )(zi −z )]

2π j

n

k+ = k⊥ + k1z (k⊥ ) z,

2

d k⊥ , (53) k1 k1z ( k ⊥ ) (54)

to obtain



M3mn (k1 ri )ejKb ·Ri d Ri 3

D

=−

(48) (49)



3.1. The coherent field reflected by the layer

 N3mn

3

D−Da (r )

(46) (47)

N V − V0  ≈ n0

p(Ri ; r ) Escti (r ) d Ri =

In Eqs. (43) and (44), the conditional configuration averages of the fields are given by



(50)

2π 1 − ej(bKz +k1z )H jk1R ·r e mmn ( sR ), n+1 k k (bK + k ) z 1 1z 1z j

(55)

N3mn (k1 ri )ejKb ·Ri d Ri 3

D

=−

2π 1 − ej(bKz +k1z )H jk1R ·r e jnmn ( sR ), n+1 k k (bK + k ) z 1 1z 1z j

(56)

Fig. 3. For an internal observation point P, the exclusion domain Da (r) is a complete sphere of radius a if a ≤ z ≤ H − a (left), and a truncated sphere of radius a if 0 ≤ z < a (middle) or H − a < z ≤ H (right).

A. Doicu and M.I. Mishchenko / Journal of Quantitative Spectroscopy & Radiative Transfer 230 (2019) 86–105

with k1R = k1 sR = k1s⊥ − k1z z. Inserting Eqs. (24), (25), (55) and (56) in Eq. (51), we find that the coherent reflected field is computed according to

8π 2 n0 jk1R ·r  b ERη ( r ) = j e [1 − ej(bKz +k1z )H ] k1 k1z Kz + bk1z b=±  ( × [mmn ( sR )  m−mn ( s b )] · η sb )Tn1 x1bηn mn



(57)

M3mn (k1 p ) =



1 2π j

n

   8π 2 n0  b 1 − ej(bKz +k1z )H χn n(n + 1 ) k1 k1z Kz + bk1z n b=±

× [Mn ( sR ·  sb )Tn1 x1bθ n − Nn ( sR ·  sb )Tn2 x2bθ n ]

R0ϕ = − j



8π n0  b 1 − ej(bKz +k1z )H k1 k1z Kz + bk1z 2

(59)



b=±

8π k1 k1z

×

(

b=±

b Kz + bk1z

)

[Mn  sR ·  sb Tn1 x1bθ n

R0ϕ = − j



−

N3mn (k1 p ) =

(67)

2

d k⊥ , k1 k1z ( k ⊥ ) (68)



1

mmn ( k− )ej[k⊥ ·p⊥ −k1z (k⊥ )(zi −z )]

2π j

n

jnmn ( k− )ej[k⊥ ·p⊥ −k1z (k⊥ )(zi −z )]

2

d k⊥ , (69) k1 k1z ( k ⊥ )

)

Nn  sR ·  sb Tn2 x2bθ n ],

M3mn (k1 ri )ejKb ·Ri d Ri 3

2π 1 − ej(bKz −k1z )H jk1s ·r e mmn ( s ), n+1 k k (bK − k ) z 1 1z 1z j

(71)

(61)  D

b=±

(62)

and in a second step,

   8π 2 n0 1 − e2jKz H χn n(n + 1 ) k1 k1z (Kz + k1z ) n

× [Mn ( sR ·  sT )Tn1 x1+θ n − Nn ( sR ·  sT )Tn2 x2+θ n ],



= −

  8π 2 n0  b χn n(n + 1 ) k1 k1z Kz + bk1z n

× [Nn ( sR ·  sb )Tn1 x1bϕ n − Mn ( sR ·  sb )Tn2 x2bϕ n ],

(70)

yielding

D

χn n(n + 1 )

= −

N3mn (k1 ri )ejKb ·Ri d Ri 3

2π 1 − ej(bKz −k1z )H jk1s ·r e jnmn ( s ), n+1 k k (bK − k ) z 1 1z 1z j

(72)

along with Eqs. (24) and (25) to obtain

(63)

   8π 2 n0 R0ϕ = − j 1 − e2jKz H χn n(n + 1 ) k1 k1z (Kz + k1z ) n

( ETη ( r ) = η s )ejk1s ·r + j ×

 b=±

(64)

For a θ -polarized incidence, we plot in Fig. 4 the reflection coefficient |R0θ | together with the reflection coefficient |R0θ EFA | for a homogeneous layer with thickness H and wavenumber K placed in a background medium with wavenumber k1 . The latter corresponds to the effective field approximation, and is computed as

(1 − e2jKz H )rη R0ηEFA = , 1 − (rη )2 e2jKz H

3

(60)

n

× [Nn ( sR ·  sT )Tn1 x1+ϕ n − Mn ( sR ·  sT )Tn2 x2+ϕ n ].



k− = k⊥ − k1z (k⊥ ) z,



(

Tn1 e1bηmn M3mn (k1 ri )

χn n(n + 1 )

n

for a ϕ -polarized incidence. Furthermore, from Eqs. (34)–(36) of the generalized Ewald–Oseen extinction theorem, we find in a first step

n0 

D



× [Nn ( sR ·  sb )Tn1 x1bϕ n − Mn ( sR ·  sb )Tn2 x2bϕ n ]

2



For z ≥ H + a > H ≥ zi , we use the integral representations (p = Ri − r = −ri )

(58)

and using the addition theorem for vector spherical harmonics (cf. Eqs. (26) and (27) and Eqs. (28) and (29)), we obtain

for a θ -polarized incidence, and



+ Tn2 e2bηmn N3mn (k1 ri ) ejKb ·Ri d Ri , z ≥ H + a.

( ( ERη ( r ) = η sR )R0η ejk1R ·r = η sR )R0η ejk1s⊥ ·r e−jk1z z ,

R0θ = j

The coherent field transmitted by the layer is

b=± mn

Defining the reflection coefficient R0η through the relation

R0θ = j

3.2. The coherent field transmitted by the layer

ETη ( r ) = E0η ( r ) + n0

( + [nmn ( sR )  n−mn ( s b )] · η sb )Tn2 x2bηn .

R0θ = j

91

×

 mn

8π 2 n0 jk1s ·r e kk1z



b 1 − ej(bKz −k1z )H Kz − bk1z



( {[mmn ( s )  m−mn ( s b )] · η sb )Tn1 x1bηn

( + [nmn ( s )  n−mn ( s b )] · η sb )Tn2 x2bηn }.

(73)

Defining the transmission coefficient Tη2 by

( ( ETη ( r ) = η s )Tη2 ejk1s ·r = η s )Tη2 ejk1s⊥ ·r ejk1z z ,

(74)

(65)

where

and employing the addition theorem for vector spherical harmonics (cf. Eqs. (26) and (27) and Eqs. (28) and (29)), we get

rη =

Tθ2 = 1 + j

⎧ K cos θ0 − k1 cos θ ⎪ , η = θ, ⎪ ⎨ K cos θ0 + k1 cos θ ⎪ ⎪ ⎩ k1 cos θ0 − K cos θ , η = ϕ . k1 cos θ0 + K cos θ

(66)

The results show that (i) in general and in agreement with the findings of Ref. [11], the reflection coefficients are small; (ii) they increase with the volume concentration; and (iii) in the low frequency limit, the effective field approximation is valid.





8π 2 n0  b 1 − ej(bKz −k1z )H k1 k1z Kz − bk1z b=±  

 × χn n(n + 1 ) Mn ( s · sb )Tn1 x1bθ n − Nn ( s · sb )Tn2 x2bθ n n

(75) for a θ -polarized incidence, and

92

A. Doicu and M.I. Mishchenko / Journal of Quantitative Spectroscopy & Radiative Transfer 230 (2019) 86–105

Fig. 4. Reflection coefficients |R0θ | (solid curve) and |R0θ EFA | (solid curve with circles) for a θ -polarized incidence as functions of the size parameter k1 a and for two values of the volume concentration f and the layer thickness H. The incidence angle is θ0 = 30◦ , the wavenumber of the background medium is k1 = 10 μm−1 , the relative refractive index of the particles is m = 1.33, and the maximum expansion order in the expansion (61) is Nrank = 15.





8π 2 n0  b 1 − ej(bKz −k1z )H k1 k1z Kz − bk1z b=±  

 × χn n(n + 1 ) Nn ( s · sb )Tn1 x1bϕ n − Mn ( s · sb )Tn2 x2bϕ n

Tϕ2 = 1 − j

×

 n



 χn n(n + 1 ) Nn ( s · sb )Tn1 x1bϕ n − Mn ( s · sb )Tn2 x2bϕ n , (78)

and then

n

(76) for a ϕ -polarized incidence. Furthermore, invoking Eqs. (33)–(35) of the generalized Ewald–Oseen extinction theorem, we find

×

 n

8π 2 n0  b Tθ = − j ej(bKz −k1z )H k1 k1z Kz − bk1z b=±  

 × χn n(n + 1 ) Mn ( s · sb )Tn1 x1bθ n − Nn ( s · sb )Tn2 x2bθ n ,



Tθ2 = e−j(Kz +k1z )H 1 + j



8π 2 n0 1 − e2jKz H k1 k1z (Kz − k1z )





 χn n(n + 1 ) Mn ( s · sT )Tn1 x1+θ n − Nn ( s · sT )Tn2 x2+θ n , (79)

2



Tϕ2 = e−j(Kz +k1z )H 1 + j

n

(77)

×





8π 2 n0 1 − e2jKz H k1 k1z (Kz − k1z )





 χn n(n + 1 ) Nn ( s · sT )Tn1 x1+ϕ n − Mn ( s · sT )Tn2 x2+ϕ n .

n

(80) Tϕ2 = j

n0 

8π k1 k1z 2

b=±

b ej(bKz −k1z )H Kz − bk1z

Note that although the representations (77) and (78) are mathematically equivalent with the representations (79) and (80),

A. Doicu and M.I. Mishchenko / Journal of Quantitative Spectroscopy & Radiative Transfer 230 (2019) 86–105

93

Fig. 5. Transmission coefficients |Tθ2 | (solid curve) and |Tθ2EFA | (solid curve with circles) for a θ -polarized incidence as functions of the size parameter k1 a and for two values of the volume concentration f and the layer thickness H. The incidence angle is θ0 = 30◦ , the wavenumber of the background medium is k1 = 10 μm−1 , the relative refractive index of the particles is m = 1.33, and the maximum expansion order in the expansion (77) is Nrank = 15.

respectively, the exponential term exp(−jKz H ) shows that the last two are numerically unstable for large H. Therefore, in numerical implementations, the representations (77) and (78) should be used. The plots in Fig. 5 show that the transmission coefficient |Tθ2 | agrees with the transmission coefficient |Tθ2EFA | corresponding to a homogeneous layer, where

Tη2EFA =

[1 − (rη )2 ]ej(Kz −k1z )H . 1 − (rη )2 e2jKz H

(81)

Moreover, for higher frequencies, the transmission coefficients decrease with increasing (i) the volume concentration, and (ii) the geometrical thickness of the layer. To justify the result Tη2 ≈ Tη2EFA , we consider Eq. (33) of the generalized Ewald–Oseen extinction theorem, in which we employ the assumption |x1−,η2n |  |x1+,η2n |; we then get

1= −j



  8π 2 n0 1 χn n(n + 1 ) k1 k1z Kz − k1z n



× Mn ( s · sT )Tn1 x1+θ n − Nn ( s · sT )Tn2 x2+θ n ,

(82)

so that from Eqs. (79) and (80), we find

Tη2 ≈ ej(Kz −k1z )H .

(83)

On the other hand, it is obvious from Eq. (81) that the same approximation is valid for Tη2EFA under the assumption |rη2 |  1. It should be pointed out that this result was exploited by Gustavsson et al. [11] to compute the effective wavenumber in the framework of the integral equation method developed by Kristensson [12] to analyze the coherent scattering by a discrete random layer. From the results of Figs. 4 and 5, we deduce that only the reflection data contain information of the particulate structure of the layer (except for low frequencies, when the effective field approximation applies) [11].

94

A. Doicu and M.I. Mishchenko / Journal of Quantitative Spectroscopy & Radiative Transfer 230 (2019) 86–105

3.3. The coherent field inside the layer In this section we present several methods for computing the coherent field inside the layer. 3.3.1. The method of Fikioris and Waterman We consider the coherent field Ec (r ) given by Eq. (44). From Eqs. (45)–(47) it is apparent that the computation of Ec (r ) requires the calculation of the integrals





M3mn (k1 ri )ejKb ·Ri d Ri , 3

D−Da (r )

and



M1mn (k1,2 ri )ejKb ·Ri d Ri , 3

Da ( r )

N3mn (k1 ri )ejKb ·Ri d Ri 3

D−Da (r )



N1mn (k1,2 ri )ejKb ·Ri d Ri . 3

Da ( r )

(84)

(85)

These volume integrals are transformed into surface integrals by means of the Green’s theorem of the second kind. For the scalar ,3 function f = exp(jKb · Ri ) and the vector function G = M1mn (kri ) or ,3 N1mn (kri ), solving the Helmholtz equation in the domain D0 with the wavenumbers K and k, respectively, the Green’s theorem yields

 f G dV =

D0

K2

1 − k2

 ∂G ∂f f −G dS, ∂ n ∂ n S0





D−Da (r )

dV →

Sz0 ∪SzH

 n+ z dS −



Sa ( r )

 n+ a dS,

(87)

while the volume integral over Da (r) will be transformed into one surface integral, i.e.,





Da ( r )

dV →

Sa ( r )

 n+ a dS.

(88)

 n+ a

Here, is the outward-pointing unit vector normal to the spheri+   cal surface Sa (r) of radius a centered at r, and  n+ z = −z and nz = z are the outward-pointing unit vectors normal to the plane surfaces Sz0 and SzH , respectively. Let us denote by 1. IbzMmn (r ) and IbzNmn (r ) the integrals over the union of the planes Sz0 ∪ SzH involving M3mn (k1 ri ) and N3mn (k1 ri ), respectively; and by b αb 2. Iα aMmn (r ) and IaNmn (r ) the integrals over the spherical surface α Sa (r) involving Mα mn (kα ri ) and Nmn (kα ri ), respectively. In the latter case, the vector spherical wave functions Mα mn (kα ri ) and Nα mn (kα ri ) are defined as follows: 1 α 1 1. for α = 1, we have kα = k1 , Mα mn = Mmn , and Nmn = Nmn ; 1 , and Nα = N1 ; 2. for α = 2, we have kα = k2 , Mα = M mn mn mn mn and 3 α 3 3. for α = 3, we have kα = k1 , Mα mn = Mmn , and Nmn = Nmn .

With this notation, we express the coherent field as

Ecη ( r ) = Ecη 1 ( r ) + Ecη 2 ( r ),

(89)

2π j

n

α =1



1

N3mn (k1 p ) =

2π j

n

d k⊥ , k1 k1z ( k ⊥ )

jnmn ( k− )ej[k⊥ ·p⊥ +k1z (k⊥ )z]

d k⊥ , k1 k1z ( k ⊥ )



1

M3mn (k1 p ) =

2π j

n



1

N3mn (k1 p ) =

2π j

n

2

d k⊥ , k1 k1z ( k ⊥ )

jnmn ( k+ )ej[k⊥ ·p⊥ +k1z (k⊥ )(H−z )]

d k⊥ , k1 k1z ( k ⊥ )

where kb = k⊥ + bk1z (k⊥ ) z for b = ± and k1z (k⊥ ) = obtain

2



k21 − p2⊥ . We

2π ejk1s ·r mmn ( s) k1 k1z (bKz − k1z ) 2π + n+1 ejk1s⊥ ·r⊥ ejk1z (H−z ) ejbKz H mmn ( sR ) j k1 k1z (bKz + k1z ) (93)

IbzMmn (r ) = −

j

n+1

and

2π j ejk1s ·r nmn ( s) k1 k1z (bKz − k1z ) 2π j + n+1 ejk1s⊥ ·r⊥ ejk1z (H−z ) ejbKz H nmn ( sR ), j k1 k1z (bKz + k1z ) (94)

IbzNmn (r ) = −

j

n+1

so that the field Ecη1 (r ) computes as

( Ecη 1 ( r ) = η s )ejk1s ·r − ×

(90) ×

(91)

2

mmn ( k+ )ej[k⊥ ·p⊥ +k1z (k⊥ )(H−z )]

 1

2π n0 jk1s ·r  1 e k1 k1z bKz − k1z b=±



Tn1 e1bηmn mmn n+1

j

( s ) + jTn2 e2bηmn nmn ( s)



2π n0 jk1s⊥ ·r⊥ jk1z (H−z )  ejbKz H e e k1 k1z bKz + k1z  1 n+1

j



b=±



Tn1 e1bηmn mmn ( sR ) + jTn2 e2bηmn nmn ( sR ) .

Next we proceed as follows.

(−1 )α Eαcη2 (r ),

2

mmn ( k− )ej[k⊥ ·p⊥ +k1z (k⊥ )z]

while for z ≤ H − a and in the plane pz = H − z (zi = H), we use

mn 3 



1

M3mn (k1 p ) =

+

and

(92)

3.3.1.1. Computation of Ecη1 . To compute IbzMmn (r ) and IbzNmn (r ), we employ the integral representations of the radiating vector spherical wave functions in terms of plane electromagnetic waves; in particular, for z ≥ a and in the plane pz = −z (zi = 0), we use (p = Ri − r = −ri )

mn

b=± mn

Ecη 2 ( r ) =

b=± mn



Tα1n e1bηmn IαaMbmn (r ) + Tα2n e2bηmn IαaNbmn (r ) ,

respectively. In Eq. (92), we used the convention that (i) Tα1n,2 = 1 1,2 for α = 1 (exciting field), (ii) Tα1n,2 = Tint n for α = 2 (internal field), 1,2 1,2 and (iii) Tα n = Tn for α = 3 (scattered field). In the following, we will calculate the components Ecη1 (r ) and Ecη2 (r ) of the coherent field separately. In fact, we will show that Ecη1 (r ) = 0, so that Ecη (r ) = Ecη2 (r ). In view of Eq. (90), the result Ecη1 (r ) = 0 means that the waves produced by the particles situated at the lower and upper boundaries of the medium extinguishes the incident wave.

where Ecη1 (r ) and Ecη2 (r ) are given by

Ecη 1 ( r ) = E0η ( r )  

 + n0 Tn1 e1bηmn IbzMmn (r ) + Tn2 e2bηmn IbzNmn (r ) ,

 

(86)

where S0 is the boundary surface of D0 , and  n is the outwardpointing normal unit vector to S0 . Actually, the volume integral over D − Da (r ) will be transformed into three surface integrals (we neglect the integrals over the cylindrical surfaces at infinity), i.e.,



Eαcη2 (r ) = n0

(95)

A. Doicu and M.I. Mishchenko / Journal of Quantitative Spectroscopy & Radiative Transfer 230 (2019) 86–105

1. Let E1cη1 (r ) be the cumulative contribution of the first two terms on the right-hand side of Eq. (95). By means of Eqs. (24) and (25), this contribution can be written as

( E1cη1 (r ) = η s )ejk1s ·r + j ×

8π 2 n0 jk1s ·r  1 e k1 k1z bKz − k1z

mn

( {[mmn ( s )  m−mn ( s b )] · η sb )Tn1 x1bηn

( + [nmn ( s )  n−mn ( s b )] · η sb )Tn2 x2bηn }.

(96)

Using Eqs. (26)–(29) along with Eqs. (33)–(35) of the generalized Ewald–Oseen extinction theorem, we find that E1cη1 (r ) = 0.

2. Let E2cη1 (r ) be the contribution of the third term on the right-hand side of Eq. (95). Again, by means of Eqs. (24) and (25), we have

E2cη1 (r ) = − j

8π 2 n0 jk1s⊥ ·r⊥ jk1z (H−z )  ejbKz H e e k1 k1z bKz + k1z b=±

 ( × sR )  m−mn ( sb )] · η sb )Tn1 x1bηn [mmn ( mn

(97)

Thus, the field Ecη1 (r ) = have Ecη (r ) = Ecη2 (r ).

+ E2cη1 (r )



2 β vαmn ( p ) · vm n ( p) d  p = δ pq δm ,−m δn n ,

vanishes, and we

 p × lmn ( p ) = 0,  p × nmn ( p ) = −mmn ( p ),  p × mmn ( p ) = nmn ( p ). (103) Note that the relations (99)–(101) have been derived by using the representations of the vector spherical wave functions in terms of vector spherical harmonics (cf. Appendix A of Ref. [1]) and the differential equation satisfied by the spherical Bessel and Hankel functions zn (x), i.e.,

IαaMbmn (r ) = 4π j InMα (a )mmn ( sb ),

(105)

IαaNbmn (r ) = − 4π j

(106)

n+1

where

InNα (a ) =





(Kp ) ·

∂

Mα

(k p ) ∂ (kα p) mn α  ∂ 2 − KMαmn (kα p ) · M1m n (Kp ) d  p ∂ (K p )   = kα  p · [M1m n (Kp ) × Nαmn (kα p )]  2 + K p · [N1m n (Kp ) × Mαmn (kα p )] d  p, kα M1m n

∂ Nα ( k p ) ∂ (kα p) mn α  ∂ 2 − KNαmn (kα p ) · N1m n (Kp ) d  p ∂ (K p )   = kα  p · [N1m n (Kp ) × Mαmn (kα p )]  2 + K p · [M1m n (Kp ) × Nαmn (kα p )] d  p, 1 K 2 − k2α





kα L1m n (Kp ) ·

∂ Nα ( k p ) ∂ (kα p) mn α



K

a − k2α

kα

(107)

jn (Ka )[kα aznα (kα a )]



kα α z (kα a )[Ka jn (Ka )] , K n

(108)

a  n(n + 1 )znα (kα a ) jn (Ka ), Kkα



M1mn (kα p )ejKb ·p d p and 3

(109)

N1mn (kα p )ejKb ·p d p, 3

Da

(110)

where α = 1, 2 and Da is the domain occupied by a spherical particle of radius a centered at the origin of the coordinate system, reduces to the calculation of the following integrals:



 Da

 Da

N1m n



Da

L1m n (Kp ) ·



(111)



M1mn (kα p ) 3 d p, N1mn (kα p )

(Kp ) ·

and



M1mn (kα p ) 3 d p, N1mn (kα p )

M1m n (Kp ) ·



(100)

jn (Ka )[kα aznα (kα a )]

zn1 (x ) = zn2 (x ) = jn (x ), zn3 (x ) = hn (x ), and jn (x) and hn (x) are the spherical Bessel and Hankel functions, respectively. Thus, the fields Eα (r ) result from Eq. (92) along with Eqs. (105)–(109). cη 2 Before continuing our analysis, we make a short comment. In view of Eq. (98), the volume integrals

Da

kα N1m n (Kp ) ·

and

InLα (a ) =



(99)

K2 −

+ jm−m n ( sb )  M1m n (Kp ) + n−m n ( sb )  N1m n (Kp )]; (98) (ii) the dyadic identities a = I · a and (ab ) · c = a(b · c ); (iii) the relations



a K 2 − k2α

[InLα (a )lmn ( sb ) + InNα (a )nmn ( s b )] ,

− znα (kα a )[Ka jn (Ka )] ,

m n



(104)

We obtain

InMα (a ) =

IejKb ·p = IejKsb ·p  n +1 = − 4π j [l−m n ( sb )  L1m n (Kp )

(102)

and finally, (vi) the identities

b αb 3.3.1.2. Computation of Ecη2 . To compute Iα aMmn (r ) and IaNmn (r ), we make the change of variable p = Ri − r = −ri , and use (i) the expansion (cf. Appendix A)



α , β = 1, 2, 3;

n

Using Eqs. (26)–(29) along with Eqs. (34)–(36) of the generalized Ewald–Oseen extinction theorem, we find that E2cη1 (r ) = 0. E1cη1 (r )

(101)

x2 zn (x ) + 2xzn (x ) + [x2 − n(n + 1 )]zn (x ) = 0.



( + [nmn ( sR )  n−mn ( sb )] · η sb )Tn2 x2bηn .



∂

2 − KNαmn (kα p ) · L1   (Kp ) d  p ∂ (K p ) m n  1 2 = − Mαmn (kα p ) · [ p × L1m n (Kp )]d  p; kα

(iv) the representation of the vector spherical wave functions in terms of the vector spherical harmonics v1mn ( p ) = lmn ( p ), v2mn ( p) = nmn ( p ), and v3mn ( p ) = mmn ( p ); (v) the orthogonality relations (cf. Appendix A of Ref. [1])

b=±



95

(112)



M1mn (kα p ) 3 d p. N1mn (kα p )

(113)

Instead of using the Green’s theorem (86), we can compute the integrals (111)–(112) by means of the second vector Green’s theorem.

96

A. Doicu and M.I. Mishchenko / Journal of Quantitative Spectroscopy & Radiative Transfer 230 (2019) 86–105

For P(r ) = M1m n (Kr )/N1m n (Kr ), Q(r ) = M1mn (kα r )/N1mn (kα r ), and a spherical domain, we use



[ P · ( ∇ × ∇ × Q ) − Q · ( ∇ × ∇ × P )] d r  2 = a2  r · [ Q × ( ∇ × P ) − P × ( ∇ × Q )] d  r.

we find that the coherent field can be written in the form

( ( Ecη ( r ) = η sT )Tη1FW ejKsT ·r + η sTR )R1ηFW ejKsTR ·r

3

( ( =η sT )Tη1FW ejk1s⊥ ·r ejKz z + η sTR )R1ηFW ejk1s⊥ ·r e−jKz z ,

Da

(114)

On the other hand, by taking into account that ∇ × L1m n = 0, the integral (113) can be computed by means of the first vector Green’s theorem. For P(r ) = L1m n (Kr ), Q(r ) = M1mn (kα r )/N1mn (kα r ), and a spherical domain, we use



where the transmission and reflection coefficients are given by

Tη1FW

= 2π n0

 +

[ ( ∇ × P ) · ( ∇ × Q ) − P · ( ∇ × ∇ × Q )] d r  2 = a2  r · [ P × ( ∇ × Q )] d  r. 3

Iα b

(115)

ejKb ·r

b=±

  3  ( × [mmn ( sb )  m−mn ( s b )] · η sb ) (−1 )α Tα1n InMα (a ) x1bηn α =1

mn

( + [nmn ( sb )  n−mn ( s b )] · η sb )

 3 α =1

( + [lmn ( sb )  n−mn ( s b )] · η sb )

 3

α =1



m



(−1 )α Tα2n InLα (a ) x2bηn .

∂θ

(116)

(117)

where  k= k ( θ , ϕ ),  k =  k ( θ  , ϕ ), x =  k · k = cos(θ − θ  ),

∂ Pn (x ) = ∂θ 

n(n + 1 )Ln ( k,  k ),

(118)

Ln ( k,  k ) = sin(θ − θ  )πn1 (x ),

(119)

and Pn (x ) = Pn0 (x ) are the Legendre polynomials. Since for  k= k we have Ln ( k,  k ) = 0, it is readily seen that

 ( [lmn ( sb )  n−mn ( s b )] · η sb ) = 0.

(120)

m

Thus, the last term in Eq. (116) vanishes. Substituting Eqs. (26)– (29) in Eq. (116), using the special values of the functions Mn and Nn at x = 1 (cf. Eqs. (31) and (32)),

Mn (1 ) = − Nn (1 ) = ωn , with

ωn

1 = 2



n(n + 1 )(2n + 1 ) , 2

α =1

RηFW = 2π n0

 +



(−1 )α Tα1n InMα (a ) x1+ηn



(−1 )α Tα2n InNα (a ) x2+ηn

 mn

3 

α =1

 ( 2n + 1 )

3 

α =1

(124)

 (−1 )α Tα1n InMα (a ) x1−ηn



(−1 )α Tα2n InNα (a ) x2−ηn

,

(125)

respectively. From Eq. (123), we see that the coherent field Ecη (r ) is a superposition of upwelling and downwelling plane electromagnetic waves propagating in an effective medium with wavenumber K. If we extend the representations of the coherent fields reflected and transmitted by the layer in the critical domains −a < z ≤ 0 and H ≤ z < H + a, respectively, as well as the representation of the coherent field in the critical domains 0 ≤ z < a and H − a < z ≤ H, we may expect that the boundary conditions for the electric fields,

 z × Ecη ( r ) =  z × [ERη (r ) + E0η (r )] at z = 0  z × Ecη ( r ) =  z × ETη (r ) at z = H,



∂ Pn (x )    lmn ( k )  n−mn ( k ) = χn k  θ ( k ),  

3 

3 

(126)

and

(−1 )α Tα2n InNα (a ) x2bηn

Next, we use the following result. The addition theorem for the vector spherical harmonics lmn ( k ) and n−mn ( k ) in the case ϕ = ϕ  is (cf. Appendix B of Ref. [1])



( 2n + 1 )

mn

1

Iα b

3.3.1.3. Computation of Ecη . We are now in a position to obtain a complete representation for the coherent field. Inserting Eqs. (105) and (106) together with Eqs. (24) and (25) in Eq. (92), and the result in Eq. (91), we get





and

The final expressions for aMmn (r ) and aNmn (r ) are as in Eqs. (105) and (106) along with Eqs. (107)–(109); however, the proof is now more straightforward and does not make use on the identities (99)–(101).

Ecη (r ) = 16π 2 n0



α =1

Da

(123)

(121)

(122)

(127)

are satisfied. Fikioris and Waterman [3,4] proved that the zerothorder fields satisfy the boundary conditions (126) and (127) in the low frequency limit k1 a  1. Unfortunately, at higher frequencies, these boundary conditions are not satisfied. This can be inferred from Fig. 6 illustrating the relative errors in the reflection and transmission coefficients

εRθ =

||R1θ | − |R1θ FW || ||Tθ1 | − |Tθ1FW || , εTθ = , 1 |Rθ | |Tθ1 |

where |R1θ | and |Tθ1 | correspond to a model in which the boundary conditions for the electric fields are matched (cf. Eqs. (129) and (130) below). Therefore, for particles with sizes comparable to or larger than the wavelength, higher-order approximations to the fields should by computed by means of the iteration scheme described in Ref. [1] and the validity of the boundary conditions (126) and (127) for each order of approximation should be analyzed [2–4]. 3.3.2. The method of Tsang and Kong The method of Tsang and Kong [6] is a simplified approach that avoids an explicit computation of the coherent field inside the layer, and so, of higher-order approximations. In this approach, 1. the representations of the (zeroth-order) coherent fields reflected and transmitted by the layer are extended in the critical domains; 2. the coherent field inside the layer is assumed to be a superposition of plane electromagnetic waves propagating in an effective medium with wavenumber K, that is,

( ( Ecη ( r ) = η sT )Tη1 ejk1s⊥ ·r ejKz z + η sTR )R1η ejk1s⊥ ·r e−jKz z ;

(128)

A. Doicu and M.I. Mishchenko / Journal of Quantitative Spectroscopy & Radiative Transfer 230 (2019) 86–105

97

Fig. 6. Relative errors in |R1θ | and |Tθ1 | for a θ -polarized incidence as functions of the size parameter k1 a and for two values of the volume concentration f. The layer thickness is H = 10a, the incidence angle is θ0 = 30◦ , the wavenumber of the background medium is k1 = 10 μm−1 , the relative refractive index of the particles is m = 1.33, and the maximum expansion order is Nrank = 15. Note that in the case H = 100a, the relative errors have virtually the same dependency on k1 a.

3. the reflection and transmission coefficients R1η and Tη1 are computed from the boundary conditions for the electric fields (126) and (127). We obtain



 cos θ

1 0 R1θ = (1 − R0θ )e2jKz H − Tθ2 ej(Kz +k1z )H , cos θT 1 − e2jKz H Tθ1 = (1 − R0θ )

cos θ0 + R1θ cos θT

for a θ -polarized incidence, and

Tϕ1 =



(131)

(132)

for a ϕ -polarized incidence. We emphasize once again that the coefficients R1η and Tη1 are computed from the boundary conditions for the electric fields, while in the effective field approximation, the coefficients R1ηEFA and Tη1EFA , given by

R1θ EFA = −

1 + rη k1 rη (1 + rη )e2jKz H k1 , Tθ1EFA = , K 1 − (rη )2 e2jKz H K 1 − (rη )2 e2jKz H

(133)

Rϕ EFA

1 + rη rη (1 + rη )e2jKz H =− , Tϕ1EFA = , 1 − (rη )2 e2jKz H 1 − (rη )2 e2jKz H



3

D

d Ri

(135)

is typical of sparse media and is referred to as the sparse-medium approximation for the integration domain. Also note that in this case, the standard definition of the single particle probability density function p(Ri ) = 1/V is maintained. The integrals of the radiating vector spherical wave functions are computed in the sense of Cauchy’s principal value by excluding from the integration domain D a layer of thickness 2ε around the plane z, and then by letting ε → 0. We obtain



M3mn (k1 ri )ejKb ·Ri d Ri = −

(134)

are computed from the boundary conditions for both the electric and magnetic fields. For the magnetic field, the discontinuity in the coefficients R1η and Tη1 can be used to predict the effective magnetic permeability of the discrete random layer [4]. 3.3.3. Sparse-medium approximation Let us compute the coherent field Ec (r ) from Eq. (44) by neglecting the integrals over Da (r). In other words, in the course of

2π

3

D

2π

n+1

j

k1 k1z (bKz − k1z )

ejk1s ·r mmn ( s)

ejk1s⊥ ·r⊥ ejk1z (H−z ) ejbKz H mmn ( sR ) k1 k1z (bKz + k1z )  m (  mmn ( 2π sR ) mn s ) + n+1 ejk1s⊥ ·r⊥ ejbKz z − bKz − k1z bKz + k1z j k1 k1z +

n+1

j

(136)

and



and 1

3

d Ri ≈

(129)



R1ϕ = 1 + R0ϕ − Tϕ1

 D−Da (r )

(130)

1 1 + R0ϕ − Tϕ2 ej(Kz +k1z )H , 1 − e2jKz H

evaluating the integrals in Eq. (44) we treat the particles as point scatterers. Note in this respect that the approximation

D

N3mn (k1 ri )ejKb ·Ri d Ri 3

2π j ejk1s ·r nmn ( s) k1 k1z (bKz − k1z ) 2π j + n+1 ejk1s⊥ ·r⊥ ejk1z (H−z ) ejbKz H nmn ( sR ) j k1 k1z (bKz + k1z )  n (  nmn ( 2π j sR ) mn s ) + n+1 ejk1s⊥ ·r⊥ ejbKz z − . bKz − k1z bKz + k1z j k1 k1z

= −

j

n+1

(137)

Consequently, we find that the coherent field Ecη (r ) is given by

98

A. Doicu and M.I. Mishchenko / Journal of Quantitative Spectroscopy & Radiative Transfer 230 (2019) 86–105

( Ecη ( r ) = η s )ejk1s ·r − ×

×

×

[T 1 e1 m n+1 n bηmn mn

j

( s ) + jTn2 e2bηmn nmn ( s )]

( η s )Tη2 ejk1z H = − j

2π n0 jk1s⊥ ·r⊥ jk1z (H−z )  ejbKz H e e k1 k1z bKz + k1z Tn1 e1bηmn mmn n+1

j

( sR ) + jTn2 e2bηmn nmn ( sR )

 ( × [mmn ( s )  m−mn ( s b )] · η sb )Tn1 x1bηn



mn

 1 mn

b=±

Tn1 e1bηmn

n+1

j

+ jTn2 e2bηmn

 m ( mn s ) bKz − k1z

 n ( mn s ) bKz − k1z

−

−

mmn ( sR ) bKz + k1z

so that, by setting z = H in Eq. (139), using Eq. (142) and the result E2cη1 (r ) = 0, we obtain



( Ecη ( r⊥ , z = H ) = η s )Tη2 ejk1s⊥ ·r⊥ ejk1z H = ETη (r⊥ , z = H ). (143)

!

 nmn ( sR ) . bKz + k1z

(138)

As before, the field E1cη1 (r ) summing the contributions of the first two terms on the right-hand side of Eq. (138), as well as the field E2cη1 (r ) given by the third term on the right-hand side of Eq. (138) vanish. Hence, by taking into account the representations (24) and (25) for the coefficients e1bηmn and e2bηmn , respectively, we obtain

Ecη (r ) = ERη (r ) + E0η (r ) at z = 0

mn

( + [nmn ( s )  n−mn ( sb )] · η sb )Tn2 x2bηn

and

Ecη (r ) = ETη (r ) at z = H.



where

Aη =

b=±

 ( × sR )  m−mn ( sb )] · η sb )Tn1 x1bηn [mmn ( mn





(139)

We are now concerned with the computation of the field Ecη (r ) at the lower and upper boundaries, that is, at z = 0 and z = H. For this purpose, we extend the representations (57) and (73) for the coherent fields reflected and transmitted by the layer in the critical domains. 1. Case z = 0. Using the expression of the coherent reflected field ERη (r ) as given by Eq. (57) in conjunction with Eq. (58), and the result E2cη1 (r ) = 0, where E2cη1 (r ) is given by Eq. (97), we find the identity

8π n0  1 k1 k1z bKz + k1z 2

Cη =



mn

(150)

( Ecη ( r ) = [η sT )Tη1 + Tη1 z]ejk1s⊥ ·r ejKz z ( + [η sTR )R1η + R1η z]ejk1s⊥ ·r e−jKz z ,

(151)

where

Tϕ1 = R1ϕ = 0

(152)

1 1 − e2jKz H







1 + R0θ − Tθ2 ej(Kz +k1z )H sin θ0





− 1 − R0θ − Tθ2 ej(Kz +k1z )H sin θT , (140)

Setting z = 0 in Eq. (139), using Eq. (140) and the fact that the field E1cη1 (r ) as given by Eq. (96) vanishes, we obtain

( +η sR )Rη e

0 jk1s⊥ ·r⊥

= E0η ( r⊥ , z = 0 ) + ERη ( r⊥ , z = 0 ).

(149)

and then as



( + [nmn ( sR )  n−mn ( sb )] · η sb )Tn2 x2bηn .

(148)

e2jKz H R0 , 1 − e2jKz H η

Tθ1 = −

 ( sR )  m−mn ( sb )] · η sb )Tn1 x1bηn [mmn (



1 e2jKz H − Tη2 ej(k1z +Kz )H , 1 − e2jKz H

and

b=±

(147)

1 R0 , 1 − e2jKz H η

Dη = −

(141)

(146)



1 1 − Tη2 ej(k1z +Kz )H , 1 − e2jKz H

Bη = −

( + [nmn ( sR )  n−mn ( sb )] · η sb )Tn2 x2bηn .

( Ecη ( r⊥ , z = 0 ) = η s )e

(145)

( +η sR )ejk1s⊥ ·r⊥ (Cη ejKz z + Dη e−jKz z ),

8π 2 n0  ejKb ·r +j k1 k1z bKz + k1z

jk1s⊥ ·r⊥

(144)

( Ecη ( r ) = η s )ejk1s⊥ ·r⊥ (Aη ejKz z + Bη e−jKz z )

b=±

 ( × s )  m−mn ( sb )] · η sb )Tn1 x1bηn [mmn (

×

Thus, from Eqs. (141) and (143) it is readily seen that the coherent field Ecη (r ) satisfies the continuity conditions for the electric fields

From Eqs. (139), (141), and (143), it is straightforward to show that Ecη (r ) can first be written as

8π 2 n0  ejKb ·r Ecη ( r ) = − j k1 k1z bKz − k1z

( η sR )R0η = j



( + [nmn ( s )  n−mn ( s b )] · η sb )Tn2 x2bηn , (142)

2π n0 jk1s⊥ ·r⊥  jbKz z e e k1 k1z

8π 2 n0  ejbKz H k1 k1z bKz − k1z b=±

b=±



 1 mn

+

2. Case z = H. Using the expression of the coherent transmitted field ETη (r ) as given by Eqs. (73) and (74), and the result E1cη1 (r ) = 0, we find the identity

b=±

 1 mn

+

2π n0 jk1s ·r  1 e k1 k1z bKz − k1z

R1θ =

 (1 + R0θ )e2jKz H − Tθ2 ej(Kz +k1z )H sin θ0    + (1 − R0θ )e2jKz H − Tθ2 ej(Kz +k1z )H sin θT . 1 1 − e2jKz H

(153)



(154)

From the above analysis the following conclusions can be drawn.

A. Doicu and M.I. Mishchenko / Journal of Quantitative Spectroscopy & Radiative Transfer 230 (2019) 86–105

99

Fig. 7. The ratio |Tθ1 |/|Tθ1 | as a function of the volume concentration f for three values of the size parameter k1 a. The layer thickness is H = 50a, the incidence angle is θ0 = 30◦ , the wavenumber of the background medium is k1 = 10 μm−1 , the relative refractive index of the particles is m = 1.33, and the maximum expansion order is Nrank = 15. The ratio |R1θ |/|R1θ | has a similar dependency on f.

1. The electric fields are continuous across the boundaries (not only their tangential components). 2. For normal incidence, as well as for a ϕ -polarized incidence, the representations (128) and (151) for the coherent field coincide, i.e., the coherent field is a superposition of plane electromagnetic waves with the wavenumber K. 3. For an oblique θ -polarized incidence, the coherent field satisfies the vector Helmholtz equation but is not divergence free, i.e., the coherent field is not a superposition of plane electromagnetic waves with a wavenumber K. However, as it can be infered from Fig. 7, Tθ1 and R1θ are small compared to Tθ1 and R1θ , respectively; in the case k1 a = 5, we have |Tθ1 | < 2 × 10−3 |Tθ1 | for f = 0.01 and |Tθ1 | < 0.01|Tθ1 | for f = 0.05, while in the case k1 a = 10, we have |Tθ1 | < 1.5 × 10−3 |Tθ1 | for f ≤ 0.01 and |Tθ1 | < 8 × 10−3 |Tθ1 | for f = 0.05. Thus, for small volume concentrations, the coherent field can be assumed to be approximately divergence free.

while for a ϕ -polarized incidence, Eqs. (35) and (36) become

1= −j −

0=



Mn (−1 ) = Nn (−1 ) = −(−1 )n ωn ,

(155)

we find that for a θ -polarized incidence, Eqs. (33) and (34) of the generalized Ewald–Oseen extinction theorem become

 1 π n0  1= −j 2 ( 2n + 1 ) (T 1 x1 + Tn2 x2+θ n ) K − k1 n +θ n k1 n

0=

K + k1



 (Tn1 x1−θ n − Tn2 x2−θ n ) ,

( 2n + 1 )

n −jKH

 (−1 )n ejKH K + k1

(156)

−



e + (T 1 x1 + Tn2 x2−θ n ) , K − k 1 n −θ n

(157)



1 (T 1 x1 + Tn2 x2+ϕ n ) K − k1 n +ϕ n

 (Tn1 x1−ϕ n − Tn2 x2−ϕ n ) ,

( 2n + 1 )

 (−1 )n ejKH K + k1

n

(158)

(Tn1 x1+ϕ n − Tn2 x2+ϕ n )



e−jKH (T 1 x1 + Tn2 x2−ϕ n ) . K − k 1 n −ϕ n

(159)

Substituting Eq. (155) in Eqs. (63) and (64) gives

R0θ = − j ×



π n0

k21



( K + k1 )

1 − e2jKH



(2n + 1 )(−1 )n (Tn1 x1+θ n − Tn2 x2+θ n ),



π n0

R0ϕ = j

1 − e2jKH

(160)



k21 (K + k1 )  × (2n + 1 )(−1 )n (Tn1 x1+ϕ n − Tn2 x2+ϕ n ),

(161)

n

while substituting Eqs. (121) and (155) in Eqs. (77) and (78) gives

Tθ2 = − j +

π n0  k21

(−1 )n K + k1

Tϕ2 = − j −

(Tn1 x1+θ n − Tn2 x2+θ n )

( 2n + 1 )

n

K + k1

n

Let us particularize the above findings to the case of normal incidence. Using the special values of the functions Mn and Nn at x = 1 as given by Eq. (121), as well as their special values at x = −1 given by

(−1 )n

k21

(−1 )n

3.4. Normal incidence

+

π n0 

(−1 )n K + k1



n

1 ej(K−k1 )H (Tn1 x1+θ n + Tn2 x2+θ n ) K − k1



e−j(K+k1 )H (Tn1 x1−θ n − Tn2 x2−θ n ) ,

π n0  k21

( 2n + 1 )

n

( 2n + 1 )



(162)

1 ej(K−k1 )H (Tn1 x1+ϕ n + Tn2 x2+ϕ n ) K − k1



e−j(Kz +k1z )H (Tn1 x1−ϕ n − Tn2 x2−ϕ n ) .

(163)

The reflection and transmission coefficients for the coherent field inside the layer are then obtained by inserting Eqs. (160)–(163) in Eqs. (129)–(132).

100

A. Doicu and M.I. Mishchenko / Journal of Quantitative Spectroscopy & Radiative Transfer 230 (2019) 86–105

Note that in view of the result |x1−,η2n |  |x1+,η2n |, Eqs. (156)–(159) yield

x1+,θ2n

=

x1+,ϕ2n

=

x1+,n2 .

and R1θ = R1ϕ = 0; and furthermore,

Hence, as in the effective field approxi-

( Ecη ( r ) = η sT )Tη1 ejKsT ·r .

mation, we have R0ϕ = −R0θ and Tϕ2 = Tθ2 , as well as R1ϕ = −R1θ and Tϕ1 = Tθ1 .

Note that the coherent field computed under the sparsemedium approximation for the integration domain is (cf. Eq. (151))

3.5. Semi-infinite discrete random medium In the limit H → ∞, the above formalism reduces to that of a semi-infinite discrete random medium. Actually, because Kz > 0, we have ejKz H e+ → 0 as H → ∞. Employing this result in Eq. (15) and taking into account the representations of the elements of the matrix Jb2zH given by Eqs. (133) and (134) of Ref. [1], we find that e−jKz H e− → 0 as H → ∞. Consequently, e− → 0 as H → ∞, and the equations of the Lorenz–Lorentz law (13) and the generalized Ewald–Oseen extinction theorem (14) become + e+ = n 0 ( J+ 1a + J2a )Te+

(164)

and e0 + n 0 J+ 2z0 Te+ = 0,

( Ecη ( r ) = [η sT )Tη1 + Tη1 z]ejKsT ·r , Tθ1 = (1 − R0θ ) sin θT − (1 + R0θ ) sin θ0 .

extinction

In principle, expressions for the effective wavenumber K and the coherent field Ec (r ) can be obtained by particularizing the results for a semi-infinite discrete random medium in the case of low volume concentration f. The procedure is as follows. 1. In Eq. (168), we approximate x1+,n2 ≈ 1 for f  1, and find that the effective wavenumber K can be computed as

K = k1 − j

theorem

(166)

for a θ -polarized incidence, and

1=j

  8π 2 χn n(n + 1 ) k1 k1z (Kz − k1z ) n

× [Nn ( s · sT )Tn1 x1+ϕ n − Mn ( s · sT )Tn2 x2+ϕ n ]

(167)

π n0  k21

(2n + 1 )(Tn1 x1+n + Tn2 x2+n ).

(168)

2. For the coherent field reflected by the layer (cf. Eq. (57)),

mn

},

(169)

the reflection coefficients are given by (cf. Eqs. (63) and (64))

  8π 2 n0 Rθ = j χn n(n + 1 ) k1 k1z (Kz + k1z ) n 0

and

R0ϕ = − j ×



  8π 2 n0 χn n(n + 1 ) k1 k1z (Kz + k1z ) n

(

)

Nn  sR ·  sT Tn1 x1+ϕ n

−

(

)

Mn  sR ·  sT Tn2 x2+ϕ n

(170)



.

cos θ0 , Tϕ1 = 1 + R0ϕ , cos θT

(177)

π n0 k21

 1 ej(K−k1 )H (2n + 1 )(Tn1 x1+n + Tn2 x2+n ), K − k1 n (178)

so that from Eq. (176) and the approximation 2

Tθ = Tϕ = e

j(K−k1 )H

.

(171)

(172)

x1+,n2

≈ 1, we find

(179) − K

Hence, by using the result K = jn0Cext , where Cext is the extinction cross section of a spherical particle, we are led to the celebrated Bouguer exponential attenuation law [13] for the transmissivity T of the layer [11]:  T = (|Tθ2 | )2 = (|Tϕ2 | )2 = ej(K −K )H = e−n0Cext H .

3. Setting x1−,η2n = 0 in Eqs. (77) and (78), we find in the limit H → ∞ that Tθ2 = Tϕ2 = 0, that is, the coherent field transmitted by the layer vanishes. Consequently,

Tθ1 = (1 − R0θ )

Tθ2 = Tϕ2 = −j

2

× [Mn ( sR ·  sT )Tn1 x1+θ n − Nn ( sR ·  sT )Tn2 x2+θ n ]

(176)

It is interesting to mention that in the case of a layer of finite geometrical thickness H and normal incidence, the assumption |x1−,η2n |  |x1+,η2n | implies

8π 2 n0 ejk1R ·r k1 k1z (Kz + k1z )  ( × {[mmn ( sR )  m−mn ( s T )] · η sT )Tn1 x1+ηn

ERη ( r ) = j

( ) + [nmn ( sR )  n−mn ( s T )] · η

(2n + 1 )(Tn1 + Tn2 ).

n

( ( where E 0 ( sT ) = E0θ θ sT ) + E0ϕ ϕ sT ). Under the same approximations, Eq. (175) yields Tθ1 = 0, and the coherent field computed under the sparse-medium approximation for the integration domain is as in Eq. (177). Thus, the inconsistency between the two states of polarization of this field disappears.

n

 sT Tn2 x2+ηn

k21

Ec (r ) = ejKsT ·r E 0 ( sT ),

for a ϕ -polarized incidence. For normal incidence, x1+,θ2n = x1+,ϕ2n = x1+,n2 , and the above equations reduces to (cf. Eqs. (156) and (158))

K = k1 − j

π n0 

In view of Eqs. (18) and (19) and Eqs. (24) and (25), these approximations imply e1+ηmn ≈ e10ηmn and e2+ηmn ≈ e20ηmn ; this means that for sparse media, the conditional configurationaveraged exciting field can be considered to be approximately equal to the incident field [6]. 2. In Eq. (172), we approximate θT ≈ θ0 and R0η ≈ 0; then Tη1 ≈ ( 1, and so, from Eq. (173), Ecη (r ) = η sT ) exp(jK sT · r ). This result implies that in the case of an incident wave of arbitrary polarization, we have

  8π 2 1= −j χn n(n + 1 ) k1 k1z (Kz − k1z ) n × [Mn ( s · sT )Tn1 x1+θ n − Nn ( s · sT )Tn2 x2+θ n ]

(175)

4. Coherent field in a sparse medium

respectively. Further simplifications are listed below. Ewald–Oseen

(174)

with Tϕ1 = 0 and

(165)

1. The generalized (165) gives

(173)

(180)

In the following we intend to derive a self-contained theory for a sparse medium by starting from the initial formulae (11) and (44). Since the volume concentration is low, we (i) assume that the positions of the particles are uncorrelated (g(Ri j ) = 1), and (ii) use the far-field approximation for the fields. Note that these assumptions are standard in deriving the vector radiative transfer equation. With these simplifications, the solution method is as follows.

A. Doicu and M.I. Mishchenko / Journal of Quantitative Spectroscopy & Radiative Transfer 230 (2019) 86–105

101

1. We consider the integral equation (11) with g(Ri j ) = 1 and the sparse-medium approximation for the integration domain, that is,



 

ei = e0i + n0 i

 

D

Q ( k1 Ri j ) e j d R j , 3

(181)

j

where the incident field coefficients are given by 

e0i = ejk1 s·Ri e0 ,

(182)

e0 = 4π x ( s ) · E 0 ( s ) = 4π [xθ ( s )E0θ + xϕ ( s )E0ϕ ]

(183)

 ( with E 0 ( s ) = E0θ  θ (s ) + E0ϕ ϕ s ). Employing the far-field representation for the matrix Q(k1 Ri j ) (cf. Eq. (45) of Ref. [1]) Q ( k1 Ri j ) =

1 g0 (Ri j )Q∞ ( Ri j ), k1

(184)

Fig. 8. Scattering geometry for computing the coherent field in the case of a sparse discrete random medium.

where g0 (Ri j ) = exp(jk1 Ri j )/Ri j and (cf. Eq. (46) of Ref. [1]) Q∞ ( Ri j ) = − 4π jx ( Ri j ) · xT ( Ri j )T

= − 4π j[xθ ( Ri j )xTθ ( Ri j ) + xϕ ( Ri j )xTϕ ( R i j )] T , 

(185)

we compute the conditional configuration-averaged exciting field coefficients from the integral equation

 

ei = e0i + i

n0 k1



 

D

3 g0 (Ri j )Q∞ ( Ri j ) e j d R j .

(186)

j

2. In the integral representation (cf. Eq. (44) with the sparsemedium approximation for the integration domain)

Ec ( r ) = E0 ( r ) + n0



 

D

XT3 (k1 ri )T ei d Ri , 3

(187)

i

we use the far-field approximation for the radiating vector spherical wave functions (cf. Eq. (11) of Ref. [1])

X3 ( k1 ri ) = −

Our objective is to compute each integral in Eq. (190) and then to sum up the resulting series. Consider the first integral in Eq. (190) and denote it by I1 (Ri ). Making the change of variable (cf. Eq. (189)) R j = Ri + R ji and using the asymptotic spherical-wave expansion of the plane wave exp(jk1 s · R ji ) [14],

ejks·R ji = j

 2π

I1 ( Ri ) = j

 2π

 

g0 (ri )x ( ri )T ei d Ri . T

D

3

Ri = r + p, R j = Ri + R ji , . . .

(189)

The integration domain is D without any exclusion sphere. As shown in Fig. 8, for direction  p, p ranges from zero at the observation point to the corresponding value at the point C (where the straight line with the direction vector  p crosses the lower plane boundary); for direction  R ji , Rji ranges from zero at the origin of particle i to the corresponding value at the point Ci (where the straight line with the direction vector  R ji crosses the lower plane boundary), etc. 4.1. Conditional configuration-averaged exciting field coefficients

ei = e0i + i

+

n20 k21

n0 k1  D

D

k21



g0 (Ri j )g0 (R jk )Q∞ ( Ri j )Q∞ ( R jk )e0k d Rk d R j + · · · . 3

3

(190)

 n0 e



s 0

jk1 s·Ri



e2jk1 R ji dR ji



s 0

dR ji



sQ∞ ( s) −



e2jk1 s − 1 Q∞ (− s ) e0 . 2jk1

(192)

i

where RAi and RA are the points where the straight lines with the i direction vectors − s and  s going through Oi cross the lower and the upper boundary of the layer, respectively. By assuming that the particles are separated by distances much larger than the wavelength, i.e., k1 s  1 (except for points in the immediate vicinity of the boundary), and taking account of

" 2jk1 s " − 1" "e " 2j " =

"

"

1 " 2jk1 s " 1 − 1" ≤ "e 2 2

" "  " 2jk1 s " "e " + 1 = 1,

(193)

we find

s

"

"

1 " e2jk1 s − 1 " ≥" ". k1 2jk1

(194)

Thus, the second term on the right-hand side of Eq. (192) is much smaller than the first term, and the result is I1 ( Ri ) = j

3 g0 (Ri j )Q∞ ( Ri j )e0 j d R j

(191)

Here, s = s(Ri , − s) =  s · (Ri − RAi ) and s = s(Ri , s) =  s · ( RA − Ri ),

 2π

The iterated solution of Eq. (186) is



= j

(188)

i

The procedure of configuration averaging requires the computation of integrals over particle positions. To integrate over all positions of particle i we use a local coordinate system with the origin at the observation point; to integrate over all positions of particle j we use a local coordinate system with the origin at particle i; and so on. In other words, we make the changes of variables

 

k21



n0 ejk1s·Ri Q∞ ( s )e0

− Q∞ (− s )e0

and compute the coherent field according to

j Ec ( r ) = E0 ( r ) − n0 k1

  δ ( s+ R ji )e−jk1 R ji − δ ( s− R ji )ejk1 R ji

as k1 Rji → ∞, where δ ( s ) is the solid-angle delta function, we obtain

j g0 (ri )x( ri ), k1 

2π k1 R ji

k1



n0 s ejk1s·Ri

1 k1



Q∞ ( s ) e0 .

(195)

For the second integral, denoted by I2 (Ri ), we change variables according to R j = Ri + R ji and Rk = R j + Rk j , employ the asymptotic spherical-waves representations for the plane waves exp(jk1 s · R ji ) and exp(jk1 s · Rk j ), and hence derive

102

A. Doicu and M.I. Mishchenko / Journal of Quantitative Spectroscopy & Radiative Transfer 230 (2019) 86–105

n20

ejk1s·Ri 2

I2 ( Ri ) =

k1 ×



D

 2π k1

g0 ( Rk j ) e

2 n0



ejk1s·Ri

1 2π = j n0 s 2 k1

Using the incident field representation

g0 (R ji )ejk1s·R ji Q∞ (− R ji ) jk1 s·Rk j

D

= j



2 e

E0 (r ) = ejk1 s(r,−s) E0 (rA ),

Q∞ (− Rk j )e0 d Rk j d R ji

1 k1

jk1 s·Ri



3

Q∞ ( s)

1 k1

2  s 

s

e0

Q∞ ( s)

2

0

R ji





Ec (r ) = exp j k1 +

dRk j dR ji

e0 .

(196)

Next, we use Eqs. (183) and (185), the expression of the elements of the amplitude matrix given by (cf. Eq. (13) of Ref. [1])

Sημ ( s,  s) = −

4π j T x ( s )Txμ ( s ), k1 η

η, μ = θ , ϕ ,

(197)

and the well-known identity for spherically symmetric particles

η  = μ,

Sημ ( s,  s ) ≡ 0,

(198)

to derive

(199)

where

j  S(0 ) = Sηη ( s,  s) = − (2n + 1 )(Tn1 + Tn2 ), 2k1 n

η = θ , ϕ . (200)

k1

Q∞ ( s)

n

e0 = S n ( 0 )e0 , n = 1, 2, . . . ,

(201)

Proceeding similarly for all terms in the series (190), and accounting for Eq. (201), we find that the sum of the series (190) is



 

ei = exp j k1 s · Ri + i

2π n0 S(0 )s(Ri , − s) k1



e0 .

(202)

To compute the coherent field, we consider the integral representation (188). Using the series expansion (cf. Eq. (202)) ei = e

jk1 s·Ri

i





2π 1+j n0 S(0 )s(Ri , − s ) + · · · e0 . k1

(203)

D

K = k1 +

2π π n0  n0 S ( 0 ) = k1 − j 2 (2n + 1 )(Tn1 + Tn2 ). k1 k1 n

(210)

and so, the following equivalent representation for the coherent field

K0 = k1 s+

 2π z n0 S ( 0 ) , k1 cos θ0

2π sn+1 (r, − s ) jk1s·r T =j e x ( s )Te0 , n = 0, 1, . . . , k1 n+1

K0 = k1 s + ( K − k1 )

 z . cos θ0

j T 4jπ T x ( s )Te0 = − x ( s )Tx ( s ) · E 0 ( s) k1 k1 4jπ T  ( = − [x ( s )θ ( s ) + xTϕ ( s )ϕ s )]T[xθ ( s )E0θ + xϕ ( s )E0ϕ ] k1 θ = S(0 )E 0 ( s ),

(204)

−

 

k1

(214)

(215)

respectively. Thus, for a layer with a sparse concentration of particles, the effective wavenumber is computed from Eq. (209), the conditional configuration-averaged exciting field coefficients from Eq. (214), and the coherent field from Eq. (215). Some comments are in order. 1. The effective wave vector K, defined by K = K sT , is related to K0 by the relation

(205)

K0 = K + Kz zˆ ,

(216)

where



n0 S(0 )s(r, − s ) E0 (r ).

(213)

Using the relation s(r, − s) =  z · r/ cos θ0 , we express Eqs. (202) and (208) in terms of the wave vector K0 as (compare with Eq. (12))

Ec (r ) = ejK0 ·r E 0 ( s ),

and the identity (cf. Eqs. (183), (197), and (198))

Ec (r ) = exp j

(212)

and

3 g0 (ri )ejk1s·Ri sn (Ri , − s )xT ( ri )Te0 d Ri

 2π

(211)

Thus, in contrast to a layer with densely packed spherical particles, the coherent field in a layer with a sparse concentration of particles is an upwelling wave. The upper boundary does not produce a downwelling wave, and the upwelling coherent field is uniquely determined by the boundary condition (210). It is for this reason that in the introductory part of the section, we computed the effective wavenumber and the coherent field by particularizing the results for a semi-infinite discrete random medium. Let us define the wave vector K0 by

i

we end up with

(209)

Then, from Eq. (208) we deduce that the coherent field is analytically equivalent to a plane electromagnetic wave with the wavenumber K, propagation direction  s, and amplitude exp(−jK s· rA )E0 (rA ). Therefore, K can be interpreted as the effective wavenumber in the particulate medium. Setting r = rA in Eq. (208) and taking into account that s(rA , − s ) = 0, we find the boundary condition

ei = ejK0 ·Ri e0

the result



(208)

implying that (cf. Eq. (209))

4.2. Coherent field

 



Let us define the wavenumber K by the relation (compare with Eq. (176)) [16]

Ec (r ) = ejKs(r,−s) Ec (rA ).

From Eq. (199), it is readily seen that in general,

1



2π n0 S(0 ) s(r, − s ) E0 ( rA ). k1

Ec ( rA ) = E0 ( rA ),

Q∞ ( s )e0 = 4π k1 [xθ ( s )Sθ θ ( s,  s )E0θ + xϕ ( s )Sϕϕ ( s,  s )E0ϕ ]

= k1 S ( 0 )e0 ,

(207)

with E0 (rA ) = exp(jk1 s · rA )E 0 ( s ), we express Eq. (206) as

3

(206)

Here, s(r, − s ) is defined through the relation s(r, − s) =  s · ( r − rA ), where rA is the point where the straight line parallel to the incidence direction and going through the observation point crosses the boundary of the medium (see Fig. 8).

Kz =

1 − cos(θT − θ0 ) K. cos θ0

(217)

Approximating θT ≈ θ0 , yields Kz ≈ 0, and so, K ≈ K0 . Therefore, K0 can be interpreted as the effective wave vector in a sparse medium. In view of the approximations K ≈ K0 and θT ≈ θ0 implying sT ≈  s, we deduce that the representations (177) and (215) for the coherent field are equivalent.

A. Doicu and M.I. Mishchenko / Journal of Quantitative Spectroscopy & Radiative Transfer 230 (2019) 86–105

2. Setting r = Ri in the representation of the coherent field as given by Eq. (215) yields Ec (Ri ) = exp(jK0 · Ri )E 0 ( s ), so that from Eq. (214) and the representation (cf. Eq. (183)) e0 = 4π x ( s ) · E 0 ( s ), we find

 

ei = 4π x ( s ) · Ec ( Ri ).

Inserting Eq. (218) in Eq. (188) and taking into account the representation for the far-field scattering dyadic (cf. Eq. (15) of Ref. [1])

4π j T x ( r )Tx ( s ), k1

we infer that for sparse media, the coherent field satisfies the integral equation

Ec ( r ) = E0 ( r ) + n0



g0 (ri )A( ri ,  s ) · Ec ( Ri ) d Ri . 3

(219)

In the following parts of this series, we will compute the second-order moment of the fields by approximating the conditional configuration-averaged exciting field coefficients of a layer with densely packed particles by those of a semi-infinite medium with densely packed particles. In this context, with K being calculated from the generalized Lorenz–Lorentz law for a dense semiinfinite discrete random medium at normal incidence, Tishkovets and Jockers [15] assumed that the sparse-medium approximation (214) is also valid in this case. The benefit of using Eq. (214) to gether with Eq. (213) is that the vector e+ , which determines ei i according to the generalized Lorenz–Lorentz law and the generalized Ewald–Oseen extinction theorem for a dense semi-infinite discrete random medium, does not need to be computed.

In this section, we will show that the Twersky and Foldy approximations are implicit in our derivation. Consider a configuration with N − 1 particles iN−1 in which particle i is removed from the group. According to the Twersky approximation, the field exciting particle i at a point ri near particle i is the total electric field that would exists at that point if the particle i were removed from the group, i.e.,

(220)

Here, the superscript “(i)” indicates that the fields are written in the coordinate system of particle i. Taking the configuration average of Eq. (220) with the position of particle i held fixed, and using the configuration average rule



f (r, N ) =



i

yields



(i )



(i )

f (r, N ) p(iN−1 |Ri ) diN−1

Eexci (ri ) = Ec (ri | i

i N−1

).

(221)







 

i

i

(i ) Eexci (r ) = Eexc (r ) = XT1 (k1 ri ) ei . i i i



 

D

XT3 (k1 r j )T e j d R j . 3

(223)

j

Choosing ri in the neighborhood of particle i before the removal, we use the translation addition theorem X3 (k1 r j ) = T31 (k1 Ri j )X1 (k1 ri ) with r j = ri + Ri j , the incident field representation E0 (r ) = E0(i ) (ri ) = exp(jk1 s · ri )E0 (Ri ) = XT1 (k1 ri )e0i , and the integral equation for the conditional configuration-averaged exciting field coefficients (cf. Eq. (181)),



 

ei = e0i + n0

D

 

T T31 ( k1 Ri j )T e j d R j , 3

j

to obtain

 

Ec (r|iN−1 ) = Ec(i ) (ri |iN−1 ) = XT1 (k1 ri ) ei .

(224)

i

From Eqs. (222) and (224), the Twersky approximation (221) readily follows. The Foldy approximation states that the coherent field at a point ri near particle i is the conditional configuration average of the field exciting particle i at that point, or equivalently and in view of Eq. (221), that the coherent field at a point ri near particle i is the coherent field that would exist at that point if the particle i were removed from the group:





(i ) Ec(i ) (ri ) = Eexc (r ) = Ec(i) (ri |iN−1 ). i i

(225)

i

 

To prove this result we use Eq. (218), that is, ei = 4π x ( s) · i Ec (Ri ), where (cf. Eqs. (213) and (215))

Ec (Ri ) = ej(K−k1 )s(Ri ,−s) E0 (Ri ),



(i )



(226)





(i ) Eexc (r ) = ejk1s·ri Ec (Ri ) = ejk1s·ri ej(K−k1 )s(Ri ,−s) E0 (Ri ). i i i

(227)

Apart from that, Eqs. (213) and (215) also give

Ec (r ) = Ec(i ) (ri ) = ejk1s·ri ej(K−k1 )s(r,−s) E0 (Ri ).

(228)

If the point ri is near particle i, and the boundary is in the farfield region of this particle, we approximate s(Ri , − s ) ≈ s(r, − s ); hence from Eqs. (227) and (228), we obtain the Foldy approximation (225). It is interesting to note that the representation Ec(i ) (ri ) = exp(jk1 s · ri )Ec (Ri ), which follows from Eq. (227), shows that the coherent field near a particle can be analytically approximated by a plane electromagnetic wave with wavenumber k1 and propagation direction  s. 5. Discussion

Thus, in the Twersky approximation, the conditional configurationaveraged field exciting particle i at a point ri near particle i is the coherent field that would exists at that point if particle i were removed from the group. Let us prove this conjecture. The conditional configuration average of the field exciting particle i at the field point r = Ri + ri is



= E0 ( r ) + n0

and Eq. (222) to conclude that Eexci (ri ) is analytically equivalent i to the plane electromagnetic wave

4.3. The Twersky and Foldy approximations

(i ) Eexc (r ) = E(i) (ri |iN−1 ). i i

  Esct j (r ) j=i

i

This equation is the so-called Foldy integral equation for the coherent field.



Ec (r|iN−1 ) = E0 (r ) +

(218)

i

A( r,  s) = −

(N − 1 )/V )

103

(222)

On the other hand, upon removing particle i from the group, the coherent field at the same field point r produced by the remaining N − 1 particles, denoted by Ec (r|iN−1 ), computes as (n0 ≈

The analysis performed in this paper is based on the assumption that the conditional configuration-averaged exciting field co  efficients are of the form ei (Ri ) = b=± exp(jKb · Ri )eb , and is rei stricted to the computation of the zeroth-order fields without a special treatment of the critical domains. In this setting we have calculated the coherent fields reflected and transmitted by the layer and the coherent field inside the layer, and found that these fields are analytically equivalent to plane electromagnetic waves. The main problem which arises is that if the representations of the zeroth-order fields are extended in the critical domains, the boundary conditions for the electric fields at the layer interfaces are not satisfied. A rigorous approach dealing with this problem is to compute higher-order approximations to the fields (which are valid in the critical domains), while a pragmatic and approximate approach is to compute the coherent

104

A. Doicu and M.I. Mishchenko / Journal of Quantitative Spectroscopy & Radiative Transfer 230 (2019) 86–105

field inside the layer from the boundary conditions for the electric fields. Another approximate method relies on the computation of the coherent field inside the layer under the sparse-medium approximation for the integration domain. In this case, we have found that (i) the boundary conditions for the electric fields are satisfied; (ii) for normal incidence and an oblique ϕ -polarized incidence, the coherent field is a superposition of plane electromagnetic waves; while (iii) for an oblique θ -polarized incidence, the coherent field is not divergence free. The inconsistency between the two states of polarization disappears for sparse media. This approximate method will be used as the starting point in a forthcoming paper to derive a simplified radiative transfer equation for dense media. In Ref. [1] we discussed an approach proposed by Kristensson [12] to solve the integral equation (11) without assuming a special form for the conditional configuration-averaged exciting field coefficients. For an oblique η-polarized incidence, we assumed that the solution is of the form

  e iη

(Ri ) = ejk1s⊥ ·Ri⊥ eη (zi ), 0 ≤ zi ≤ H,

i

(229)

and found that eη (zi ) satisfies the integral equation



eη ( z i ) = e

jk1z zi

×e

e0η + n0

jk1s⊥ ·R ji⊥

D−D2a (Ri )

eη (z j )g(R ji ) d R j

with k1z = k1z (k1s⊥ ) =



(230)

k21 − k21s⊥ and k1s⊥ = k1 sin θ0 . After solv-

ing the integral Eq. (230) for eη (zi ) = [e1ηmn (zi ), e2ηmn (zi )]T (eventually by using the approach described in Ref. [1]) the coherent field reflected by the layer is computed as

ERη ( r ) = n0



D

mn

1 Tn1 eηmn

( )

zi M3mn

( k1 ri )



+ Tn2 eηmn (zi )N3mn (k1 ri ) ejk1s⊥ ·Ri⊥ d Ri , 2

3

(231)

where the integrals over the particle positions are given by

 D

=

2π n

j k1 k1z 

ejk1R ·r mmn ( sR )

d Ri



H 0

ejk1z zi eηmn (zi )dzi , 1

(232)

N3mn (k1 ri )eηmn (zi )ejk1s⊥ ·Ri⊥ d Ri 3

2

2π n

j k1 k1z

ejk1R ·r jnmn ( sR )



H 0

ejk1z zi eηmn (zi )dzi . 2

(233)

As in the calculation of the integrals (55) and (56), the integrals (232) and (233) are evaluated by using the integral representations (52) and (53). Similarly, the coherent field transmitted by the layer is computed as

ETη ( r ) = E0 ( r ) + n0





mn

+

2 Tn2 eηmn

( )

zi N3mn

1

( k1 ri ) e

jk1s⊥ ·Ri⊥ 3

d Ri ,

(234)

=

M3mn (k1 ri )eηmn (zi )ejk1s⊥ ·Ri⊥ d Ri 3

1

D

2π n

j k1 k1z

ejk1s ·r mmn ( s)

 0

H

e−jk1z zi eηmn (zi )dzi , 1

n

j k1 k1z



ejk1s ·r jnmn ( s)

H

0

e−jk1z zi eηmn (zi )dzi . 2

(236)

Acknowledgements We thank two anonymous reviewers for very insightful and stimulating comments. The authors are pleased to acknowledge very fruitful discussions with Gerhard Kristensson during the completion of this work. Adrian Doicu acknowledges the financial support from DLR programmatic (S5P KTR 2 472 046) for the S5P algorithm development. Michael I. Mishchenko was supported by the NASA Radiation Science Program and the NASA Remote Sensing Theory Program.

The regular vector spherical wave functions can be expressed as integrals over vector spherical harmonics [10]

 2π  π

1

L1mn (kr, θ , ϕ ) =

n−1

4π j

M1mn (kr, θ , ϕ ) =

N1mn (kr, θ , ϕ ) =

0

0



lmn (θ  , ϕ  )ejkrr ·r sin θ  dθ  dϕ  , (237)

 2π  π

1 4π j

n

0

0



mmn (θ  , ϕ  )ejkrr ·r sin θ  dθ  dϕ  , (238)

 2π  π

1 n−1

4π j

0

0



nmn (θ  , ϕ  )ejkrr ·r sin θ  dθ  dϕ  , (239)

where (θ  , ϕ  ) are the polar angles of the direction  r . From Eqs. (237)–(239), and the orthogonality relation of the vector spherical harmonics

 2π  π 0

β

vαmn (θ , ϕ ) · v−m n (θ , ϕ ) sin θ dθ dϕ = δαβ δmm δnn , (240)

where α , β = 1, 2, 3, v1mn (θ , ϕ ) = lmn (θ , ϕ ), v2mn (θ , ϕ ) = nmn (θ , ϕ ), and v3mn (θ , ϕ ) = mmn (θ , ϕ ), we find the following expansion of the dyadic exp(jp · r )I [6]:

Iej prp·r = − 4π



n+1

j

[l−mn (θ p , ϕ p )  L1mn ( pr, θ , ϕ )

mn

+ jm−mn (θ p , ϕ p )  M1mn ( pr, θ , ϕ ) + n−mn (θ p , ϕ p )  N1mn ( pr, θ , ϕ )],

(241)

References

with



2π

where (θ p , ϕ p ) are the polar angles of the direction  p.

Tn1 eηmn (zi )M3mn (k1 ri )



3

An interesting feature of the method proposed by Kristensson is that the effective wavenumber K is computed from the equation Tη2 (K ) = Tη2EFA (K ), and not by finding the roots of a determinant as in Ref. [1].

0

D

=

(k1 ri )eηmn (zi )e

jk1s⊥ ·Ri⊥ 3

1

M3mn

=

N3mn (k1 ri )eηmn (zi )ejk1s⊥ ·Ri⊥ d Ri 2

D

Appendix A

Q(−k1 R ji ) 3



(235)

[1] Doicu A, Mishchenko MI. Electromagnetic scattering by discrete random media. I: The dispersion equation and the configuration-averaged exciting field. J Quant Spectrosc Radiat Transf, in press. [2] Waterman PC, Truell R. Multiple scattering of waves. J Math Phys 1961(2):512–37. [3] Fikioris JG, Waterman PC. Multiple scattering of waves. II. Hole corrections in the scalar case. J Math Phys 1964;5:1413–20. [4] Fikioris JG, Waterman PC. Multiple scattering of waves. III. The electromagnetic case. J Quant Spectrosc Radiat Transf 2013;123:8–16. [5] Tsang L, Kong JA. Scattering of electromagnetic waves from a half space of densely distributed dielectric scatterers. Radio Sci 1983;18:1260–72.

A. Doicu and M.I. Mishchenko / Journal of Quantitative Spectroscopy & Radiative Transfer 230 (2019) 86–105 [6] Tsang L, Kong JA. Scattering of electromagnetic waves: advanced topics. New York: Wiley; 2001. [7] Mishchenko MI. Vector radiative transfer equation for arbitrarily shaped and arbitrarily oriented particles: a microphysical derivation from statistical electromagnetics. Appl Opt 2002;41:7114–34. [8] Mishchenko MI, Travis LD, Lacis AA. Multiple scattering of light by particles: radiative transfer and coherent backscattering. Cambridge, UK: Cambridge University Press; 2006. [9] Mishchenko MI. Electromagnetic scattering by particles and particle groups: an introduction. Cambridge, UK: Cambridge University Press; 2014. https://www. giss.nasa.gov/staff/mmishchenko/publications/Book_4.pdf. [10] Boström A, Kristensson G, Ström S. Transformation properties of plane, spherical and cylindrical scalar and vector wave functions. In: Varadan VV, Lakhtakia A, Varadan VK, editors. Field representations and introduction to scattering. Amsterdam: Elsevier; 1991. p. 165–210.

105

[11] Gustavsson M, Kristensson G, Wellander N. Multiple scattering by a collection of randomly located obstacles — numerical implementation of the coherent fields. J Quant Spectrosc Radiat Transf 2016;185:95–100. [12] Kristensson G. Coherent scattering by a collection of randomly located obstacles — an alternative integral equation formulation. J Quant Spectrosc Radiat Transf 2015;164:97–108. [13] Bouguer P. Essai d’optique sur la gradation de la lumiere. Paris: Claude Jombert; 1729. [14] Saxon DS. Lectures on the scattering of light. Science Report No. 9, Department of Meteorology. Los Angeles: Unviversity of California; 1955. [15] Tishkovets VP, Jockers K. Multiple scattering of light by densely packed random media of spherical particles: dense media vector radiative transfer equation. J Quant Spectrosc Radiat Transf 2006;101:54–72. [16] Martin PA. Multiple scattering: interaction of time-harmonic waves with N obstacles. Cambridge, UK: Cambridge University Press; 2006.