## Abstract

Daylight, or sky light, is sunlight Rayleigh scattered by the atmosphere onto the ground. This random scattering propagation through clear air is governed by ‘radiative transfer’. Beyond the single-scattering approximation, the famous virtuoso analysis of Chandrasekhar formulated the problem and offered exact, but rather involved, and ultimately numerical, algorithms for its solution. However, there is no real difficulty in writing down directly the exact Rayleigh scattering series in integrals. Its practical utility is limited to fairly small thicknesses *T* of atmosphere (compared with the mean free path), but the Earth has just such. Here even the next order beyond single scattering (error order *T*^{2}) supplies a formula for the brightness and partial polarization of daylight across the sky, which captures the essential topology of the polarization pattern, and also remains uniformly valid in the small thickness limit, for all elevations of the Sun and viewing angles. The status of the mathematical polarization direction pattern invented by Berry, Dennis and Lee as the simplest fit to the required topology is clarified.

## 1. Introduction

The standard idealization in the treatment of scattering in a clear atmosphere, which is well justified, is that of an infinite flat horizontal Earth with an atmosphere above it with density depending on height only (‘stratified’). The incident light is from an infinitely remote point source in some direction, i.e. it is an infinite parallel beam, oblique to the vertical in general. In consequence, the problem has horizontal translational invariance; there is no horizontal dependence of any property of the light. Because the air molecules have almost uncorrelated positions (or equivalently the density fluctuation correlation function is zero), scattered intensities add, no coherence or interference is involved and the ‘radiative transfer’ description is appropriate. The light can be considered as a beam of independent rays passing like particles through a continuous scattering medium. In addition, because the air molecules are isotropically polarizable by the electric field of the light, the scattering is Rayleigh scattering. This is ‘elastic’—there is no absorption, a beam of light is attenuated only by scattering out of it.

Radiative transfer has a remarkable independence, or equivalence property explained in appendix D: the vertical density profile of the atmosphere is immaterial. An exponential profile atmosphere, for instance, cannot be distinguished from a uniform slab atmosphere (if they have the same total number of air molecules *N* above each unit area of ground). Every feature of the emergent daylight is the same. All that matters is the total attenuation thickness *T* of the atmosphere, a dimensionless number, defined by the transmission fraction exp(−*T*) of the whole atmosphere for a beam of normally incident light. (*T*=*Na*, where *a* is the scattering cross section of an air molecule).

The value of *T* for the Earth's atmosphere is quite small, approximately 0.15 for the blue light which predominates (because the molecular scattering strength increases rapidly with frequency; as the fourth power, derived by Rayleigh (1871)). The refinement, to examine the colour spectrum of daylight by treating each infinitesimal band of frequency separately with its own value of *T*, will not be considered here, though it is probably necessary for extremes like sunset. The small value of *T* suggests that the familiar single-scattering approximation should be good. Although it is often good enough, sensitive features such as the direction of partial polarization are not adequately described. Single scattering also fails, even for robust features such as brightness, if the Sun and viewing direction are both close to the horizon—single scattering is not the ‘uniformly valid limit’ as *T* tends to zero.

Since there is no loss of generality in specializing to a particular model and the uniform slab seems the simplest, that model will be adopted here. An interpretation of *T* is then the ratio of the thickness of the slab to the mean free path of light in an infinite medium of that uniform density.

There are two equivalent approaches to the mathematical formulation: continuous or discrete. In the continuous formulation, extensively used and developed by Chandrasekhar (1960), integro-differential equations express the propagation of parallel beams in all directions. In the discrete formulation (which can be considered in its own right, or as an iterative solution of the continuous formulation, as described in appendix F), each individual ray of the beam is imagined (figure 1) to undergo a succession of scatterings at points (molecules) with straight line propagation of appropriately attenuating intensity in between. All rays of the incident beam reaching an observation point on the ground with all possible numbers of scatterings must be accounted for, and distinguished according to their direction of arrival. This discrete formulation will be used here.

