MONTE CARLO PROCEDURES. A PROBLEM SCATTERING AND NON-SPHERICAL PARTICLES

FOR

MULTIPLE

FIRST REPORT ON THE WORK CARRIED OUT AT THE DEPARTMENT OF PHYSICS, UNIVERSITY OF FLORENCE, ITALY FOR THE CONTRACT INTAS 05-1000000-8024. 30June 2007 P.Bruscaglioni, S.Lolli. Department of Physics, University of Florence S.Del Bianco, IFAC CNR. Florence sect 1 . sect.2 2a. sect.3 3a. 3b. sect 4 sect.5 5a 5b sect.6

The kind of problem addressed by the research The "modes" of polarization Polydispersion of particles. Simple checks of the formalism for the modes Small dielectric ellipsoids Small spherical particles with chiral properties. Orthogonality and division of fields into modes Simple particles. Small rotational ellipsoids A case of polarization change. The field scattered in the forward direction Possibility for Monte Carlo procedures.

REFERENCES

page 2 page 3 page 6 page page page page page page page

7 11 12 15 19 20 22

page 24

Appendix 1 Propagation and scattering for the two modes page Appendix 2 Two modes with different extinction coefficients page Appendix 3. A line for a Monte Carlo code page Appendix 4 A note on the phase propagation page Appendix 5 Small dielectric ellipsoids: the scattered field. page

25 27 28 29 33

1

MONTE CARLO PROCEDURES. A PROBLEM FOR MULTIPLE SCATTERING AND NON-SPHERICAL PARTICLES P.Bruscaglioni, S.Lolli. Department of Physics, University of Florence S.Del Bianco, IFAC CNR. Florence 1) The kind of problem. The linear extinction coefficient of an otherwise isotropic non-absorbing medium containing non-spherical or chiral particles depends on the polarization of propagating radiation. Moreover for such type of suspended particle, in general polarization changes during propagation while coherent radiation is attenuated by scattering and absorption. A brief discussion of this point is shown later of this paper. This point is of interest for the lidar technique based on the Radiative Transfer scheme. This is based on the considering cross sections of particles and the corresponding linear extinction coefficient of the medium. A number of models of lidar returns relevant to this scheme, both analytic or Monte Carlo have been considered (papers by the INTAS group: Munich group. Florence group, Minsk group, Quebec group), Next section of this paper deals with the mentioned change of polarization. However, it is known [1] [2] that there are generally two modes of polarization for which polarization is maintained during the propagation. A simple way to obtain these modes can be related to the properties of the scattering matrix connecting the components of the incident and scattered radiation. A couple of examples relevant to simple shape particles and small spherical with chiral properties are presented Suggestions for taking into account the effect if a Monte Carlo code has to be developed are indicated at the end.

2

2. The "modes" of polarization. In the case of a medium containing one type of particle the relationship between the incident (i) and scattered (s) fields (transverse components in the far field) via the amplitude scattering M matrix M can be written as (*): (1) Es = M EI with:

M=

exp(ikr ) ⎛ f11 f12 ⎞ ⎜⎜ ⎟⎟ . r ⎝ f 21 f 22 ⎠

(*) In this report bold upper case characters are used for vectors and matrices, bold lower case characters for unit vectors. In general the elements fmn depend on the direction of arrival and scattering. Given the incident field propagating in the z direction:: Ei = E0 u exp (ikz), with the unit vector u = ax x + ay y defining polarization, we show that the Optical Theorem for the particle extinction cross section Ce can be written as (2)

2 2 ⎛ 4π ⎞ Ce = ⎜ ⎟ Im⎧⎨ f11 ax + f12axa y + f21axa y* + f22 a y ⎫⎬ ⎭ ⎝k⎠ ⎩

Where ax has been taken as reference for phase, and the matrix elements fmn are taken in the forward direction. For non-spherical particles C e depends on the direction of incidence on the particle. Eq.2 and the consequent linear extinction coefficient of the medium can be obtained as follows from the general relationship 3) Ce = (4π / k) Im ( f · u*) Where the scattering amplitude f is taken in the forward direction, and u is the unit vector indicating polarization. From the matricial relationship 1, one has : f = (f11 ax + f12 ay) x + (f21 ax + f22 ay) y (fmn in the forward direction). Let us define as x, y (x, y) axes of the reference system those along two directions for which the 2x2 scattering matrix J has the elements f11, f12, f21, f22, which are functions of the direction s with respect to the incident wave direction z. Let the incident wave field Ei (unitary amplitude) have the complex components a x, a y along the x and y axes respectively. (| ax|2 + | ay|2 = 1. Thus polarization is defined). Let us first consider the case of a the medium containing equal (and equally oriented) particles, with N number of particles per unit volume. Consider an elementary layer (width dz) where the scatterers are contained, and the scattered field on an x y plane at a distance D from the layer. 3

If we keeps the same reference axes as for the incident field, in a generic point on this plane the scattered field has also a component parallel to the z axis. This component, will now be neglected since we are interested is the component parallel to the x y plane. This part of the scattered field element dEs can be added to the incident wave, and is given as d Es= exp (ik (ρ2+D2)1/2) / (ρ2 + D2)1/2 ((aX f11 + aY f12) x +( aX .f21 + aY f22)y) N dx dy dz ρ = (x2 + y2)1/2 Ν number of particles per unit volume, and the parameters fmn are linear coefficients depending on the positions x,y and D. If the layer is extended laterally one can apply the principle of stationary phase to obtain by integration the field scattered by the whole layer : (4) d Es = (2 i π / κ) exp (ikD) ((ax f11 + ay f12) x +( ax .f21 + ay f22)y) Ndz where the fmn parameters are the elements of the scattering matrix taken in the forward direction. Eq.4 show that there is a uniform scattered field on the plane at D. Thus dEs can be coinsidered as a plane wave travelling in the z direction, and can be added to the original wave field. Let us consider power W per unit area at the distance D. Apart from a constant factor 1 / (2Z), by taking into account the differential, one has : (5) W = |EI + dEs|2 = |EI |2 + 2 RE{ EI . dEs*} By considering eq. 4 one has the increment per unit distance (dz = 1 ) (6) dW : = = - (4π/k ) N Im {f11 |α|2 +f12 aX ∗ aY + f21 aX aY ∗ + f22 | aY |2} dz (fmn in the forward direction).since the assumed incident power is unitary, one obtains the linear extinction coefficient of the medium This corresponds to the expression of the extinction cross section of eq, 2. 2⎫ 2 ⎛ 4π ⎞ ⎧ * = + + + C Im f a f a a f a a f a ⎜ ⎟ ⎨ 11 x e 12 x y 21 x y 22 y ⎬ (2) ⎭ ⎝ k ⎠ ⎩ As for polarization, let us still consider the case of equal scatterers. In order to obtain the two mode of polarization not changing during propagation, one can observe that, according to Eq.5, the x and y components of the scattered field in the forward direction are given as ax f11 +ay f12 and axf21 + ay f22 apart from a constant factor. One can impose that the the polarization of the scattered field in the forward direction is equal to that of the incident field. That is the ratio of complex amplitude components is the same as that of the incident field.

