15 compare scattering to radiation for dielectric

Radiation Approach vs Multipole Scattering Approach to Dielectric Sphere PhL 2.6,7.03 Conclusions arrived at in this little paper: (1) Reminder that the multipole expansion can give you the full exact answer to a problem as a sum of terms that you can evaluate (not integrate!) and you can reasonably hope that only some finite number of terms will be needed to get a reasonable answer. (2) It is just fine, for a dielectric radiator, to use J = dP/dt in the general full solution for A as in (9.3). (3) In the Appendix below, we use this J but in the simplified (9.13) where only the zeroth order term is kept in the exponential expansion. This gives the E1 dipole field. We compute in exact form A, , E and B for a dipole field, see the Appendix at the end of this paper. (4) For slightly extended dielectric sources, in our last paper we reviewed the Portis-like method of computing radiation patterns. These patterns involve probably all multipole orders and we suggest how much of each could order is involved. But these solutions are the first term in an expansion in smallness parameters (d/) and (d/r) and thus require both these parameters to at least be < 1, if not << 1. Thus, although all multipole orders are "hit" by the Portis-type solution, the amount of hit is not the correct amount for the full solution without the above smallness requirements. But his results are much better than the raw dipole-only result. (5) We study various realizations of the exact integrals for radiation by a dielectric, basically (9.3) in several coordinate systems. The integrals are of course not doable in closed form, even in cases with azimuthal symmetry. (6) We outline a plan for numerical integration of these formulas, and note that very many calculations are needed to get a result even at a particular observation plane z. Maybe 16 billion evaluations of eikR/R to get just a crude 100 x 100 result at one value of z. No doubt there are clever ways to speed this up using like Gaussian quadrature or FFT's. In the scattering approach, we do a boundary condition match at the sphere boundary, and this tells us the amount of scattering in each partial wave. No doubt the number of activated multipoles increases as the sphere radius increases relative to . This gives a complete and exact solution to the problem, there are no integrals to do, you just mechanically add up partial waves. You get the scattering at any angle you want, and at any distance from the sphere. For sure we have more than electric dipole action for a sphere that is not super-small compared to . Also for sure, we have more than electric dipole stimulation coming in as our plane wave. In the expansion for this wave, it keeps going up to all values of . This is most clear in the plane wave expansion on page 569. We see that in general in a given partial wave we have field components in all three directions: r, , and r x . The first and last appear in the second term. I don't think the magnitude of the components decreases at all as increases. You can see this better on page 567 where the amplitudes are the j functions in each order times P . So, the point is that we expect to get lots of action in the higher partial waves. Now, having said all that, let's look at the radiation approach of Jackson Chapter 9. We have a complete and exact formula on page 269 for A as an integral of the current density. Now what happens if we put into this thing the J shown on page 196, that is, J = dP/dt = -iP = -iE. This is the polarization current inside the dielectric. It is all over the volume, because the surface charge is running back and forth between the surfaces, so to speak. I did not think about this current until recently. You could then do some kind of expansion of the exponential such as in 9.8 page 270 for the far zone. Experiment #1. What happens if we just put J =- iP =-iE into Jackson's 9.13? The full formula applicable in this problem is, A = -ik P(r') d3r' (1) The normal lowest-order approximation that people usually do at this point is this, A = -ik (2) where one sets P(r') = E0 exp(ik1r') // exact, but in the King Smile approximation ik(R - r) ik(- r') = -i kr' if k k // the normal lowest order approximation Now you could talk about an even lower order approximation, which would be just this A = -ik P(r') d3 r' (3) Now equation (2) appears in my notes "Radiation by a Dielectric and Metal Objects", section 10. You see two reasons that the phase is going to vary. The first reason is due to the usual phase eikR in the formula, which reflects the different Huygens distances to a fixed observation point from different locations in the dielectric. The second reason is that the dielectric is being illuminated by a plane wave which causes the dielectric to have different phases at different locations even before any part of it radiates. So is there a situation where it would be reasonable to use the bare-bones approximation (3) ? Yes. That situation is when the dielectric object is very, very small compared to . In that case, you can ignore the phase exp(ik1r') of P(r'), and then you don't need the other phase either for the same reason, and that is why we have the form (3). Furthermore, if P is uniform over this little source, we can take it out of the integral to get A = -ik P V (4) where V is the volume of the source. Now we know that P arises in a macroscopic dielectric as the sum of many small dielectric objects, so P = n p where n is the density of the p objects, and P is the amount of p stuff per unit volume. We can consider a dielectric that consists of just one "large" object (but everything is still small relative to ). In that case, we just replace P V with p, the E1 dipole moment of that object. Recall that p = (r') r' d3r'. Thus, we arrive precisely at Jackson's 9.16, yet we started with Jackson's 9.13 using J = dP/dt. Tentative Conclusion. Therefore, in the general consideration of radiation from any dielectric object of any size, we can imagine an expansion of the combined exponential eikr' * eik(R-r) in some organized manner. The lowest term will have the combined exponential as unity, and this component of the radiation will give the E1 portion exactly as shown in (4) above. This is true regardless of the shape of the dielectric object, or how large it is relative to . If the object is more than a very tiny fraction of , then we expect there to be action in the higher moments. What moments are active in the Expansion (2) above? This expansion is, A = -ik (2) We have already said that when we approximate eik(R-r) ~ exp(ikr) = 1 + all remaining terms, it is the 1 term which gives the E1 contribution. Now in our Portis-like analysis, we combine the phase eik(R-r) with the phase in P, and the combined phase we then put into a form factor that was part of the answer in the case of uniform P, g(k)= (1/V) where k1 is the wavevector of the incoming plane wave, and k = r , the location of the observation point. It is important to understand the phase in this form factor is exactly the combination of the two phases discussed above. Now, if we write this phase now as 1 + something, then we can interpret the "1" term as the phase factor in the forward direction, which gives just g=1. Therefore, when we solve a problem using the above approximation, and we get an answer like this, B = k2 (1/V) x E0 g(k), // far field then if we just set g(k) = 1, what we get is the E1-only result! By the way, here is our far-field E1 result from the appendix below B = k2 x p So, the entire modulating effect of the form-factor on the dipole result comes from "higher terms". This is pretty obvious if you think about it, the form factor is unity if the phase does not vary over the volume of the radiating source. OK, I think we have made that point. Now, which higher order terms are we getting when g is present? We know that g is in general a function g(,). We can see that the B field has an extra external sin dependence. If we select a z axis along k1 as is usually done, then that is the axis to which , refer, and we know that sin = sqrt( 1 - sin2 cos2( - ') ) where ' is the azimuth of p , the direction of the assumed linear polarization of the incoming plane wave. Therefore, we could take this messy function g(,) * sqrt( 1 - sin2 cos2( - ') ) and we could integrate it against Ym(,) d to see which partial waves are active. If we ignore g and look only at the dipole sin term, we will find that Y1,1 are active due to the direction choice. But when we include the strange g(,) factor, I suspect that ALL partial waves will be affected, since the form factor is generally a very messy function like g = 2 / thin dielectric disk or like g = 3 [ - ] / dielectric sphere These messy results arise from very simple geometries. There is no dependence because the geometries chosen were symmetric about the beam axis. Normally this would restrict us to m = 0 partial waves, but the sin factor changes that. I think we would be facing this integral: Ym (,)* [ A Y11 (,) + B Y1-1 (,)] g(,) d so if g(,) has no dependence, we can think of g as a sum of Y's with m=0, and then the final results will have only m = 1 and -1. But in general, if g has an term, then the result can have +1, and -1 terms due to the sin mixing. Conclusion: Solutions to problems using (2) above generally have ALL partial waves represented. The expansion done to get (2) is not the multipole partial wave expansion, which is a rotation group thing. Rather, it is an expansion in both (d/) and (d/r), as we discussed in another paper -- both these ratios are smallness parameters of the expansion. See page 11, Section 6 of that document "radiation by dielectric...". We would draw a similar conclusion looking at the other phase factor (the one from P). Therefore, if a dielectric sphere has radius d = 4 and = 1/2 , we are not in an appropriate limit for the solution of the dielectric sphere based on equation (2). Our expansion parameter is (d/) = 8 which is not even less than 1! Having said that, we still have an exact result for the solution of our problem: A(r) = -ik P(r') d3r' = -ikP0 exp(ikz') d3r' where R = sqrt[ (x-x')2 + (y-y')2 + (z-z')2 ] . We can define the complicated factor f as, A(r) = -ikP0 f(r, k) where f(r, k) = exp(ikz') d3r' Option #1. In the integration, we can switch to spherical coordinates and choose the integration z axis to be in the plane wave k1 direction (here ) and. We would then write r' as r',',', and r as r,,, so we get R = | r - r'| = z' = r' cos' with the result f(r) = f(r,,) = exp(ikz') d3r' = r'2 dr' d(cos') d' exp(ikr' cos') Option #2. An alternative choice is to line up the integration z-axis with r, then we get R2 = r2 + r'2 - 2 r r' cos ' R = " kz' "= k r' = r' k[ sin' sin" cos('-") + cos' cos " ] In this axis choice, we think of r as being fixed, and k as moving around and being at k, ", ". It seems that this choice simplifies the messier terms of the integral. Then f(k)= r'2 dr' d(cos') * d' exp[ ir' k[ sin' sin" cos('-") + cos' cos " ] = r'2 dr' d(cos') exp[ ir' k cos' cos "] * d' exp[ ir' k[ sin' sin" cos(') ] The ' integration is a standard Bessel function, 2 J0(r' k sin' sin"), leaving a double integral, f(k) = f (k,",") = 2 r'2 dr' d(cos') exp[ ir' k cos' cos "] J0(r' k sin' sin"), The good news is that now we only have a double integral that we cannot do, but the bad news is that we get the result in coordinates that we don't really like. Note that although " = -', the relation between " and ' is not simple. That is to say, we are pretending that r (observation point) is fixed and we study the result as a function of the different angles of k. Option #3. What about cylindrical coordinates relative to the plane wave direction? A(r) = -ikP0 f(r, k) where f(r, k) = exp(ikz') d3r' In this case we get R2 = (2 + z2) + ('2 + z'2) - 2 [ ' cos (-' ) + z z' ] = (z - z')2 + 2 + '2 - 2 ' cos (-' ) R = f(z,,) = dz' exp(ikz')' d' d' Now, if the volume has cylindrical symmetry relative to the beam direction, then the limits on the d' integration are 0 to 2 and we can formally do this integration. We can shift out the dependence and replace the integral with twice 0 to . The result then becomes f(z,) = 2 dz' exp(ikz')' d' d' I don't think this integral can be done in closed form. The main point, however, is that the result of this undoable integral will be independent of . Numerical integration? To do this, you would probably just stay in Cartesian coordinates and go with this rendition having k1 along the lab z-axis: A(r) = -ikP0 f(r, k) where f(r, k) = exp(ikz') d3r' R = f(r) = f(x,y,z) = dx'dy'dz' exp(ikz') You would decide on some mesh size first of all, then set up the integral to only cover the volume of your selected object. Each triple integration over the mesh would give you one data point for f. If you had a mesh size of 100 x 100 x 100 = a million points, then you are summing a million function evaluations. Then you might pick a distance z for your observation screen, and then do 100 x 100 computations to develop the pattern on this screen. But this would only give you A. You need calculations on a few adjacent mesh planes in order to be able to compute the curl at that plane, in order to have the B field. But then you need to do the curl of THAT to get the E field. So this is a major calculation effort! Perhaps you need something like 16 planes, so we have 160,000 * 100,000 function evaluations, and this just gets you the pattern at one z. Of course if the volume is symmetric about the beam axis, it is probably much better to start with the cylindrical result above, f(z,) = 2 dz' exp(ikz')' d' d' Now if you compute this for a few planes near z , you can use the cyl curl formula without any angle terms to get the B field as B(z,), and then you need a few planes of that to get the E field. You still have a 3D integral to evaluate at each point, but we don't waste time computing angular stuff. Apppendix A. Some Exact E1 Field and Potential Calculations Computation of the E1 Scalar Potential We have, A = -ik p where p = PV // exact E1 potential Notice again that we are ignoring the spatial dependence of the incoming plane wave because we have put our little tiny object right at the origin of our coordinate system. What is the scalar potential? Since we are in the Lorentz Gauge using the above formulas for A, we know that = (1/ik)A = - ( p) = - ( ) p - ( p) = - (ik-1/r) p - ( p) = = - (ik-1/r) p Now think of this in terms of as the full integral in the equation parallel to (9.3) which has no 1/c according to page 180. The lowest order term there is going to be times the total charge. Although the second term above, - ( p), is really zero for a fixed dipole p, this term does mimic the = - P term we see on page 196, and we could no doubt derive this as Jackson did. For uniform P, though, the total charge integrates to 0, so the lowest level term gives us nothing at all! Now the leading term in the first term above comes not from the lowest term in the integral, but in the next term up where you expand ik(R - r) ik(- r')in the phase of eikR. That is to say, we do this: exp[ik(R - r)] exp[ik(- r')] 1 - ik(- r')]. The spatial integral of the non-1 term gives just p. The non-leading term above (down 1/r) is coming from expanding 1/R as 1/r (1 + r' /r) . Note that in general, the Cartesian power series expansion of the exponential is where you get the various Cartesian moments of things (as opposed to the multipole moments in the world ). So my point here is just this: although A comes entirely from the lowest zero-phase level in the integral over J, the corresponding comes from the first two phase levels of the eikR/R expansion. In fact for a constant P or a point p, the lowest zero-phase level in the integral over gives 0, and that of course would violate the Lorentz Gauge. So having the Lorentz gauge has somehow mixed in different expansion levels in the world relative to the A world. That is something I was not aware of. The Lorentz gauge applies to the total A and the total , so we are not guaranteed that they match in each level of some expansion we choose to do. They probably would match in a proper orthogonal expansion, however. Computation of exact E1 B and E Fields From the potential above that A = -ik p , we get at once B = x A = -ik x p = -ik(ik-1/r) x p and |B| = k2 | 1+i/kr| |p| sin // exact E1 results Notice that the vector x p is a vector transverse to r which is perp to the r,p plane. This is a TM mode, so we expect to find B is a vector transverse to r. The electric field is then given by -ikE = x B from Ampere's Law with displacement current, E = (i/k) x B = { (ik-1/r) } x (x p) + { (ik-1/r) } x (x p) We can do the gradient in the first term to get, E = r{ (ik-1/r) } x (x p) + { (ik-1/r) } x (x p) Now we can use these two vector results x (x p) = [ ( p ) - p ] x (x p) = -(1/r) [ ( p ) + p ] Notice that both these vectors r x (x p) and x (x p) lie in the r,p plane. The first one is an exact transverse vector to r in this plane, while the second is not a transverse vector. This gives the result E = r{ (ik-1/r) } [ ( p ) - p ] - { (ik/r-1/r2) } [ ( p ) + p ] But now Maple will do the first term derivative for us to get r{ (ik-1/r) } = [ -k2 -2ik/r + 2/r2 ] leaving us with E = [ -k2 -2ik/r + 2/r2 ] [ ( p ) - p ] - { (ik/r - 1/r2) } [ ( p ) + p ] We can now combine terms of the same vector character to get E = [ ( p ) {-k2 -2ik/r + 2/r2 - ik/r + 1/r2)} - p { -k2 -2ik/r + 2/r2 + ik/r - 1/r2} ] which we can simplify to be E = [ ( p ) {-k2 - 3ik/r + 3/r2 )} - p { -k2 - ik/r + 1/r2 } ] and then we do a rewrite as E = k2 [ ( p ) {-1 - 3i/kr + 3/k2r2 )} + p { +1 + i/kr - 1/ k2r2 } ] This is the full and exact E1 result. Let's summarize the two results right here: B = k2 (1 + i/kr) x p // exact E1 fields E = k2 [ ( p ) {-1 - 3i/kr + 3/k2r2 )} + p { +1 + i/kr - 1/ k2r2 } ] Now in the last step, we can take the large-r results: B = k2 x p // far E1 fields E = k2 [ p - ( p ) ] = - k2 x (x p) which looks like the fields of a distant outgoing plane wave, both fields are now transverse. Recompute the E1 E field using A and Another way to compute the E field would be from 6.31 on page 179: E = - + ikA where we had from above A = -ik p = - (ik-1/r) p We can compute - = { (ik-1/r) } p + { (ik-1/r) } ( p) = r{ (ik-1/r) } ( p ) + { (ik-1/r) } ( p) Now we know that ( p) = (p ) ( ) + p x ( x ) = (p ) ( ) = (1/r) [ p - ( p ) ] which then gives - = r{ (ik-1/r) } ( p ) + { (ik/r-1/r2) } [ p - ( p ) ] We then have as our final result E = - + ikA = r{ (ik-1/r) } ( p ) + { (ik/r-1/r2)} [ p - ( p ) ] + k2 p = [ -k2 -2ik/r + 2/r2 ] ( p ) + { (ik/r-1/r2)} [ p - ( p ) ] + k2 p and now we group terms = [ ( p ) { -k2 -2ik/r + 2/r2 - ik/r +1/r2 } + p { k2 + ik/r-1/r2 } ] = [ ( p ) { -k2 -3ik/r + 3/r2 } + p { k2 + ik/r - 1/r2 } ] = k2 [ ( p ) { -1 -3i/kr + 3/k2r2 } + p {+1 + i/kr - 1/ k2r2 } ] which we can compare to an intermediate result above from the other method, E = k2 [ ( p ) {-1 - 3i/kr + 3/k2r2 )} + p { +1 + i/kr - 1/ k2r2 } ] As hoped, the results are the same either way you compute them. The method was not any easier than the first method. At least this confirms our expression for .