There is no real difficulty in writing down the exact Rayleigh scattering series, and for the Earth's small atmosphere thickness *T* truncating this series can supply simple approximate formulae. In contrast, for the equations of the continuous formulation, rather involved manipulation casts the exact solution in terms of infinite products, requiring, ultimately, numerical evaluation. Much of the early work (King 1913; van de Hulst 1948), as well as most of Chandrasekhar's book (1960), ignored polarization and treated a simplified ‘scalar scattering’ situation more appropriate to sound waves than light. This remains a valuable initial exercise (§2) because the integrals for the scattering series are largely the same, and the visualization is much more straightforward. Indeed visualization for Rayleigh scattering is usually suppressed, but since it plays a key role in the present analysis, the next paragraph presents a rather fanciful (and optional) preliminary interpretation. Thereafter, the algebraic description begins.

Rayleigh scattering takes into account the partial polarization of the light. A ray line can be imagined to be an infinitely narrow tube like a drainpipe tube, but with an elliptical cross section (which represents its ‘partial polarization’—unpolarized light would have circular cross section). At a scattering, the incoming drainpipe is cut off perpendicular to its axis (this cut plane is the plane of electric field fluctuation at the scattering molecule due to the ray in question). The outgoing section of drainpipe has a new elliptical cross section depending on its direction—that unique cross section which will join on perfectly to the cut end of the old drainpipe, with an oblique cut which is also fixed by the leakproof perfect joining condition (physical drainpipes would be restricted to in–out angle >*π*/2, but waive that restriction). The fraction of the ray's light intensity which is polarized is the difference in area of the exscribed and inscribed circles of the ellipse divided by their sum. The manner in which contributions of the (continuum of) different ray paths with their ellipses combine is via the polarization matrix described below and in appendix A.

The partially polarized light seen can be described by a 2×2 symmetric real matrix (reduced from a 3×3 one in which the analysis is naturally carried out) which depends on observation direction (elevation angle *θ*, and azimuth *ϕ* with respect to the Sun). The basis directions for the matrix are along the local *ϕ* and *θ* coordinate directions (latitude line and longitude line of the observation direction). The trace of the matrix is the brightness seen (the sum of the eigenvalues). The polarized fraction is the difference in the eigenvalues divided by their sum, and the polarization direction is the eigendirection of the greatest eigenvalue.

A striking general consequence emerges from the fact that at every point in the atmosphere the electric field fluctuations are symmetric about a mirror plane parallel to the Sun's meridian plane. The consequence of this fact (within Rayleigh scattering), derived in §4, is that if the 2×2 polarization matrix is known for any one direction (*θ*, *ϕ*) of some elevation *θ*, then it can be directly deduced for all directions of that elevation; the *ϕ* dependence is explicit from the beginning. One has exactly(1.1)where *α*, *β*, *γ* and *δ* are *independent* of *ϕ*, but depend on *θ* as well as the two parameters of the problem, the elevation angle of the Sun and the attenuation thickness *T* of the atmosphere. The more scatterings that are taken into account, the more accurate are the values for *α*, *β*, *γ* and *δ*. Pictorially, the resulting polarization and brightness can be displayed as a field of crosses (figure 3 for the main formula (4.5)). Noticeable features are the brightening towards the horizon, and the two directions in the sky of unpolarized daylight (‘neutral’ directions), rather than a single one in the direction of the Sun (Coulson 1988).

This last feature of daylight was discovered observationally in the 1800s and was taken as a starting point by Berry *et al*. (2004) for the invention of a mathematical function (a complex quartic), which, given only the two neutral directions, surrounded them with the simplest natural pattern of polarization directions across the sky. There seemed to be good agreement with observed field of directions, and with a fairly natural extension of the quartic function, a polarized brightness (total minus unpolarized) field across sky was proposed. (The total brightness field was beyond the scope of the invention.)

To try to account physically for this mathematical function, and to complete it with a specification of the total brightness field, I proposed (Hannay 2004) a simple ‘canopy atmosphere’ model with values of *α*, *β*, *γ* and *δ* in (1.1) which are ‘constants’—independent of *θ* as well as *ϕ*, depending only on the Sun elevation and *T*. As I showed, the constancy is all that is needed to reproduce exactly the Berry–Dennis–Lee pattern of directions. In addition, if the 1/sin *θ* prefactor in (1.1) is disregarded, the constant *α*, *β*, *γ* and *δ* model reproduces exactly the polarized brightness predicted by the Berry–Dennis–Lee extension of the quartic. Unfortunately, the 1/sin *θ* prefactor, being present for a compelling physical reason, cannot be disregarded, and leads to an unphysical infinite brightness on the horizon. This was already lamented in Hannay (2004), but the underlying flaw in the model was then unclear. Appendix E explains this, and supplies instead an exactly solvable ‘mean field-like’ model, with no blazing horizon. The new model retains the Berry–Dennis–Lee pattern of polarization directions (but sacrifices the extension to polarized brightness). Albeit rather artificial, this model is not unphysical like the canopy atmosphere.