4

(ax f11 + ay f12) / ( ax f21 + ay f22) = ax / ay We thus obtain a second order equation for the ratio ax / ay, and the Ex and Ey component of the incident field. Apart from constant factors one has the two modes : (8) a,b E1

=

x - 2 y f21 / ( f22 - f11 + R)

E2

=

y + 2 x f12 / ( f22 - f11 + R)

with: R = ((f22-f11)2 + 4 f12 f21)1/2, and the elements fmn of the scattering matrix are taken in the forward direction.

5

2a. Polydispersion of particles. If there is a polydispersion of particles the expressions for the modes and the linear extinction coefficient of the medium are obtained from those for the monodispersions, with the following substitutions (i : index for the kind of particle):

with:

(9)

g11 , g12 , g 21 , g 22

(10)

g mn =

∑N f

i mni

i

in the place of

,

and Ni number of particles of kind i per unit volume. becomes: (11)

σe =

∑N C i

i

f 11 , f 12 , f 21 , f 22

ei

The extinction coefficient

.

(Different indexes i for particles with different particles and orientation. All the quantities depend on the direction of propagation.) And the polarization modes are obtained from Eqs. 8a, 8b by substituting the quantities gmn in place of fmn One can see that Eqs.8a,b for the polarization modes correspond to those presented in Ref.1, Sect.3.8, where they were obtained by considering extinction matrix for Stokes vectors. Our introduction of modes and linear extinction coefficient of the medium starts directly from considering the scattering matrix 1.

6

3. Some simple checks for the MODES formalism. 3a. Small dielectric ellipsoids. We consider some simple cases, in order to check some points of the relationships derived in Sect.2. We first refer to dielectric rotational ellipsoids. A particular case is now presented in details. A list of other cases is presented in Ref.3 The polarization modes are obtained for a couple of dielectric (non-chiral) ellipsoids. One ellipsoid has its mayor axis along x and the second, of equal shape, has its mayor axis in the x y plane, with an angle of 45°with the axes. We use the expressions for the elements of the amplitude matrix, which can be obtained as indicated in Sect. 5. The forms shown below correspond to those of a paper by Ablitt et al.[4], where the case of small dielectric ellipsoids is dealt with, also taking into account possible chiral properties.. For a general orientation of the ellipsoid, the elements of the matrix A in the forward direction, with the parameter α determining the elongation, are given as: f11 = 1 +( α -1) sin 2β cos2 φ f22 = 1 + ( α -1) sin2β sin2 φ f12 = ( α -1) sin2β sin φ cos φ f21 = ( α -1) sin2β sin φ cos φ, where the angles are defined in Fig.1, and α is a parameter depending on the ratio of the mayor (a) to minor axes (b) (see Sect. 5 for an explanation) ( α = 1 for a sphere).

7

z

z' θ A β y

φ p x

Fig.1 The propagation direction of the incident radiation is along z. The scattering direction is along z', and is on the plane y z (scattering plane). The vector A is along the mayor axis of the dielectric ellipsoid. It makes the angle β with the direction z. The unit vector p is along the projection of A on the plane perpendicular to z. p makes the angle φ with the axis x.

8

Then, for the case under examination, in the forward direction, one has, by summing the analogue elements of the two matrices (Sect.2a): 3 (α − 1) 2 1 g 22 = 2 + (α − 1) 2 1 . g12 = (α − 1) 2 1 g 21 = (α − 1) 2 g11 = 2 +

(12)

By using Eqs. 8a.b of Sect.1, and Eqs. 5 and 6 of Sect.3, one obtains the two vectors x / ( 205 - 1) + y

(13a)

and: . x - y / ( 205 - 1)

(13b)

They are two linear polarization states, along the bisector line between the two mayor axes (8a), and perpendicular to this direction (8b). They correspond to the unit vectors with components: Ex = cos 22.5°, Ey = sin 22.5°, Ex = cos 67.5°, Ey = - cos 22.5°.

and:

Although this conclusion seems satisfactory, we need to check if radiation with one of these polarization states does not change the relevant Stokes Vector, apart from a common attenuation factor for the four elements. That is, considering the I, Q, U, V formulation of the vector, the relative attenuation must be equal for the four elements. We take into account the extinction matrix for the Stokes Vector (Sect.1.VIII of Ref.5 - Consider the present definition of the Stokes vector elements). Let us consider Eq. 13a. For us

9

I = E x ⋅ E x ∗ + E y ⋅ E y∗ = 1

( )

Q = E x ⋅ E x ∗ − E y ⋅ E y ∗ = cos 45o =

{

}

( )

U = 2 Re E x ⋅ E y ∗ = sin 45o =

{



}

1

1 2

.

2

V = −2 Im E x ⋅ E y = 0

(14) And one obtains for a unitary increment of z:

(

)

− dI = 4 + 2 + 2 (α − 1) − dQ = − dU =

(4 + 2(α − 1)) + (α − 1) 2 (4 + 2(α − 1)) 2

+ (α − 1)

− dV = 0

