Polarization of sky light from a canopy atmosphere

Light from the clear sky is produced by the scattering of unpolarized sunlight by molecules of the atmosphere and is partially linearly polarized in the process. Singly scattered light, for instance, is fully polarized in viewing directions perpendicular to the sun direction and less and less so towards the parallel and antiparallel directions, where it is unpolarized. The true, multiple, scattering is much less tractable, but importantly different, changing the polarization pattern's topology by splitting the unpolarized directions into pairs. The underlying cause of this ‘symmetry breaking’ is that the atmosphere is ‘wider’ than it is deep. Simplifying as much as possible while retaining this feature leads to the caricature atmosphere analysed here: a flattened sheet atmosphere in the sky, a canopy. The multiple scattering is fully tractable and leads to a simple polarization pattern in the sky: the ellipses and hyperbolas of standard confocal ellipsoidal coordinates. The model realizes physically a mathematical pattern of polarization in terms of a complex function proposed by Berry, Dennis and Lee (2004 New J. Phys.6 162) as the simplest one which captures the topology.


Introduction
Rayleigh scattering of unpolarized sunlight by molecules of the atmosphere renders the sky light partially linearly polarized. Standard simple geometrical arguments suffice to show that the singly scattered light is fully polarized in viewing directions perpendicular to the sun direction and unpolarized parallel to it. To account for light that is scattered twice or more is considerably harder, but important in that it produces new topological features in the pattern of polarization direction in the sky. No longer is the sun direction (and its opposite) unpolarized; instead there are twin directions of unpolarized light, above and below the sun direction, as was observed by Arago and Babinet and Brewster in the early 1800s. The underlying cause of this symmetry breaking is that the atmosphere is (much) less deep than it is 'wide'.
Here, to capture this feature in caricature, the atmosphere is imagined to be flattened to a plane sheet in the sky, or canopy (figure 1). This model is analysed, including all numbers of scatters. The pattern of line segments covering the sky hemisphere of viewing directions, indicating the orientation of maximum polarization at each point, is shown to generate the curves of standard confocal ellipsoidal coordinates. Equivalently, the model realizes physically the complex function form for the pattern of polarization proposed by Berry, Dennis and Lee [1] as the simplest mathematical one which captures the topology. As they point out, the same pattern applies to a related problem in crystal optics: finding the polarization eigendirections for a biaxial (transparent, non-chiral) crystal as a function of the propagation direction.
Multiple scattering in a three-dimensional atmosphere, for example one with exponential decay of density with height, and having absorbtivity proportional to density, is complicated, even for a flat Earth [2,3]. Although hard to solve, the mathematical problem of 'radiative transfer' is easy enough to formulate as follows. A light ray from the sun to the observer's eye is a sequence of straight segments between discrete scatters by air molecules, while subject to continuous absorption on its route. The initial segment has fixed direction from the sun, treated as a point source of unpolarized light at infinite distance, and the final segment has its direction fixed to be in the direction in which the observer is looking (or rather, the opposite direction). All possible scattered rays are allowed.
Partial polarization of the light and its Rayleigh scattering is also easy to visualize geometrically and formulate mathematically. The excitation of an atom in the sun is a statistically isotropic vibration of dipole moment in three dimensions and the electric field is likewise. The Geometry of scattering from a top view. The bent 'drainpipe' shown is a schematic way of visualizing the geometrical (part of the) reduction of electric field amplitude of a single scattered ray. The three middle legs lie in the canopy scattering sheet, and the first and last are the sun and observer legs. The sun leg has a circular cross section representing the unpolarized electric field fluctuations and terminates at the first scatter with a perpendicular cut. The new leg of pipe has the necessary elliptical cross section (partial polarization) and oblique cut to join the sun leg. It is cut perpendicularly at the second scatter and so on. (At the third scatter shown, the new pipe emerges somewhat backwards, passing through the old: unphysical for a drainpipe, but perfectly acceptable here.) sunlight arriving at the Earth's atmosphere is a projection of this isotropic ensemble along the propagation direction, that is, an ensemble of transverse electric field vectors, isotropic in two dimensions, representing unpolarized light. It can be specified mathematically by a 3 × 3 tensor: the time average of the product E i E j of electric field components, E x , E y , E z (or of |E E| in the notation to be used below).
Ignore absorption for the present. At first scatter (figure 2), this ensemble is simply projected along a new direction parallel to the line to the second scatter, making a new anisotropic ensemble in the two transverse dimensions. There is also a reduction factor from the scatter. This is repeated until the ray enters the eye. A consequence of the entirety of scattered rays is that the sky light is partially polarized. For each viewing direction, looking through a transverse polaroid, rotatable in its own plane, there is a brightest orientation, and orthogonally a darkest one. These are the two 'eigendirections' of the 3 × 3 matrix associated with the light entering the eye which are perpendicular to the viewing direction, the third being parallel. The absorption factor for the whole ray is exp(−constant × density of the atmosphere integrated along the ray path). All this can be cast into a local, differential equation format, which is Chandrasekhar's method, although with less explicit geometry. I shall not follow this here, but retain the global picture, and apply it to a flat sheet, canopy atmosphere. Having two-dimensional geometry and constant absorbtivity renders the scattering calculation quite straightforward, but there is a separate reason, independent of the details, why the answer obtained has a simple mathematical form. This is explained below, before the scattering calculation in the next section.
For the flat Earth model with a general atmosphere, either three-dimensional or canopy, the atmosphere is necessarily horizontally uniform at each height level. Also, the sun illuminates the atmosphere uniformly. Therefore the statistical electric field (3 × 3 tensor) experienced by any molecule at a particular height level, resulting from all rays arriving there with all numbers of scatters, is also necessarily uniform (and so is its consequent statistical dipole excitation). Now, with the canopy atmosphere model, the observer, looking up in different directions, sees different molecules, but all with the same excitation tensor. It is the different projections of a fixed 3 × 3 tensor which is responsible for the simple analytic form of the partial polarization observed for a canopy atmosphere. One might imagine that one could use the same argument for a three-dimensional atmosphere, applied to each layer separately. The feature spoiling the simplicity in that case is that there is a different absorption of the light through the layers of atmosphere below the contributing molecule (the final scatter), so even molecules in the same layer are not only observed via different projections (as for the canopy case) but with different strengths.
As a matter of notation preparatory to the scattering theory, it is convenient to denote column and row vectors differently, as |a and a|, respectively, each being the same set of three separate real numbers (no quantum connotation is intended; this Dirac bracket notation is practical more widely). Then the symbol a | b denotes the usual scalar dot product of the vectors, and the symbol |a b| denotes the 3 × 3 matrix of product entries a i b j . Importantly, if a is a unit vector, then |a a | b is the column vector formed by resolving b along the a direction, and so |b − |a a | b is the vector formed by resolving b perpendicular to the a direction. This latter can be written [I − |a a|] |b with I as the 3 × 3 identity matrix. The bracket term is thus the matrix for projecting along the a direction, that is, yielding a vector perpendicular to a.