The main task, though, is the extraction of the exact Rayleigh scattering series, and its truncation to produce an approximate formula for *α*, *β*, *γ* and *δ*, i.e. for the brightness and partial polarization of daylight. It will indeed turn out that for elevation angles of Sun and viewing (*θ*) that are not too small (bigger than about *T*), *α*, *β*, *γ* and *δ* turn out to be nearly constant, explaining the success of the Berry–Dennis–Lee pattern of polarization directions over most of the sky at most times of day. By the same token, though, their extension to polarized brightness is too small by the factor sin *θ*.

## 2. Scalar scattering series

As a simplified preparatory analysis in this section, the polarization of the light is ignored temporarily so that instead of the correct ‘Rayleigh scattering’, there is the scalar isotropic three-dimensional scattering of rays mentioned in §1. The procedure for the true Rayleigh scattering later will closely parallel that for the scalar scattering, and very little further calculation will be required. As described, attention can be restricted to a uniform slab atmosphere without loss of generality. The slab is fully characterized by its attenuation depth *T* introduced earlier (normal transmission fraction exp(−*T*)). Let the parallel incident rays of light from the Sun have a direction cosine with the downward vertical *μ*_{s} (the sine of the elevation angle of the Sun). There follow then a number of simple results which will be needed.

Considering one such ray, the probability that it remains unscattered by the atmosphere of attenuation thickness *T* is exp(−*T*/*μ*_{s}). The probability that it remains unscattered at least to an attenuation depth *t*(<*T*) is correspondingly exp(−*t*/*μ*_{s}). The probability that it is scattered for the first time between *t* and *t*+d*t* is therefore . The vertical bar denotes that *μ*_{s} is a given parameter, in the manner of conditional probabilities, and the symbol *P* will be used for various different probability densities distinguished by their types of argument. Given such a scattering at depth *t*, the probability that the scattered ray hits the ground without further scattering and with the sine of its elevation angle *θ* between *μ* and *μ*+d*μ* is a different probability density (any component of a uniformly random direction unit vector being itself uniformly random between ±1, with, therefore a probability density 1/2). This singly scattered light as well as all subsequently scattered light shares the property of azimuth independence, meaning that the sky brightness from this scalar model (unlike the true Rayleigh scattering later) is rotationally symmetric about the zenith, wherever the Sun is. The corresponding probability that instead the light exits the top of the atmosphere (*μ*<0) without further scattering is . The sum of the integrals over *μ* of these two expressions (each over the range 0<|*μ*|<1) is the total probability, given a first scattering at *t*, that there will be no further scattering. Thus, for instance, the probability of at least one scattering and the probability of exactly one scattering are respectively(2.1)The latter probability can also be evaluated analytically, but is rather involved and is not needed.

To proceed to obtain the exact scattering series (2.6) for the distribution of scattered light the last quantity required is the probability *P*(*t*′|*t*)d*t*′. This is the probability that, given a scattering (which may or may not be the first) at attenuation thickness *t*, there is a next scattering between attenuation thickness *t*′ and *t*′+d*t*′. (Each scattering may be anywhere horizontally.) It is easily evaluated exactly, recalling azimuth independence. The probability density d*μ*/2, for *μ*, should be integrated over the interval −1 to 0, or 0 to 1, depending on the sign of (*t*′−*t*), but equivalently by sign change of the integration variable,(2.2)where *E*_{1} is the standard ‘exponential integral’ function. The related functions *Ei* and *E*_{n} (*n*=positive integer) will be required later, the definitions being(2.3)It is worth noting for later that the series expansion of *E*_{1} and *Ei* lead with a log term followed by a normal power series. The exact expression for the scattered (scalar) light as a scattering series of integrals is as follows. The probability that the incident ray emerges at the bottom of the atmosphere with the sine of its elevation angle between *μ* and *μ*+d*μ* is *P*(*μ*|*μ*_{s})d*μ*, where(2.4)By independence of azimuth *ϕ* mentioned earlier, the probability of exiting in the solid angle *μ* to *μ*+d*μ*, *ϕ* to *ϕ*+d*ϕ* is just *P*(*μ*|*μ*_{s})d*μ* d*ϕ*/2*π*.