The relative variations of the four Stokes parameter are the same. The same conclusion can be obtained for the 13b) polarization. Checks analogous to that of Sect.3a have been also carried out for single ellipsoids with different orientation. The results can be found in Ref.3. By summarizing the results, we show in Ref.3 that the two modes obtained by applying Eqs. 8a,b of Sect.1 are with linear polarization parallel or perpendicular to the projection on the x y plane of the direction of the ellipsoid mayor axis. No condition is found for the case where the axis is along the propagation direction z. For a case where the directions of the ellipsoid axes are distributed homogeneously within a solid angle, the two modes are with linear polarization perpendicular and parallel to a central direction of the angle. All the cases show that the Lambert Beer law with a constant extinction coefficient is not valid for a homogeneous medium containing such type of particle, apart from some special cases of polarization and particle attitude. One can see, however, that the Lambert Beer law re-assumes its validity for a random orientation of the ellipsoids. This can be easily observed by averaging fmm over a complete solid angle: no particular polarization modes result. These conclusions look rather trivial. However they can be extended to particles, or group of particles, where a symmetry allows one to choose two directions in the plane perpendicular to the incidence direction for which g12 = g21 = 0.

10

3b Small spherical particles with chiral properties. For small dielectric spheres with chiral properties the scattering matrix can be written in the form, apart from a constant factor: cosϑ iγ (1 + cosϑ )⎞ ⎛ ⎟⎟ , M = exp(ikr) / r ⎜⎜ 1 ⎝ − iγ (1 + cosϑ ) ⎠

(15)

with γ real parameter defining chirality, and ϑ scattering angle. Eqs 8ab of Sect.2 indicate the polarization vectors (apart from a normalizing factor): x+iy

and,

x-iy

(16)

they are two modes with circular polarization. The form of the Matrix is obtained by using the formalism by Bohren and Huffman's book [6]where the scattered field for a small chiral sphere is given as an expansion in spherical functions, as an extension of Mie theory. The extinction cross section and the scattering matrix for small chiral dielectric particles can be obtained by using the formalism of Ref.6. (The elements of the amplitude matrix f11, f12, f21, f22 correspond to the elements S2, S3, S4, S1 of Ref.5, Eq.3.12 divided by -ik.) For small chiral spheres ((2π / λ) a << 1, a radius) the two relative refractive indexes for the two circular polarizations are indicated as Mr and Ml, with average M. Let us write : Mr = M + Δ and Ml = M -Δ . By using the first order approximation of the spherical Bessel functions, one obtains the first order approximation for the functions W1, V1, A1, B1 of Ref.5, Sect.8.3.. One observes that the scattering cross section σs can be approximated by its zero order approximation in Δ , according to the expression of Csca in Ref.5 Sect. 8.3.1. For both circular polarizations; one has (k = 2 π / λ) : It corresponds to the value for a small non-chiral sphere with refractive index M. (having neglected the first order term in Δ). For applying a Monte Carlo procedure, the formalism indicated in Sect 6 would need the two different cross sections for the two circular polarizations. One can calculate the average cross section , and then consider Eq.2 of Sect. 1 to evaluate the relative differences between the two cross sections.

11

4. Orthogonality of the MODES, Decomposition of the incident field. 4 a. Orthogonality It has been shown how to obtain the two “Modes” propagating with maintained polarization. Given the modes (from Eqs. 8ab of Sect.2, and next Eq.4 ), one can consider the conditions for their orthogonality from the point of view of power propagation. In general, if the vectors : E1 = Ax + By and the vectors E2 = Cx + Dy , represent the modes for radiation propagating in the z direction, the summed propagating power is given as: ⏐E1⏐2+ ⏐E2⏐2 + 2 Re{E1•E2*} = ⏐E1⏐2+ ⏐E2⏐2 + 2 Re {A C* +BD*}

(17)

In order to check the modes orthogonality, the third term of Eq.2.must be zero. Re {A C* +BD*} = 0

(18)

This is obvious if the Modes represent two perpendicular linear polarizations in the plane x y. normal to the z direction of propagation of the incident radiation. MODES (Apart from a normalization factor) : E1 = x - 2 y (f21 / (f22 -f11 +R)

( 19)

E2 = y + 2 x (f12 / (f22 -f11 +R) With : R = ((f22-f11)2 + 4 f12 f21)1/2, and the fmn are the elements of the Jones matrix M connecting incident Ei and scattered Es fields. . In Eq. 19 the fmn elements are taken in the forward direction, and a factor in order to obtain unitary modulus is subintended for each mode. Let us indicate the relevant factor as F1 and F2 respectively for E1 and E2: F1 = 1 / ( 1 + │f21│2 /│D│2)1/2

(20)

F2 = 1 / ( 1 + │f12│2 /│D│2)1/2 12

With :

D =: = f22 -f11 + R = f22 -f11 + ((f22-f11)2 + 4 f12 f21)1/2.

Eq.18 requires : Re { F1 F2* ( f12* / D* - f21 / D) } = 0

(21)

In Sect.3 the explicit forms for the fmn elements of the Jones matrix M are shown for small ellipsoids , and Eq.21 is verified. Eq.18 can also be verified for cylinders and small chiral spheres 4b. Decomposition of the incident polarization into two orthogonal modes. Given he Jones matrix and the modes of Eq.4, with forms according to Eqs.4,5, the incident vector Ei can be divided into the components Ei1 and Ei2 along the polarization modes E1 and E2 of Eq.1: Ei1 = (Ei • E1*) E1

(22)

Ei2 = (Ei • E2*) E2 One can verify the decomposition correctness. Let us consider a component of Ei1 , say the x component. After decomposition (22) the x component of the sum of Ei1 and Ei2 is : Ex = (Eix A* + Eiy B*) A + (Eix C* + Eiy D*) C =

23)

Eix (A* A + C C*) + Eiy (B* A + D* C)

Where the suffixes x and y indicate the components of the incident field. In Eq. (23) Ex can be seen to be equal to Eix for the case of small chiral spheres, as well as cylinders, rotational ellipsoids. It can be verified that the first parenthesis at second member in (23) becomes equal to 1, while the second parenthesis becomes equal to zero. In an analogous way one can find that the y component of the incident and decomposed fields are equal. 13

In summary, one can verify that orthogonality and correct decomposition are found directly from the mode expressions for cylinders of any orientation, small spherical chiral particles, dielectric rotational ellipsoids, as well as other particles with an axis of cylindrical symmetry . However, for details of an example, let us consider the case of a small dielectric rotational ellipsoid.

14

5. Simple particles. Small dielectric rotational ellipsoids. Some considerations about the scattering by small dielectric ellipsoids, and the cross sections dependence on polarization. As reminded in Sets. 1,2, it can be shown that for a general situation of a non spherical scatterers (ref.1) there are two modes which maintain polarization during propagation. If the electrostatic approximation is adopted, the cross section of a dielectric particle cannot be simply obtained from the optical theorem, since scattered radiation in the forward direction is in phase with the exciting radiation. The cross section , however, can be obtained by considering the electric moment M internal to the particle, and the relevant radiated power W. (analogously to the case of a small dielectric sphere, where the result obtained in this way corresponds to that obtained by the first opportune term in the expansion of the scattered field in spherical functions) Given a reference system x' y' z' for incident radiation propagating in the z ' direction, we now consider the case of a rotational ellipsoid with mayor axis A in the plane x' z', making an angle θ with the x' axis. The two equal (minor) axes are of length B. E' now denotes the incident field, with components E'x, E'y An x y z system is now chosen with the x direction along the major axis of the ellipsoid and the z axis in the plane x' z'. In the electrostatic approximation the expressions for the electric moments components Mx, My, Mz internal to the particle can be found in basic textbooks. They are given as : 24) Mx = V ε0 (εr-1) Ex / D1 Mz = V ε0 (εr-1) Ez / D2

