20 second Mie scattering attempt

Second Attempt at Mie Scattering Coefficients PhL 2.15.03 updated 2.16.03 Plane Wave Expansion. See details in Appendix A of previous document, but this is the basic idea. Write down Jackson's exact plane wave results 16.139. E = (i) [ j(kr) X,1 (1/k) x { j(kr) X,1 } ] B = (i) [ ∓ i j(kr) X,1 - (i/k) x { j(kr) X,1 } ] Then take the appropriate linear combination of his results that corresponds to E = (E+ + E- )/2 = eikz B = (B+ + B- )/2 = eikz // see 16.130 and fact that B = x E His plane wave formulas then read, E = (1/2)(i) [ j(kr) [ X,+1 + X,-1] + (1/k) x { j(kr) [ X,+1 - X,-1] } ] B = (1/2) (i) [ - i j(kr) [ X,+1 - X,-1] - (i/k) x { j(kr) [ X,+1 + X,-1] } ] Then make the following definitions for o = odd, and e = even: Mo1 j(kr) r x o1 where o1 P1(cos) sin Me1 j(kr) r x e1 where e1 P1(cos) cos (note that sin is "odd" in ) with corresponding formulas for the N functions No1 (1/k) x Mo1 = (1/k) x [j(kr) r x o1 ] Ne1 (1/k) x Me1 = (1/k) x [j(kr) r x e1 ] Notice that the "pilot functions" are not Ym's. They are in fact the sum and difference of two Ym's of opposite m values, allowing us to make this replacement in Jackson's formula, [ Y,+1(,) + Y,-1(,)] = 2 i o1 [ Y,+1(,) - Y,-1(,)] = 2 e1 Then Jackson's plane wave formulas shown above boil down to this: Ep = e ikz = i [ M(1)o1 - i N(1)e1 ] Bp = e ikz = i [ - M(1)e1 - i N(1)o1 ] In the B equation, the first term is a (- i) from the Jackson formula, but then we divide by i because the difference of Y's does not have the i that the sum did have, and we change o to e as well. For the second term, we get -i times the N function that matches the first term in the E equation. In general, the rule as stated below is that you cross multiply coeffs by -i and change the nature odd/even. This result agrees with (1.6) in all respects except one: in the derivation above, we clearly see that the m-index should have the value 1, as shown in the subscript labels above. But in (1.6), this subscript appears as the general m (which they happen to call ). This is an error in (1.6). In the Krugel PDF paper, the above formula appears as equation 2.44, exactly as I have stated it. The superscript (1) implies that we are using the j radial functions. Recall that Jackson's "first term" in the E equation 16.139 is the "magnetic" or TE term, and this is still true when we combine the two helicity solutions. Thus, we can associate "magnetic" with "odd", and of course "electric" with "even". That is a nice label to have on the M and N functions. So if the first term is odd, we always expect the second term to be even. Rule for getting B from E: Suppose you only had the Ep equation above, how could you at once deduce the Bp equation? Look at the general 16.47 expansion. If we choose to superpose a few different values of the Xm, the form should stay the same. So if we have the E equation, we get the B by these rules: (1) including the (1/k) factor as part of the derived term structure, in the B expansion the coefficient of the derived term is -i times the coefficient of the E main term. (2) the main term coefficient in the B expansion is -i times the derived term coefficient in the E expansion. (3) in the B expansion, the main term matches the "nature" of the E expansion derived term, and vice versa. This nature might be even/odd, or bessel type or what have you. So let's start with: Ep = e ikz = i [ 1 M(1)o1 - i N(1)e1 ] and identify M as the main term and N as the derived term. Then the three rules would imply this for B: Bp = i [ - M(1)e1 - i N(1)o1 ] and this agrees with the form above obtained by other means. Summary of the E to B rules: (1) B derived coeff = -i * E main coeff (2) B main coeff = -i * E derived coeff (3) B main nature = E derived nature Summary: To get B from E, for each term use the nature of, and -i times the coeff of, the cross term. Aside on Krugel: I stumbled onto the contents and first three chapters on-line of the following book: E. Krugel, The Physics of Interstellar Dust, 12/02 $135 Amazon He is the one who cleared up the above plane wave mystery. Krugel is at the Max Planck Institute (MPI) for Radio Astronomy in Bonn. I notice that the PDF documents cannot be printed, nor can anything be selected for cut and paste! Attempt to set up the scattering problem in the N and M world. First, let's look back at our earlier derivation where I wrote down the four boundary conditions. My two curl conditions implied matching of tangential field components. The r x B condition is problematical for a conductor because there really is a surface current. However, If you use H instead, H does not see the possibly magnetization surface current, the way D does not see polariozation surface charge. So I guess I accept Krugel's two r x boundary conditions in the general case of arbitrary isotropic and , which is what he is talking about. Krugel then notes that each of these r x boundary conditions implie 2 equations if you use the and directions for tangentialness. He says you then get four equations for the four unknown coefficients and you can then solve. When I did this problem the first time (see separate paper) , I was indeed finding redundancy with my boundary conditions. So I guess the answer is that the two r x conditions are all you need! Here is Krugel's setup of the problem: Ep = K [ M(1)o1 - i N(1)e1 ] K = i // incident plane wave Es = K [ - b M(3)o1 + i a N(3)e1 ] // scattered field, (3) means h(1) (kr) Ei = K [ c M(1)o1 - i d N(1)e1 ] // internal field, (1) means j(kr) Notice the slightly unusual way the a and b coefficients are defined by these equations, probably this is to make the resulting formulas as simple as possible. I hope I can get that far! *** I have concluded that Krugel is missing a factor of i in front of his a . If you don't do this, a comes out imaginary. The factor of i is present in my little 4-page Mie theory paper, which makes me think this is correct. Now we write the H equations using our above rule for going from E to B noted above: " Summary: To get B from E, for each term use the nature of, and -i times the coeff of, the cross term. " H'p = (1/o) K [ - M(1)e1 - i N(1)o1 ] K = i H's = (1/i) K [ a M(3)e1 + i b N(3)o1 ] H'i = (1/i) K [ - d M(1)e1 - i c N(1)o1 ] Notice that we have to put in factors of due to the conversion from B to H. BUT, we are not quite done! As described in my notes for Jackson chapter 16, Appendix B, the fields shown above are really the H' fields such that H' = H where = 1/ = v/c. It is the primed B (or H) fields which appear in Jackson's multipole formulas, so we have to account for that here with one more restatement of the H fields. Define 1/ = as n, the (possibly complex) index of refraction. This appears as m in the Krugel work, and is called the "optical constant". So we restate one more time as Hp = go K [ - M(1)e1 - i N(1)o1 ] K = i Hs = go K [ a M(3)e1 + i b N(3)o1 ] Hi = gi K [ - d M(1)e1 - i c N(1)o1 ] / i = inside where we are now careful to label the index for i = inside the sphere, and o = outside the sphere. Here we have introduced the quantity g = n/ = / = Now the boundary conditions are of the general form Xi = Xp + Xs evaluated at r=a, because both the plane wave and the scattered wave are on the "outside", while the internal is on the "inside". So we can write our boundary conditions as r x Ei = r x Ep + r x Es r x Hi = r x Hp + r x Hs In each expansion, we can then remove the sum and ignore the constant X overall factor and we end up with inside plane scattered r x [ c M(1)o1 - i d N(1)e1 ] = r x [ M(1)o1 - i N(1)e1 ] + r x [ - b M(3)o1 + i a N(3)e1 ] r x [- gi d M(1)e1 - i gi c N(1)o1 ] = r x [ - go M(1)e1 - go i N(1)o1 ] + r x [ go a M(3)e1 + i go b N(3)o1 ] So obviously we now have to go compute the four r x products shown in these formulas. Let's do one example to see how this goes with an M function: Mo1 j(kr) r x o1 No1 (1/k) x Mo1 = (1/k) x [j(kr) r x o1 ] r x M(1)o1 = r x j(kr) r x o1 = j(kr) r x r x o1 = - j(kr) r2 o1 (appendix A1) r x N(1)o1 = r x (1/k) x [j(kr) r x o1 ] = (1/k) r x x [j(kr) r x o1 ] = (1/k) r x { g r 2 - (g + r g' ) } where g = j(kr) = - (1/k) (g + r g' ) r x = - (1/k)(rg)' r x where the k in (1/k) is specific to the medium so write it as (1/k)(1/n) where k = /c. [ My notation is weak in this respect.] Now some more definitions, o = o1 e = e1 o = odd e = even jo = j(kor) ji = j(kir) o = outside i = inside ho = h(1) (kor) then r x M(1)o1 = f * - r2 o r x N(1)o1 = (rf)' / n * (1/k) r x o // k here on right is /c but you have to put in the correct Bessel function, and the correct subscript on it! We shall now "process" each of the two BC equations. The first is this, recalling Xi = Xp + Xs, inside plane wave scattered r x [ c M(1)o1 - i d N(1)e1 ] = r x [ M(1)o1 - i N(1)e1 ] + r x [ - b M(3)o1 + i a N(3)e1 ] The M terms are all going to have the common factor - r2 o . We look at their coefficients: c ji = 1 jo - b ho The N terms are going to have the common factor - (1/k) r x e , and the coefficient match gives - i d (rji)'/ni = - i (rjo)'/no + i a (rho)'/no The second equation (magnetic match) for processing is this, inside plane wave scattered r x [- gi d M(1)e1 - i gi c N(1)o1 ] = r x [ - go M(1)e1 - go i N(1)o1 ] + r x [ go a M(3)e1 + i go b N(3)o1 ] The M terms are all going to have the common factor - r2 e . We look at their coefficients: - gi d ji = - go jo + go a ho The N terms are going to have the common factor - (1/k) r x o , and the coefficient match gives - i gi c (rji)'/ ni = -i go (rjo)'/ no + i go b (rho)'/ no So let's gather up these four equations into one place: c ji = 1 jo - b ho - i d (rji)'/ni = - i (rjo)'/no + i a (rho)'/no - gi d ji = - go jo + go a ho - i gi c (rji)'/ ni = -i go (rjo)'/ no + i go b (rho)'/ no But now let's give them one round of simplification before going further, using n ni/no and g gi/go // = n (o/i) by the way c ji = jo - b ho - d (rji)' = - n (rjo)' + a n (rho)' - g d ji = - jo + a ho - c g (rji)' = -n (rjo)' + b n (rho)' Now before continuing, we are going to make some more definitions. We are going to re-use the symbol as follows: (xo) xo j(xo) xo = ko r " Riccati-Bessel functions" The above equations contain things like jo = j(xo) where xo = ko a in the case that we are "outside". But ko = k/o = no k where k is the speed of light value /c xo = no ka = no x So let's make the replacement jo = j(xo) = (xo)/xo = (xo)/(nox) = o / (nox) Now another occurrence is this: (r jo ) ' = r [ r j (kor) ] = x [x j (x) ] = x [(x) ] = '(xo) Similarly, define the function to go with the h function. Now let's process the equations using these changes: jo = o / (nox) (r jo ) ' = 'o ho = o / (nox) (r ho ) ' = 'o ji = i / (nix) (r ji ) ' = 'i BEFORE: c ji = jo - b ho - d (rji)' = - n (rjo)' + a n (rho)' - g d ji = - jo + a ho - c g (rji)' = -n (rjo)' + b n (rho)' AFTER: c i / (nix) = o / (nox) - b o / (nox) - d 'i = - n 'o + a n 'o - g d i / (nix) = - o / (nox) + a o / (nox) - g c 'i = -n 'o + b n 'o And simplify again to get, c i = n o - n b o - d 'i = - n 'o + a n 'o - d g i = - n o + a n o - c g 'i = -n 'o + b n 'o Now finally rewrite these in more or less standard form: i c + n o b = n o 'i d + n 'o a = n 'o g i d + n o a = n o g 'i c + n 'o b = n 'o These are trivial to solve since we have two sets of equations each with only 2 variables, but let's have Maple do the solution, just for Maple practice. Change to these symbols sI*c + n*zO*b = n*sO s = "sigh" = sIp*d + n*zOp*a = n*sOp z = "zie" = g*sI*d + n*zO*a = n*sO I = inside, O = outside g*sIp*c + n*zOp*b = n*sOp suffix p = prime I have fed these into Maple and here is what I get: n (-zOp sO + sOp zO) g sI sOp - sO sIp -g sIp sO + sOp sI n (-zOp sO + sOp zO) {d = - --------------------, a = -----------------, b = ------------------, c = - --------------------} g sI zOp - sIp zO g sI zOp - sIp zO -zO sIp g + zOp sI -zO sIp g + zOp sI Notice that the two scattering coefficients a and b are functions only of g. When both media have the same , then we can identify g with the optical constant n (ie, the complex index). But in the general case, g is not the same as the "optical constant". I will now translate this into the fancy notation. For comparison, let's make these replacements: g n // this says that both media in and out have the same n m sI (mx) sIp '(mx) sO (x) sOp '(x) zO (x) zOp '(x) Now translate the above results: anum = m (mx) '(x) - (x) '(mx) aden = m (mx) '(x) - (x) '(mx) bnum = (mx) '(x) - m(x) '(mx) bden = (mx) '(x) - m(x) '(mx) cnum = m [(x) '(x) - (x) '(x)] cden = (mx) '(x) - m(x) '(mx) // same as bden dnum = m [(x) '(x) - (x) '(x)] // same as cnum dden = m(mx) '(x) - (x) '(mx) // same as aden So we can now construct our results a = // agrees with Krugel a (he has num and den negated) // agrees with (3.9) for a in 3-page clip b = // agrees with Krugel b (he has num and den negated) // agrees with (3.10) for b in 3-page clip c = // see below for Wronskian version d = // see below for Wronskian version Now, the c and d have Wronskians in the numerator. Amazingly, page 445 of AS has some data on this: W[ z jn(z), z yn(z) ] = 1 // note that y = n, two notations for the second kind functions! But we know that h1 = j + i n, therefore z h1 = z j + i z n W[ z jn(z), z h1n(z) ] = W[ z jn(z), i z n] = i W[ z jn(z), z yn(z) ] = i So we know now that W[(x), (x)] = i = (x) '(x) - '(x) (x) = i So we can now simplify the results for c and d, c = = (-m) * d(3.12) d = = (-m) * c(3.11) A comment in this 3-page clip say that BH's results are swapped (the way mine are) and are larger by factor m (the way mine are). But I also have an extra minus sign in both my results. Of course people really care about the a and b results since these are the scattering results. The c and d are only for the internal fields that usually no one cares about. The results quoted in the 3-page clip are those of van der Hulst. Probably this guy was just trying to simplify things. Now, can I "rescue" the results of my 4-page paper somehow? I guess when data is presented that way, I have to interpret the prime on [ ]' as differentiation with respect to x. In that case, we note that, /x [ mx j(mx) ] = m /(mx)[ same ] = m (mx). In this case we can now compare: a(4pp) = -b(me) // swap and - signs 4pp = "four page paper" b(4pp) = -a(me) c(4pp) = (1/m) c(me) // factor 1/m d(4pp) = (1/m) d(me) So now at least the results are reasonable, no quadratic powers of m anywhere in reality. However, the results for and and b are therefore WRONG for the setup the author has shown, which is the same setup as with me and with Krugel. But usually the sum of squares is what matters, so OK... Repair for general According to my Jackson modification, I am supposed to do H = H'/ = n H'. This is exactly what I did. So the results are OK for general already, I don't know what I was thinking of. Oh yes, the fact that H = B/ does not affect this scaling rule! Restatement of results for two different media. Make these reverse subs into our comparison results, i = (mx) 'i = '(mx) m = ni / no o = (x) 'o = '(x) o = (x) 'o = '(x) Also, let's keep our results fully general in terms of m and g = gi / go where gi = ni/i. anum = g i 'o - o 'i aden = g i 'o - o 'i bnum = i 'o - g o 'i bden = i 'o - g o 'i cnum = n [o 'o - o 'o] cden = i 'o - g o 'i // same as bden dnum = n [o 'o - o 'o] // same as cnum dden = gi 'o - o 'i // same as aden a = b = c = n = d = n = where n = ni /no = index inside / index outside i = sqrt(-1) g = gi / go = ni /no * o/ i Note that g = n when o = i (x) = x j(x) (x) = x h1(x) the Riccati-Bessel AS page 445 ( just regular spherical Bessels mult by x) o = outside, xo = noka i = inside, xi = nika a = sphere radius, k = /c Examples: o = (xo) = xo j(xo) = noka j( noka) 'o = '(xo) = { x (x) }|x = xo = x [x j(x)] |x = xo = { x j'(x) + j(x) } |x = xo = { xo j'(xo) + j(xo) } = noka j'( noka) + j(noka) Test Limits: (1) If m = 1, we find that a = b = 0, so there is no scattering. Also, we get that c = d = 1. This is very reasonable looking back at our formulas: Ep = E [ M(1)o1 - i N(1)e1 ] E = i // incident plane wave Es = E [ - b M(3)o1 + i a N(3)e1 ] // scattered field, (3) means h(1) (kr) Ei = E [ c M(1)o1 - i d N(1)e1 ] // internal field, (1) means j(kr) In this limit, Ei is identical to Ep, the "internal field" is just the plane wave! (2) Suppose we have a conductor for our sphere. In this case, according to Krugel 1.116 for "low" frequencies we can say ni = For a perfect conductor, this goes to infinity along the angle ei/4. So let's take m and see what happens. We need some facts first for the c and d limits only, i = (x) = x j(x) K x+1 'i = '(x) (+1) K x see jackson page 540 These both blow up since xi = nika, and this means that c = d = 0. Meanwhile the scattering coefficients are finite, determined by a and the outside index no : a = = = (1/2) [ 1 + ] b = = = (1/2) [ 1 + ] and amazingly enough, we have exactly duplicated Jackson's page 570 result for the conducting sphere. Notice that he has (1/2) out in front of 16.141. We do have to think about how his X is scaled compared to our M and N stuff, but I think it comes out right. So these were two good limits to take I think! Let's pause for a rest now, we finally have the Mie Scattering coefficients figured out right! Appendix A. 1. r x r x (,) = (r) r - (rr) = - r2 2. x [g(r) r x (,) ] = g x r x + g x r x = 1 + g*2 1 = g x r x = (g ) r - (g r) = 0 - g' r = - g' r 2 = x r x = r 2 - (1 + r r ) = r 2 - x [g(r) r x (,) ] = - g' r + g r 2 - g = g r 2 - (g + r g' ) Appendix B. Regarding the rule for finding B from E Assume we have an expansion for E, such as this Es = E [ - b M(3)o1 + i a N(3)e1 ] // scattered field, (3) means h(1) (kr) where we are in the outer medium having o and o . We use (7.1) as our Maxwell's equations in a medium like this, and we see that: B = -i/kf x E where k = /c. f = free space Now we know from our work on M and N functions that: x M = kN x N = kM and here k is what is in the Bessel function, so is media-specific. Therefore, applying these rules, we get Bs = (-ik/kf ) E [ - b N(3)o1 + i a M(3)e1 ] // scattered field, (3) means h(1) (kr) But we know that k = nkf so we get our factor of n out front, and we pick up our "-i rule", so the result is: Bs = n E [ i b N(3)o1 + a M(3)e1 ] which agrees with our "rule" for doing this given above in the text. Now if we move ahead to the field H = B/, there is simply no way to avoid that extra factor! Bs = (n/) E [ i b N(3)o1 + a M(3)e1 ] This factor then "gets into" the two magnetic equations that arise from the n x H boundary condition. Appendix C. Regarding the Poynting Vector I have not spent much time on this, but I can see it is a real mess. What you want is this: <S> = (c/8) E x H* and outside the sphere we really have to put the sum of the plane wave and the scattered wave, so we get <S> = (c/8) { Ep + Es } x { Hp+ Hs }* where Ep = K [ M(1)o1 - i N(1)e1 ] Es = K [ - b M(3)o1 + i a N(3)e1 ] Hp = go K [ - M(1)e1 - i N(1)o1 ] Hs = go K [ a M(3)e1 + i b N(3)o1 ] So we are faced with a barrage of nasty cross products if we want to write an expression for the Poynting vector in closed form. Examples include: M(1)o1 x M(1)e1' M(1)o1 x N(3)o1 and so on. I tried computing one of the simplest ones, the M x M, and got a semi-reasonable result that of course depends on the details of the functions used. The worst are the NxN ones which are really the full product ( x M1) x ( x M2). And did I mention all the cross terms between different orders? They don't just go away because we don't have a solid angle integration here to clean things out. This seems such a mess to me that in practice if is probably easiest just to compute the four fields shown to some order in , and then directly compute the Poynting vector numerically from those field results. With this document, I think I am done now for a while with this subject.