The same formula, with the altered definition of *P*(*μ*|*t*) described, applies to upward scattering back into space. It is easiest to concentrate henceforth on the downward scattering (daylight) and state the necessary modification for backscatter, namely multiply the answer by exp(*T*/*μ*) and then replace *μ* by −|*μ*|.

There is an important ‘obliquity factor’ *μ*_{s}/*μ* by which *P*(*μ*|*μ*_{s}) should be multiplied to determine the brightness seen by an observer (in this scalar scattering model, but the same statement holds for the proper Rayleigh scattering later) for the following reason. The result *P*(*μ*|*μ*_{s})d*μ* d*ϕ*/2*π* supplies the probability of exit at the bottom surface of the atmosphere, in any element d*μ* d*ϕ* of solid angle, given one ray entering at the top surface. For an oblique parallel beam of incoming sunlight, the ‘rate of entry’ per unit area of the horizontal top surface is the rate through a smaller area *μ*_{s} *perpendicular* to the beam. In addition, after exit, what is seen is the rate per unit perpendicular area, rather than the rate per unit horizontal area supplied by *P*(*μ*|*μ*_{s}). The former rate is 1/*μ* times the latter rate. Therefore, precisely stated: the rate out per perpendicular area per solid angle is (*μ*_{s}/*μ*)(*P*(*μ*|*μ*_{s})/2*π* times the rate in per perpendicular area. This ratio is what one may call the brightness seen.

The first two terms of the scattering series (2.4), corresponding to single and double scattering can be evaluated analytically ((2.5) and (2.7)), but higher numbers of scatterings cannot, in full. They can, however, each be expanded as a separate series in the attenuation depth *T*, and the leading two orders in *T* for each can be found analytically. Apart from the single scattering, which is free from log terms, the *n*-fold scattering term starts with (log *T* being large in magnitude for *T* small). This is followed by all the terms with *n*′>*n* and *n*″<*n*. Table 1 displays this, continuing indefinitely downwards and to the right.

Blank entries are not present in the series, ticked entries can be evaluated analytically and crossed ones cannot. Single scattering as a whole is explicit, as is double scattering. Thereafter, the leading two entries for each *n* are accessible but no more. Thus, the first entry in the table, which cannot be evaluated analytically, has order *T*^{3}, so this is the error incurred by discarding it and higher-order entries (including terms of single or double scattering of order *T*^{3} or higher, at one's discretion).

The single scattering (first term of (2.4)) can be evaluated straightforwardly.(2.5)The multiple scattering integrals (*n*>1) are(2.6)The double-scattering term, *n*=2 can be evaluated in terms of exponential integrals (van de Hulst 1948)(2.7)The exponential integral *Ei*, which appears once in each brace brackets, was defined in (2.3). The triple-scattering integral cannot be evaluated analytically, but the leading two orders of it can, and this is done in appendix B (the same applies to higher numbers of scatterings). This completes the analytical approximation for scalar scattering: from table 1; the single-scattering error is *T*^{2} log *T*; including double scattering the error is *T*^{3}(log *T*)^{2} and the error including the leading two orders of triple scattering is *T*^{3}.

Provided that not both *μ* and *μ*_{s} are small (comparable to *T*), single scattering dominates (figure 2) and gives a brightness for daylight (in the scalar model) which has only a small correction from higher scatterings. In fact provided neither *μ* nor *μ*_{s} is small (more restrictive), it yields nearly a constant, *T*/2*μ*_{s} independent of *μ*. When multiplied by the obliquity factor *μ*_{s}/*μ* described after (2.4), it represents (in the scalar model) a sky brightness increasing towards the horizon (till *μ* becomes comparable with *T* and the formula no longer yields a constant). If both *μ* and *μ*_{s} are both small, comparable with *T*, single scattering does not dominate (figure 2), and therefore single scattering is not the uniformly valid limit as *T* tends to zero for all elevations of the Sun and all viewing directions. To rectify this, and obtain a uniformly valid approximation, at least the leading order of double scattering needs to be included (and this leading order suffices). The same argument applies to the case of Rayleigh scattering later.

## 3. Rayleigh scattering series

Electromagnetic radiation is not scalar in nature, but vector. Daylight is polarized, partially, by Rayleigh scattering from air molecules. Light seen by an observer looking at a point in the sky is characterized not only by brightness but also by its unpolarized fraction, together with the direction (angle about the line of sight) of the linear polarized fraction (there is no component of circular polarization in Rayleigh scattering). The intensities of the polarized part and the unpolarized part add; they are mutually incoherent.

The way in which the partially polarized light is produced by the sequence of scatterings is described in some detail in appendix A. The outcome is that as in the scalar model all possible scattering paths (from the Sun to the observer's eye along the line of sight) are integrated over, each with its probability as before, but now also each scattering path has an associated real symmetric matrix depending on its shape. (If a non-probabilistic formulation is preferred, the iteration in appendix F generates the Rayleigh scattering series.) The matrices weighted with the path probability are added (i.e. integrated over all possible scattering paths) to obtain a resulting (partial) polarization matrix ** P**. The brightness and polarization seen along the chosen line of sight are obtained from

**: its trace is the brightness, the eigendirection with the greatest eigenvalue is the direction of polarization, and the unpolarized fraction is 2(1−**

*P**greatest eigenvalue*/

*Trace*). These statements are true irrespective of which version (the 3×3 or 2×2 reduction) of the matrix