My = V ε0 (εr-1) Ey / D2

D1 and D2 parameters depending on the ratio A/B, and can be derived from to Ref..7, Sect.3.27 by numerical integration. Ex, Ey, Ez are the components of the external static field in the x y z system, and V is the ellipsoid volume. ε0 is the dielectric constant of the external medium. ε = ε0 εr is the dielectric constant of the homogeneous ellipsoid

15

In our case, since E’z = 0 Ex

=

E’x cos θ

Ez

=

- E’z sin θ

Ey

=

E’y

(25)

In (23) : D1 = 1 +(a b2) (ε−ε0) /(2 ε0) A1 D2 = 1 +(a b2) (ε−ε0) /(2 ε0) A2 and : = I [ds( (s+a2)3/2 (s+b2) )]

A1 A2

(26)

= I [ds ( (s+a2)1/2(s+b2)2 ) ]

I [ ] : integral between 0 and ∞ For a / b = 4 and ε / ε0 = 1.5 by numerical integration one obtains : : D1 = 1.0379

D2 = 1.235 ,

The scattered radiated power is proportional to the sum of the square electric moments : If the incident radiation is linearly polarized along the x' direction, one has to consider │M-x │2 + │ M-z│2 , which apart from a constant factor is equal to

0.928 for θ= 0

If the linear polarization is along the y' direction one has to consider │My │2 , which, apart from the same factor , is equal to 0.656 Then the ratio of the cross sections : σx' / σy' for linear polarization of incident radiation along the x' axis and linear polarization perpendicular to this axis) : For θ = 0

:σx'/ σy' =

1.41

For θ = 45° :

16

│M-x │2 + │ M z│2 = 0.793 │My │2

= 0.656

and the ratio σx / σy of the cross sections : σx' / σy' = 1.21 Obviously for θ = 90° σx / σy = 1 for any x’, y’ linear polarizations. In this case of incidence in a direction along the ellipsoid mayor axis any polarization is maintained during propagation The results depend on the ratio a / b and the dielectric constant of the ellipsoid. However, this simple couple of cases show that the extinction cross sections depend on polarization. Thus in general polarization changes during propagation, and the linear extinction coefficient also changes for an otherwise homogeneous medium with suspension of ellipsoids.

17

The following table shows an example of ratio between scattered powers for small dielectric ellipsoids Ratio between the scattered powers relevant to two orthogonal linear polarizations. Dielectric Rotational Ellipsoid . Ratio between major axis A and minor axis B : A/B = 4 Relative dielectric constant εr = 1.5 Propagation of the beam along the z direction Position of the axis A: In the x z plane, making an angle θ with the x axis W1 = power scattered with linear polarization of the incident beam along the x direction W2 = power scattered with linear polarization of the incident beam along the y direction D1 =1.039 (*)

D2 = 1.235 (*)

θ

W1/W2

0

1.413

π/9,

1.364

2π/9

1.242 same results for θ and π − θ

π/3

1.103

4π/9

1.012

π/2

1

(*) D1 and D2 are the parameters according to the formalism for the scattering by dielectric ellipsoids in J.Stratton. Electromagnetic Theory (Ref.7)

18

5A. A case of polarization change. Let us take a simple example. The ellipsoid major axis is in the plane XY (unit vectors x, y) of polarization of the incident field E: E = EX x + EY y. Let EX be in phase with EY : linear polarization The major axis of the ellipsoids makes an angle φ with the X axis. By a rotation by the angle -φ from X to Y one has the components suitable for applying the formalism of the previous section. The unit vectors along the major and minor axes are now indicated as x and y. One has E = Exx + Ey y with linear polarization The electrical internal moments Mx and My are Mx = M Ex/ D1 My = M Ey/ D2 with

M = V ε0(εr-1), V ellipsoid volume.

Ex and Ey are components orthogonal with respect to power propagation. The power extinction can be calculated independently. The extinction cross sections relevant to Ex and Ey are proportionally respectively to Wx and Wy Since radiated power is proportional to the square of the moment, the extinction cross sections are proportional to (1/D1)2 and (1/D2)2 If N is the number of equal ellipsoids (same axes orientation), the medium extinction coefficients σe are proportional to N (1/D1)2 and N(1/D2)2 Then let us write σex = σey =

K (1/D1)2 K (1/D2)2

After a path length Δ in the z direction the power densities propagating are reduced by different factors. Then also the field components are reduced by different factors.

19

However, a question arises : do the phases of the two components change as in the non scattering medium, or the presence of scatterers causes additional phase variation ?. This point will be considered with a NOTE in Appendix 4 and is connected to the possible choice of a proper Monte Carlo.

.5

b The scattered field in the forward direction.

With the geometry of Sect. 3, given the internal moments, the scattered field Es in the forward direction are obtained. Apart from a constant factor including exp (- ikr) / r , t he components along the x’ and y’ axes (reference system for the incident radiation) can be written as :