Scattering theory
Denote the unit vector of a ray of light from the sun by the column vector |S , those of the scattered legs of the ray by |1 , |2 , |3 , etc and, finally, that from the last scatter to the observer by |O . Let the lengths of the these legs be L S , L 1 , L 2 , L3, . . . , L O . The partial polarization of this ray arriving at the observer is described by a 3 × 3 matrix: the time average |E E| of the electric field product |E E|. The electric field |E is generated by multiplying successive projection matrices [I − |k k|] for the legs k = 1, 2, 3, etc. With α as the intensity scattering factor for each of the n scatterings, the electric field vector arriving at the observation point from the specified ray is where E S is an isotropic random vector, meaning that the time average of |E S E S | is 1 3 E 2 S I, representing an effective source at the sun.
The contribution from the ray in question to the received polarization matrix |E E| is thus (the product of n infinitesimal integration areas multiplied by) It may be noted incidentally that this matrix is of course purely real, signifying that the partial polarization is linear (that is, the light is a superposition of unpolarized intensity with fully linearly polarized intensity, not fully elliptically polarized [4]). This property is shared by any model of the atmosphere.
To account for all rays, the directions |O and |S are taken as fixed (although the latter vector is free to translate to allow any position of first scatter), and one needs to integrate the expression over all the n − 1 vector displacements (L 1 |1 ), etc lying in the scattering sheet, the x-y plane, say, each being unrestricted in length. The radial L integrations can all be done separately and each yields a factor β. This will depend on the absorbtivity specified and (logarithmically) on the inner cut-off imposed to account for the fact that scattering molecules cannot overlap. Likewise, the direction of each unit vector 1, 2, 3 . . . can be integrated in order, with each giving a result in terms of the last. So, starting with the sun in the x, z plane at elevation angle σ above the horizon, we have and using where φ k is the angle of k with the x-axis, we can calculate the integrals of the form in question (which stays the same as k changes). Thus, for any matrix element a, b, c, d (representing the outcome of all k − 1 integrations so far), For n scatters, we therefore have the polarization matrix Summing over all possible numbers of scatters n (n = 1, 2, 3, . . .) using geometric series (including the matrix series for this last 2 × 2 matrix: , we obtain the final result for the polarization matrix received by the observer. As anticipated, it is a central 3 × 3 tensor independent of the observation direction, sandwiched between projectors along the observation direction. with γ = 2παβ/8. The principal orientations of partial polarization are the two eigendirections of the result (7) (the third being |O ). For a polaroid held across the viewing direction |O , these give its orientations of maximal and minimal transmission. Importantly, there are two observation directions |O (or their opposites −|O if the observer was above the canopy) for which the eigenvalues are equal and the eigendirections consequently indeterminate: the 'unpolarized' and 'neutral' directions. For other viewing directions |O , the eigendirections of the matrix (other than |O itself) can be presented as line segments tangent to the unit hemisphere of viewing directions ( figure 3(b)). To obtain their latitude (θ) and longitude (φ) components explicitly, one rotates coordinates with a rotation matrix R has, by design, the essentially 2 × 2 form required These results can be expressed more elegantly by specifying curves on the hemisphere of viewing directions to which these eigendirections are tangent. This re-expression can be done in either of two ways. Recasting the problem in terms of a cubic equation will show that the pattern of eigendirections described is that of tangents to confocal ellipsoidal coordinate curves with consequent simple formulae for the curves. Alternatively, recasting the problem in terms of a quatic function will confirm the equivalence with the simple complex function form for these eigenvector directions, proposed as a purely mathematical model by Berry et al [1], with its consequent elliptic integral for the curves. These reinterpretations will be described in turn. For the sun at non-zero elevation, the polarization picture loses the left-right symmetry, but a suitably tilted confocal ellipsoidal coordinate system can be found with singular directions matching the new unpolarized directions, which suffices for matching everywhere.

The cubic representation
An analogous problem in mechanics was solved by Binet [5] in 1813, and it is Binet's 'remarkable theorem', as reported by Thomson and Tait [6] which I now describe. The task is to find the eigendirections of the inertia tensor of a fixed rigid object about different points in space. Binet showed that the eigendirections are exactly the three coordinate directions of confocal ellipsoidal coordinates ( figure 3(a)). For points on an infinitely large sphere centred on the object, the problem will be seen to reduce to that of (7) (with the third eigendirection radial). The coordinate system will be described first [7].
Confocal ellipsoidal coordinates λ, µ, ν are the three roots of the cubic equation in ξ, where A 2 , B 2 and C 2 are constants (which will represent the principal moments of inertia divided by the object mass). The solution relating x, y, z to λ, µ, ν is and similarly for y and z, cycling A, B and C. For a fixed, infinitely large value of the largest root λ, say, the point (x, y, z) lies on an infinitely large ellipsoid which has become spherical, x 2 + y 2 + z 2 = λ, and the curves of constant µ (next lower root) are the intersections of the sphere with hyperboloids of one sheet. Those of constant ν are the intersections of the sphere with hyperboloids of two sheets. (These curves will be tangent to the eigendirections in the applications.) Eliminating ν between (12) and its y 2 counterpart (and neglecting A 2 and B 2 compared to λ) gives where A 2 > B 2 > − µ > C 2 and the same equation applies with µ replaced by ν where A 2 > − ν > B 2 > C 2 . This concludes the description of the coordinate system. To relate this to mechanics, let G be the inertia tensor divided by the mass of the object (the 'gyration tensor'), about its centre of mass. The tensor [8] about any other point r is G + r 2 I − |r r|. If ellipsoidal coordinates based on the eigendirections of G are set up, it is straightforward to show that the coordinate directions at r are the eigendirections of the modified tensor there. To see that this problem is equivalent to (7), one recognizes that, as r → ∞, the one eigendirection with non-infinite moment of inertia eigenvalue, mass times r|G|r is radial, so The other eigendirections are thus those of [I − |r r|]G[I − |r r|], with eigenvalues equal to those of this matrix plus r 2 . This, having the form (7), verifies the equivalence of the problems. In the present sky context, the confocal ellipsoidal coordinates described and depicted in figure 3(a) can be directly matched to the polarization orientations in the sky ( figure 3(b)). The case shown is for the sun on the horizon. All that is necessary is to match the pair of singular coordinate directions with the pair of unpolarized directions in the sky determined by the central matrix in (7). If the confocal ellipsoidal coordinates are projected down onto the horizontal plane, they form the ellipses and hyperbolas of equation (13), and the polarization pattern of line segments, likewise projected, are tangent to them. As the elevation angle of the sun increases above zero, the angle between the pair of unpolarized directions changes somewhat (and therefore the angle between the singular directions needs changing similarly), and both pictures need tilting about an axis coming out of the page. The simple ellipses and hyperbolas are now formed by perpendicular projection onto the tilted equatorial plane, instead of the horizontal one.

The quartic representation
An alternative recasting of the problem in terms of a complex quartic proceeds as follows. The sky is here represented not as a hemisphere of viewing directions but as the stereographic projection of this hemisphere onto the equatorial plane (the projection point being the South Pole). Coordinates on this plane are taken as the real and imaginary parts of a complex variable ζ. The central matrix in (7) can have any multiple of the identity added to it without affecting the eigendirections of the product matrix. Therefore, to render it traceless one third of its trace multiplied by the identity is subtracted. Any traceless symmetric 3 × 3 matrix (i.e. five numbers) is naturally associated with a quartic polynomial (five coefficients) [9,10]. If the matrix is real, as in our case, the coefficients are in pairs (±complex conjugates). The association, for a quartic polynomial Q(ζ), is expressed in terms of 2 × 2 matrices (deriving from Pauli matrices using premultiplication by 0 1 −1 0 ): where the 3 × 3 symmetric traceless matrix is |a b| + |b a|− 2 3 Any 3 × 3 symmetric traceless matrix M can be represented [11,12]  With the ensemble interpretation just stated, these are evidently the unpolarized directions. As explained by Berry et al [1], and following Berry et al in earlier papers [13,14], the phase of √ Q(ς) in the complex plane supplies the orientation of (one of) the eigendirections of the received polarization (7) stereographically projected onto the plane. The vector , which is the gradient of a real function whose contours are tangent to an eigendirection. Likewise, its orthogonal is along the gradient of a companion function, the two functions being the real and imaginary parts of the complex elliptic integral Q(ς) −1 dς.

Final remarks
Of the two interpretations given, one (the cubic) supplies curves in terms of elementary functions on the sphere and the other (the quartic) does so in terms of an elliptic integral function on the stereographic plane. These need reconciling since an elliptic integral function is not elementary (and stereographic projection cannot make it so). The resolution of the apparent conflict is that area elements corresponding to confocal ellipsoidal contours at equal intervals of µ and ν are not square but rectangular. They need re-labelling to become square and hence admit description by a complex function (by the Cauchy-Riemann equations), and it is the relabelling that involves the elliptic integrals. In fact [7] if µ is replaced by the elliptic integral 1/ (A 2 + µ)(B 2 + µ)(C 2 + µ) dµ, and similarly for ν, the area elements do become square. Perhaps the most serious defect of the thin sheet atmosphere is that it gives infinite intensity of sky light from the horizon (irrespective of the elevation of the sun). This infinity is not contained in equation (7) which gives the polarized intensities observed per molecule. Indeed equation (7) gives zero intensity on the horizon since the factor 1/L 2 O is proportional to (sin θ) 2 where θ is the viewing direction elevation angle. But the number of molecules per unit solid angle is proportional to 1/(sin θ) 3 so the product, 1/sin θ, diverges unphysically. This would be exactly rectified by Lambert's cosine law, but this is unjustified within the model. As remarked earlier, one may lay the blame for the complexity of the three-dimensional atmosphere polarization pattern entirely on the absorption in the final leg which complicates the polarizations received from the different directions. But it is also this absorption which cures the problem of the unphysical infinite intensity from the horizon, so it cannot be lightly dispensed with.
As far as the canopy atmosphere is concerned, perhaps the least objectionable cure for the horizon infinity would be to invoke the curvature of the true Earth. One could reasonably argue that a curved canopy sheet atmosphere at a height of 8 km, say, was flat enough for the above scattering analysis to supply the (nearly) constant excitation tensor in the curved sheet. Only the geometry of the observer's view is different, in particular being truncated by the horizon. Light from the horizon, which lies at a distance of about 300 km for the 8 km canopy height chosen, is emitted from the canopy at a small but non-zero angle, avoiding the infinite intensity.