**is used. The expansion of each number of scatterings in the attenuation thickness**

*P**T*is governed by the scalar path probability, and the orders of the terms (table 1) are the same for Rayleigh scattering as for scalar scattering. For the present section, formally deriving the scattering series, the more natural and global 3×3 version of

**will be used. It will be reduced to the more economical 2×2 form for the formulae in §4.**

*P*The matrix describing the unpolarized unit intensity incident sunlight can be written down. Properly described, it is the ‘point coherence matrix’; the result of integrating the ** P** over the solid angle of the Sun—ideally a delta function. As described in appendix A, it is the 3×3 matrix of unit trace, (1/2)

*M*_{S}, where

*M*_{S}is the projection matrix ,

**being the identity matrix; and is the outer product (**

*I**column vector*⊗

*row vector*) of the unit vector in the Sun direction. Let

**denote the projection matrix along the line of sight, unit vector**

*M***=(cos**

*O**θ*cos

*ϕ*, cos

*θ*sin

*ϕ*, sin

*θ*) and

*M*_{j}is the projection matrix along the

*j*th intermediate segment of a scattered ray path (

*M*_{1}is written out explicitly in (3.4)). The matrix which is to be weighted by the path probability to yield the polarization matrix

**is shown in appendix A to be a product of these projection matrices: if there are**

*P**n*scatterings (

*n*−1 intermediate segments) the product is the back-to-back form,(3.1)For single scattering (

*n*=1), the matrices can be removed from the probability integration(3.2)For the double scattering contribution, the matrices

*M*_{1}cannot be removed since they depend on the segment direction being integrated (and the

*M*_{S}is stuck between them).(3.3)where the average 〈 〉

_{1}is over random scattering directions 1 on a hemisphere (upper if

*t*

_{2}>

*t*

_{1}and lower otherwise). The random azimuth part of the averaging is easy, and leads to matrix elements which are even polynomials of degree at most four in

*μ*

_{1}(the sine of the random elevation angle of direction 1).(3.4)(3.5)where the three superscripted matrices , , are constant, simply related to

*M*_{S}(3.13) in a manner specified in (3.10)–(3.12). Averaging over

*μ*

_{1}, which is uniformly distributed , and using the definition of

*E*

_{n}(2.3), yields directly the polarization matrix (3.3) for double scattering(3.6)The three integrals in (3.6) can be evaluated analytically (e.g. in the scalar context in van de Hulst (1948)), the first one being the same as (2.7). The other two are somewhat more involved and will not be presented since they are of higher order, and even for the best approximate analytical formula later (appendix C), it will suffice to replace by