E’sx = E’x ( cos2θ / (1+A1) + sin2 θ / (1+A2) or : (27) E’sx = E’x (1 / (1+A2) + cos2θ ( 1 / (1+A1) - 1 / (1+A2)) E’sy = E’y (1 / (1+A2)) Since

A2 > A1 , the factor to cos2θ in E’sx is positive.

Thus, let us write : 1 / (1+A2) = 1 ( 1 / (1+A1) - 1 / (1+A2) / (1 / (1+A2))

= α

with α = 0 when the ellipsoids becomes a sphere ( a = b). The parameter α was shown in Sect 3 for the scattering by small ellipsoids, and its form has been given explicitly here.

20

Then, apart from a constant factor, one can write for the fmn elements of the matrix M in the geometry indicated in Sect.3 : (28)

2

f11 = (α-1) cos θ + 1 f22 = 1 f12 = 0 f21 = 0 The forms in Eq. 28 correspond to those quoted from Ref.4, and reported in Sect 3., if one takes into account the correspondence : β = π/2 −θ and φ = 0 in our case More in general : when the dielectric rotationally symmetric ellipsoids is oriented with the major axis a along a direction making an angle λ with the direction z' of propagation, and the projection of the axis on the plane x' y' makes an angle φ with the x' axis, the fmn elements can be obtained by rotating the x’ y’ system by an angle φ , transforming Eqs.15, then by a contrary rotation the system by an angle - φ. One obtains : 2

2

(29)

f11 = (α-1) sin λ cos φ + 1 f22 = (α -1) sin2 λ sin2 φ + 1 f12 = (α -1) sin2 λ cos φ sin φ f21 = (α -1) sin2 λ cos φ sin φ The forms in Eq. 29 correspond to those quoted from Ref.4, and reported in Sect 3 The equations for the modes show that the two conditions for orthogonality and correct decompositions of Sect.4 are fulfilled.

21

Note that the normalization factor F1, F2, mentioned in Eq.20 of Sect.4, is the same for both Modes Also note that the angle λ in Eq. 29 corresponds to angle π/2 − θ , with θ of Fig.1 when φ = 0

6 Possibilities for a Monte Carlo. In a conventional Monte Carlo procedure a simple rule to determine the pathlength L between two subsequent scattering points is based on considering a table of random number R whose value ranges between zero and one, with equal probability for any equal interval of R. To associate an extracted number to the path length L one has to set the relationship between L and R such that the interval ΔR corresponding to ΔL is proportional to the attenuated power at L. Thus one can say that the initial power is distributed” in a correct interval ΔL . This situation is easily obtained in the case of radiation not changing polarization and extinction coefficient σ. For the case dealt in this report one has two MODES with different extinction coefficients of the medium. Thus in the case of two MODES the conventional Monte Carlo procedure needs to be modified. Some approximations are often used. A simple way would be that of considering a weighed average of the extinction coefficient (Ref.8,9). The programmed work for the Florence group considers two possible procedures, which will be used and compared in the next steps of the work. PROCEDURE A. The Stokes vector at an initial point P1 is split in the two vectors corresponding to the two modes. Then one of the two is chosen together with its extinction coefficient by using a random number, and considering the ratio of powers. The trajectory of this mode proceeds till a scattering point decided on the basis of its extinction coefficient following the common procedure based on the extraction of a random number. At the scattering point the phase matrix relevant to the mode allows the definition of the scattering direction, by means of a random number and the dependence of scatterd power 22

on the direction. Given the scattering direction, it is necessary to examine how the scattered field can be divided in the two modes, by considering the new relative directions of field and axes of the scatterers.. The vector is again split and the procedure starts again. One can use the “semianalytic “ procedure to obtain the received power from each scattering point, if the relative geometry of scattering point and receiver, given the receiver area and the field of view. One observes that, for each section of trajectory the discarded mode would really run together with the chosen one. The present procedure neglects the coherent addition of the fields. A second procedure, the one which is presented as Procedure B, is considered to take into account this point. However, the possible different changes of phase for the two modes along the considered trajectory could make the coherent addition much sensible to the value of the run distance between two subsequent scattering points. Procedure B. The distance L from a starting point P1 to the scattering point P2 is chosen on the base of one of the two extinction coefficients, say σ1 for MODE1. Then we take into account that, at an equal interval of L, MODE2 arrives with power equal to that of MODE1 multiplied by the ratio exp (- (σ2 − σ1) L ). Ιn addition, in the same interval of L, the probability of power being scattered for this mode is equal to that of MODE 1 multiplied by σ2 / σ1. From the Stokes Vector at point P1 the E fields of the two component modes are obtained with their relative phase. Then the Stokes Vector for scattering at point P2 is obtained from the E1 field at point P2 summed to the E2 field with the relative phase equal to that of point P1 (propagation with the same phase change), the latter multiplied by the real factor (σ2 / σ1)1/2 exp (- (σ2 − σ1) L/2 ). By taking amplitudes and phases into account the total Stokes Vector ready for scattering is obtained. APPENDICES 1,2,3 are added for an explanation of the procedures. APPENDIX 4 examines a possible effect due to different changes of phase for the two component modes.

23

REFERENCES 1. I.Tsang, J.A. Kong, and R.T.Shin,T. Theory of microwave remote sensing, Wiley, New York, 1988, sect.3.8. 2. A Ishimaru, C.W.Yeh. Matrix representation of the vector radiative-transfer theory for randomly distributed nonspherical particles. J.Opt.Soc. Am. A Vol.1, pp.359-364 (1984). 3. P.Bruscaglioni, Some simple checks for the formulas of sect.3 of a paper by the Florence group, INTAS contract 01-0239, project web-page http://osmf.sscc.ru/~smp/intas.html 4. B.P.Ablitt et al., Imaging and multiple scattering through media containing optically active particles,” Waves in Random Media 9, pp. 561-57, 1999. 5. M.I.Mishchenko, J.W.Hovenier, and L.D.Travis, Lightt Scattering by non-spherical Particles, Academic Press, San Diego, 2000. 6. C.F.Bohren, D.R.Huffman. Aborption and scattering of light by small particles. Wiley. N.Y.1983. 7. J.A.Stratton. Electromagnetic Theory.Mc Graw. N.Y.,1941 8. S.Wolf, N.V.Voshchinnikov, and Th-Henning, “Multiple scattering of polarized radiation by non-spherical grains. First results. ,” Astron. Astrophys. 385, pp. 365-376, 2002, eq. 7.. 9. P.G.Martin, “interstellar polarization from a medium with changing grain alignments,” Ap.J.. 187, pp. 461-472, 1974.

24

APPENDIX 1 Propagation and scattering for the two "modes" (The symbol I [f(x)dx] indicates the integral in dx between 0 and ∞) The Monte Carlo codes used for multiple scattering have been based, as usual, on evaluating the probability of scattering after a trajectory length. Let us start with one "mode" : mode 1. We consider a random number R, whose probability density p( R) is uniform between the values 0 and 1, and with integral = 1 over the extension 0,1. If σ1 is the linear extinction coefficient for power of mode 1. power scattered at the distance L in the interval dL is (apart from a constant factor) E012exp(-σ1L) dL, with E01 field at L = 0. The relevant probability must be proportional to σ1exp(-σ1L) dL. This leads to the relationship (A1) L = - (1/σ1) ln(R). that leads to R = exp(-σ1L), and the relevant interval ⏐d( R)⏐ is proportional to σ1exp(-σ1L) dL Given the mode 2, with E02 field at L = 0 and extinction coefficient σ2, the probability of scattering of power at the same L and interval dL is equal to that of mode 1 time the factor E022 / E012(σ2/σ1)exp(-(σ2−σ1)L) dL Thus, if, as usual, after extracting the value R (equiprobable between 0 and 1) and putting power of mode 1 equal to 1, power of mode 2 scattered al L in the interval dL has to be putted as (A2) W2(L) = (σ2/σ1)exp(-(σ2−σ1)L) E022 / E012 The total power scattered for field 1 is I[W2(L) p(L) dL] One has p(L) dL = σ1exp(-σ1L)dL Then: 25

I [W2(L) p(L) dL] = E022 / E012 To define the Jones matrix at L for a certain scattering direction , one has to consider the field components E1 and E2 at L active for scattering. That is for the mode 1 and mode 2 one has to assume the fields (A3) E1(L) = 1

and

E2 (L) = (E02 / E01) (σ1/σ2)1/2exp(-(σ2−σ1)L/2) However, for the case of ellipsoids, the question of the phase of the scattered field with a possible difference between the two modes will be considered in the Appendix 4..

26

APPENDIX 2 The case of propagation with two modes with different extinction coefficients. One can find the probability p1 that one mode (say mode 1) is scattered at the distance X in the interval dX. This is found by the usual rule A1) X = -(1/σ1) ln R

And the searched for probability is :

A2) p1 = - σ1 dX (1/σ1) ln R Given p1, the probability p2 of scattering in the same interval dX at the same distance X for the mode 2 is given by the conditioned probability : p2 = p1 * (p2/p1) That is A3) p2 = p1 (σ2/ σ1) exp(-(σ2- σ1)X) Thus one can extract the distance X for one of the modes by Eq.A1 Power scattered for this mode is given as A4) │E1│2 with the correct probability. The power scattered at X, dX, for mode 2 Is then given as A5) │E2│2 (σ2/ σ1) exp(-(σ2- σ1)X) If one prefer to deal with fields the scattered field of the two independent modes are obtained bi roots of powers given by A4 and A5,

27

APPENDIX 3. A line for a Monte Carlo code. On the basis of the preceding considerations, we suggest the following line to be followed for a Monte Carlo, with the following steps START propagation direction : z axis DECOMPOSITION of the incident field in the two modes according to the properties of suspended particles and the relevant elements of the j matrix. EVALUATION of the linear extinction coefficients σ1 and σ2 for the two modes DEFINITION of the length of the trajectory to the first scattering point taking into account the different coefficients DIFFERENT POWER extracted from the two modes, and form of the total field scattered in the different directions. EVALUATION of the scattered power in the different directions, with evaluation of the different probability of scattering according to the directions. CHOICE of the scattering direction from the extraction of random numbers. DEFINITION of the scattered field (after normalization) in the chosen direction. New direction of propagation THE SCATTERERS show a certain attitude with respect to the incident propagating field. NEW DEFINITION of the two modes with respect to the propagation direction , particles attitude and elements of the relevant J matrix (fmn elements in the propagation direction) NEW DEFINITION of linear extinction coefficients

:

28

APPENDIX 4. A NOTE on the phase propagation. One has to consider the phase relation between the two modes , and to what extent it is maintained during propagation or not. The question is how to connect an attenuation factor to the variation of the phase difference between the two modes at the distance where the attenuation is considered. For a simple case, here the case of small dielectric ellipsoids is considered, with the formalism of Sect.5 , where the relationship between incident field Ei and scattered field Es is given by the matricial formalism N1) Es = (exp(ikr)/R M Ei, with M Jones matrix and the the electrostatic approximation is used. One can take into account that for the case of rotational ellipsoids the matrix can become diagonal by a choice of the reference system where the x axis is along the projection of the major axis of the ellipsoid on the x y perpendicular to the propagation direction z. Then the two cross sections and the two extinction coefficients for the mode are obtained by N2) Ce = (4π/k) Im( f •u*) with f in the forward direction. For the small ellipsoid one has linear polarization modes along the x axis (mode x), and y axis (mode y),. Then in Eq.2 f is equal to f11, and f22 for linear x polarization (u = x) and y (u = y) polarization respectively. In the forward direction the two nondiagonal elements f12, f21 are zero. Then N3) Cex = (4π/k) Im(f11)

and Cey = (4π/k) Im(f22)

In the first approximation, when the electrostatic approximation is used, since the internal electric moments within the particles is taken in phase with the incident field, the scattered electric field in the forward direction for one particle is assumed to be in phase with the undisturbed field. Thus the optical theorem cannot be used for calculating cross sections and attenuation from the elements f11., or f22, assumed to be real.

29