*E*

_{3}(0)=1/2 and

*E*

_{5}(|

*t*

_{1}−

*t*

_{2}|) by

*E*

_{5}(0)=1/4, in which case the integrals are easy.

The expression for the triple-scattering polarization matrix has nine terms. If indices *i* and *j* run through the three values 1, 3 and 5, the integrals (which, however, cannot be evaluated analytically) are(3.7)and these are coefficients of nine matrices whose elements are again linear combinations of the elements of *M*_{S} given by (3.10)–(3.12) below. Thus, the result for the triple-scattering polarization matrix is(3.8)Likewise for the fourfold scattering, the expression is(3.9)and so on for higher numbers of scatterings, with the integrals defined as the obvious generalizations of (3.7), and the superscripted matrices given by iteration (one iteration for each superscript) as follows. If(3.10)(3.11)and(3.12)For instance, for which is required for the main result (4.5),(3.13)It is perhaps worthwhile to express in a coordinate independent way analogous to , namely − where ** S** is the Sun direction unit vector and

**is the vertical unit vector (and so**

*z**μ*

_{s}=

**.**

*z***). This completes the description of the exact Rayleigh scattering series.**

*S*## 4. Formulae from series truncation

As presented earlier, the contributions to the polarization matrix ** P** for the different numbers of scatterings all have the 3×3 form

*M**middle matrix*

**where**

*M***is the chosen line of sight projection matrix. It is more usual in describing polarization matrices to sacrifice the global neatness of the 3×3 form for a more economical 2×2 form. Re-expressed with respect to orthogonal axes of which the line of sight is one, the 3×3 polarization matrix**