However a correction to these diagonal elements can be obtained as follows in order to introduce an imaginary part of the diagonal f elements. One can remember the case of a small dielectric sphere. For this particle, analogously to the case of the small dielectric ellipsoid the total scattered power and the cross section are obtained by taking into account the internal electric moment and its total emitted power. On the other hand the optical theorem gives the same cross section if the first imaginary term in the expansion of scattered field in spherical wave is used for f. (Ref.2). We extend the reasoning to the case of the ellipsoid. For each mode the total scattered power is obtained from the internal moment. Then, from the field in the forward direction according to this moment, one obtains the imaginary part of f. Mode x. We indicate as W the power scattered by the ellipsoid : N4) W =1/(2Z) Ce │Ei│2, which by the optical theorem is : N5) W = 1/(2Z) (4π/k) Im ( f11) │Ei│2 From the internal electric moment M we also have : N6) W = (1/3) π Z ω 2/λ2 │M│2 With : N7) M = ε0(εr−1) V Ei / Dx where Dx is the parameter indicated in Ref(X) Sect.XX.(see Sect. XX of thispaper) Then : N8)

W = 1/3 π Z ω 2/λ2 │ε0(εr−1) V Ei / Dx│2

From the real part of the internal moment one has for the component x: N9) Re (f11) = (Z ω κ /4π) Μ / ΕI = (Z ω κ /4π) V ε0(εr−1) / Dx and for the imaginary part: N10) Im (f11) = (κ Z2 ω2 / 6 λ2) │ M / Ei│2

30

Then we have the ratio : N11) Im (f11) / Re (f11) = (2 π ω Z / 3 λ2 ) │ M / Ei│ = = (2 π ω Z / 3 λ2 ) V ε0(εr−1) / Dx = = (16 / 9) π3 ( a / λ)3(εr−1) / Dx = K / Dx Analogously we have for the ratio N12) Im (f22) / Re (f22) = (16 / 9) π3 ( a / λ)3χDy

= K / Dy

where a is the radius of the equal volume sphere. With the possible values of the relevant parameters: εr−1 = .5,

a / λ = 0.1,

Im (f11) / Re (f11) 0.022

Dx = 1.05,

= 0.027

Dy = 1.25 one has: Im (f22) / Re (f11)

=

This second approximation shows that the scattered field is in forward direction is a complex quantity, with an imaginary part small in comparison with the real part. However the scattered field from each particle, or from a small volume element, has the property of a spherical wave travelling independently from the incident field. As shown in Ref.1 and Ref.2, by considering the coherent field propagating in the forward direction one has d Ex = exp ( i (k + Re( f11) + Im (f11)) dz d Ey = exp ( i (k + Re (f22) + Im (f22)) dz

(k wavenumber)

That is : during propagation in the forward direction generally both amplitude and phase change in a different way. For the two modes By referring to the simple example considered in Sect. XXX, let us start at a point (0) where the difference of phase of the mode y with respect to the mode x is Φ. AT a distance L let the attenuation of Mode x is A : Ex(L) / Ex (0) = exp(-A) exp ( i (kL + Φx)) One has : 31

Ey(L) / Ey (0) = exp(-A Im(f22) / Im(f11)) exp ( i (kL + Φy)) with Φy - Φx = Φ + A (Re(f22) – Re(f11)) / Im(f11)

32

APPENDIX 5 Small dielectric ellipsoids.: the scattered field. Geometry : the z' axis (z') is the direction of propagation of the incident .field. Ex' and Ey' are the components of the incident field. the major axis of the rotational ellipsoid is in the plane x' z' and makes an angle θ with the x' axis this situation can be obtained by rotation of previous defined axes and is suitable for calculation of the scattered field, since ex' nd ey' are according to the two "modes". the direction of the scattered field (s) makes an angle λwith the z' axis and the projection on the x' y' plane mkes an angle φ with the x' axis. for calculating the scattered field the reference axes x, y z (x,y,z) are defined. x is along the ellipsoid major axis. thus it makes an angle θ with the x' axis. A5.1) x = x' cosθ + z' sinθ z = -x' sinθ + z' cosθ y = y' the components of the incident field are then : A5.2) Ε'x cosθ

along

- Ex' sinθ

along z

Ey'

along y' = y

x

the internal electric moment due to ex' has the components :

33

A5.3) Mx = (1/D1) Ε'x cosθ x

= (1/D1) Ε'x (x' cosθ + z' sinθ) cosθ

Mz

= - (1/D2) e'x sin z -= - (1/D2) e'x (-x' sinθ + z' sinθ) sinθ

My

=

(1/D2) ey' y

where d1 and d2 are given in sect.1. apart from the common subintended factor v ε0 (εr-1).. v volume of the ellipsoid. one then calculate the scattered field in the direction s. in the following expressions the further factor (ωkz / 4π) exp(ikr)/r is subintended. the scattered field in the diirection s is the sum of the following contributions

field due to mx. (apart from the factor e’x) : A5.4) (1/D1) cosθ ( (- cos2λcos θ + cos λ sinλcosφ sin θ -sin2λ sin2φ cos θ) x' + + (cos λ sinλ sinφ sin θ + sin2 λ sinφ cosφ cos θ) y' + +( sin λ cosφ cos λ cos θ - sin2 λ cos2φ sin θ -sin2 λ sin2φ sin θ) z' ) field due to mz (apart from the factor e’x): A5.5) (1/D2) sin θ ( (-cos2λ sin θ - sin λ cosΛ cosφ cosθ-sin2 λ sin2φ sin θ)x' + +(-cos λ sin λ sin φ cos θ + sin2 λ sinφ cosφ sin θ) y' + + (sin λ cos Λ cosφ sin θ + sin2 λ sin2φ cos θ + sin2 λ cos2φ cos θ)z' ) field due to my (apart from the factor e’y): 34