*M***has only four non-zero elements (since the projections**

*P***on either side have cut out the others). The conventional second orthogonal axis is horizontal, parallel to lines of co-latitude (constant**

*M**θ*, increasing

*ϕ*and the third one parallel to lines of longitude (constant

*ϕ*, increasing

*θ*). Then the first row and column elements of the rotated polarization matrix are all zeros. Usingone obtains for the remaining elements(4.1)

**is the rotation matrix which takes the line of sight unit vector (cos**

*R**θ*cos

*ϕ*, cos

*θ*sin

*ϕ*, sin

*θ*) to (1, 0, 0) and the horizontal vector perpendicular to it to (0, 1, 0). The middle matrix always has four zero elements (

*xy*,

*yx*,

*yz*and

*zy*, the Sun being in the

*xz*plane, with

*ϕ*=0), and naming the remaining ones

*α*,

*β*,

*γ*and

*δ*(with the obliquity factor 1/sin

*θ*separated out), gives the 2×2 reduction of the matrix

**matrix in the form (1.1) stated in §1.**

*P*(4.2)As remarked there, all the dependence on azimuth *ϕ* is explicit, the *α*, *β*, *γ* and *δ* do not depend on *ϕ*. As well as applying to the exact (unknown) ** P**, this reduction can be applied to any of the

*n*scattering contributions to

**. The brightness is the trace of the 2×2 matrix (since it equals the trace of the 3×3 one). Another useful way of writing the brightness directly from the 3×3 form is brightness=Trace(**

*P*

*M**middle matrix*

**)=Trace**

*M**middle matrix*−

**.**

*O**middle matrix*.

**, where**

*O***is the observation direction unit vector (using ).**

*O*To obtain an analytical approximate formula for daylight brightness and partial polarization, it is necessary to truncate the exact Rayleigh scattering series. The crudest formula comes from single scattering; then including double scattering rectifies its main deficiencies, and finally including the leading terms of triple scattering yields the best analytical formula.

The single-scattering formula for the polarization matrix is already explicit (3.2), and its trace gives the single-scattering brightness of daylight, namely(4.3)(using . The order in *T* of the error in this approximation is the leftmost omitted term in table 1, namely *T*^{2} log *T*. The single-scattering expressions, (3.2) and (4.3), contain an infinity of terms of higher order than this, which one could discard at one's discretion, but the expressions are simple enough as they are. The brightness formula (4.3) is fairly accurate unless the elevation sines *μ* and *μ*_{s} are both small (figure 2).

However, single scattering is inadequate for describing the pattern of partial polarization in the sky because the pattern it produces is infinitely special (Berry *et al*. 2004), the polarization directions all being tangent to circles centred on the Sun direction. The Sun direction daylight is unpolarized or neutral. This arrangement is unstable to perturbation, which splits this neutral point into two (the pair being foci of an elliptical coordinates-like pattern instead of polar coordinates-like). In addition, single scattering is formally deficient, in that, just as in the scalar case, it is not the uniformly correct limit for small thickness *T* (figure 2). To rectify these deficiencies, it is necessary to include at least the leading order of double scattering. This is simply, from (3.6),(4.4)Together with the single scattering, this yields a (uniformly valid) formula for the brightness and polarization of daylight. If *E*_{1}(*T*) is replaced by its leading order, −log *T*, the formula remains uniformly valid but contains only elementary functions and constitutes, therefore, the main useable result for the polarization matrix of daylight.(4.5)The information is displayed graphically, for a chosen Sun elevation in figure 3. The order of the error in this approximation is the leftmost omitted term in table 1, namely *T*^{2}. To reduce the 3×3 matrix to the 2×2 form (1.1) and (4.2), one may first denote by (1/*μ*)*f*(*μ*) and (1/*μ*)*g*(*μ*) the prefactor functions from (4.5) so that . Then substituting *M*_{S} and from (3.13), (4.2) gives, recalling *μ*=sin *θ*,(4.6)If neither *μ*_{s} nor *μ* is small (more restrictive than the ‘not both’ small condition, earlier), *f* and *g* are approximately constant, independent of viewing direction ( and ), and therefore so are *α*, *β*, *γ* and *δ* (as mentioned at the end of §1). In the 3×3 form, they supply an approximation of (4.5)(4.7)where the matrix in the brackets is independent of viewing direction. This is the constant *α*, *β*, *γ* and *δ* circumstance described in §1. It has the unpleasant feature in the small *μ* limit (where (4.7) is not supposed to apply anyway) that the sky brightness becomes infinite—the blazing horizon problem. However, exactly the same pattern of directions as (4.7) arises without the divergence from a definite mean field-like physical model (albeit rather artificial) described in appendix E.

To visualize the polarization pattern of (4.7) (and Berry–Dennis–Lee), across the sky one can imagine looking, from a chosen direction in the sky, down that direction towards a small elliptical plate fixed on the ground with a certain tilt; sloping somewhat down its long axis. The perceived (i.e. projected) ellipse depends, of course on the chosen direction, and the polarization pattern on the sky is everywhere along the minor axis direction of the perceived ellipse. The ellipticity and tilt of the plate are such that the perceived shape is exactly circular when viewed from a neutral, unpolarized, direction in the sky. The ‘mirror reflection’ in the plate of this neutral direction is clearly a second neutral direction from which the plate again looks circular. Because the tilt slope of the plate is down its long axis, these two are on the same meridian great circle. (If the ‘reflected’ direction is under the ground, then its opposite direction should be taken instead.) In summary (with the logic correctly reversed), the input knowledge of the two neutral directions in the sky fixes the required ellipticity and tilt of the plate, and thence the polarization pattern of (4.7) (and Berry–Dennis–Lee) across the sky.

A formula for the neutral directions (figure 4) derives directly from the form (1.1) or (4.2), setting sin *ϕ*=0 (if the neutral points did not to lie on the Sun's meridian the solutions would be imaginary, not real). Solving the quadratic trigonometric equation (converting to double angles) gives the angles of elevation of the neutral points as(4.8)The four solutions are in two antipodal pairs, so exactly two are in the sky (with their opposites ‘underground’). In the constant *α*, *β*, *γ* and *δ* approximation to (1.1) or (4.2)), these solutions are explicit; otherwise, when *α*, *β*, *γ* and *δ* depend on *θ*, they are implicit, requiring numerical evaluation.

## Acknowledgments

I am grateful to M. R. Dennis and M. V. Berry for their comments. A few years ago each of us had made some preliminary scattering calculations for daylight, unaware of the general simplification and systematization possible, but a valuable foundation, nonetheless.

## Footnotes

- Received April 3, 2007.
- Accepted July 9, 2007.

- © 2007 The Royal Society