A5.6) (1/D2) (sin2λsinφ cosφ x' -(cos2 λ+ sin2 λcos2φ)y' + sinλcosλsinφ z')

simple cheks one can see the following respected properties of the scattered field : the s direction is perpendicular to the scattered field. this is verified by considering the parts with coefficient e’x , and separately the parts proportional to 1/d1 and 1/d2. equally for what concerns the part proportional to e’y. when the ellipsoid axes are equal (sphere), d1 = d2 and the angle θ does not appear in the sum of the field components. when λ = 0 : s in the direction z' of incidence : there is no z' component of the scattered field. in this direction mx and mz (due to the x' component of the incident field) produce no component of the scattered field along y' , and my only produces a component along the y' direction. it is then verified that the modes correspond to linear polarization along the projection on the x' y' plane of the ellipsoid mayor axis, and linear polarization perpendicular to this direction. on he other hand, one can see that for scattering in a direction different from z’ ( the direction of propagation of the incident field) the scattered field due to e’x has a component trasverse to the x’z’ plane, and the scattered field due to e’y has generally a component in the x’z’ plane. thus the scattered field does not generally belong to one of the modes. when : θ = 0, and s along the y' direction ( λ = π/2, φ = π/2). there is only a component of the scattered field in the -x' διρεχτιον when θ = 0, and s along the x' direction ( λ = π/2, φ = ο).. there is only a component of the scattered field in the -y' direction due to my. .. Other simple checks are easily verified.

35

36

monte carlo procedures. a problem for multiple ...

30 Jun 2007 - Given a reference system x' y' z' for incident radiation propagating in the z ' direction, we now consider the case of a rotational ellipsoid with mayor axis. A in the plane x' z', making an angle θ with the x' axis. The two equal. (minor) axes are of length B. E' now denotes the incident field, with components E'x ...

293KB Sizes 0 Downloads 235 Views

Recommend Documents

Sequential Monte Carlo multiple testing
Oct 13, 2011 - An example of such a local analysis is the study of how the relation ... and then perform a statistical test of a null hypothesis H0 versus. ∗To whom ... resampling risk (Gandy, 2009), and prediction of P-values using. Random ...

Sequential Monte Carlo multiple testing
Oct 13, 2011 - can be reproduced through a Galaxy Pages document at: ... Then, in Section 3, we show on both simulated and real data that this method can ...

a monte carlo study
Mar 22, 2005 - We confirm this result using simulated data for a wide range of specifications by ...... Federal Reserve Bank of Kansas City and University of Missouri. ... Clements M.P., Krolzig H.$M. (1998), lA Comparison of the Forecast ...

Hamiltonian Monte Carlo for Hierarchical Models
Dec 3, 2013 - eigenvalues, which encode the direction and magnitudes of the local deviation from isotropy. data, latent mean µ set to zero, and a log-normal ...

A Non-Resampling Sequential Monte Carlo Detector for ... - IEEE Xplore
SMC methods are traditionally built on the techniques of sequential importance sampling (SIS) and resampling. In this paper, we apply the SMC methodology.

Monte Carlo Simulation
You are going to use simulation elsewhere in the .... If we use Monte Carlo simulation to price a. European ...... Do not put all of your “business logic” in your GUI.

A quasi-Monte Carlo method for computing areas ... - Semantic Scholar
Our method operates directly on the point cloud without any surface ... by counting the number of intersection points between the point cloud and a set of.

A Sequential Monte Carlo Method for Bayesian ...
Sep 15, 2002 - to Bayesian logistic regression and yields a 98% reduction in data .... posterior, f(θ|x), and appeal to the law of large numbers to estimate.

Monte Carlo simulations for a model of amphiphiles ...
Available online at www.sciencedirect.com. Physica A 319 (2003) .... many sites. The amphiphiles possess internal degrees of freedom with n different states.

Statistical Modeling for Monte Carlo Simulation using Hspice - CiteSeerX
To enable Monte Carlo methods, a statistical model is needed. This is a model ..... However, it is difficult to determine the correlation without a lot of statistical data. The best case .... [3] HSPICE Simulation and Analysis User Guide. March 2005.

Using the Direct Simulation Monte Carlo Approach for ...
The viability of using the Direct Simulation Monte Carlo (DSMC) approach to study the blast-impact ... by computing load definition for two model geometries - a box and an 'I' shaped beam. ... On the other hand, particle methods do not make the conti

accelerated monte carlo for kullback-leibler divergence ...
When using ˜Dts a (x) with antithetical variates, errors in the odd-order terms cancel, significantly improving efficiency. 9. VARIATIONAL IMPORTANCE ...

Introduction to Monte Carlo Simulation
Crystal Ball Global Business Unit ... Simulation is the application of models to predict future outcomes ... As an experimenter increases the number of cases to.

Fundamentals of the Monte Carlo method for neutral ...
Sep 17, 2001 - Fax: 734 763 4540 email: [email protected] cс 1998—2001 Alex F .... 11.3 Convergence of Monte Carlo solutions . . . . . . . . . . . . . . . . . . . . . .

Monte Carlo Null Models for Genomic Data - Project Euclid
To the rescue comes Monte Carlo testing, which may ap- pear deceptively ... order the null models by increasing preservation of the data characteristics, and we ...

Introduction to Monte Carlo Simulation - PDFKUL.COM
Monte Carlo Simulation Steps. • Following are the important steps for Monte Carlo simulation: 1. Deterministic model generation. 2. Input distribution identification. 3. Random number generation. 4. Analysis and decision making ..... perform output

Copy of CMG14 monte-carlo methodology for network capacity ...
Quality of Service (QoS) = a generic term describing the performance of a system ... We have: 1. A network a. Topology, including Possible Points of Failure b.

Particle-fixed Monte Carlo model for optical coherence ...
Mar 21, 2005 - Y. Pan, R. Birngruber, J. Rosperich, and R. Engelhardt, ..... complex partitioning schemes may be possible (like oct-tree or kd-tree [17]), they help little ..... its cluster and Dr. Tony Travis for the tip of programming on the cluste

Monte Carlo Simulation for the Structure of Polyolefins ...
unsaturated (usually vinyl) groups can be incorporated by the. CGC into a ... broaden the molecular weight distribution, since chains formed at the second site will .... to published experimental data for a lab-scale synthesis of a dual catalyst ...

Monte Carlo null models for genomic data - Semantic Scholar
Section 3 presents null model preservation hierarchies and signifi- cance orderings. Sections 4–6 ... A Galaxy Pages (Goecks et al., 2010) document allowing for ...

Dead Reckoning for Monte Carlo Localization in Low ...
Especially in scenar- ios with low seed densities, we find SA-MCL to significantly outperform MCL. S. S. S. S. S. Figure 6: Setup of demo with seeds and RC cars.

(Quasi-)Monte Carlo Methods for Image Synthesis
Previously, he was an R&D engineer at Industrial. Light and Magic where he worked on a variety of rendering problems in production. His current research interests include interactive global illumination and rendering algorithms,. Monte Carlo methods,