Radiative albedo from a linearly fibered half-space

A growing acceptance of fiber-reinforced composite materials imparts some relevance to exploring the effects which a predominantly linear scattering lattice may have upon interior radiative transport. Indeed, a central feature of electromagnetic wave propagation within such a lattice, if sufficiently dilute, is ray confinement to cones whose half-angles are set by that between lattice and the incident ray. When such propagation is subordinated to a viewpoint of an unpolarized intensity transport, one arrives at a somewhat simplified variant of the Boltzmann equation with spherical scattering demoted to its cylindrical counterpart. With a view to initiating a hopefully wider discussion of such phenomena, we follow through in detail the half-space albedo problem. This is done first along canonical lines that harness the Wiener-Hopf technique, and then once more in a discrete ordinates setting via flux decomposition along the eigenbasis of the underlying attenuation/scattering matrix. Good agreement is seen to prevail. We further suggest that the Case singular eigenfunction apparatus could likewise be evolved here in close analogy to its original, spherical scattering model. A cursory contact with related problems in the astrophysical literature suggests, in addition, that the basic physical fidelity of our scalar radiative transfer equation (RTE) remains open to improvement by passage to a (4×1) Stokes vector, (4×4) matricial setting.


Introduction
Fiber-reinforced composites 1 have, in recent decades, enjoyed a widespread penetration into the manufacture of durable structures, most noticeable among them being perhaps aircraft fuselage segments on a large scale, an arena dominated heretofore by aluminum as the material of choice. And, while the primary objectives of reduced weight and cost have indeed been achieved, subsidiary issues of somewhat lesser importance have nevertheless arisen. One among the latter is FRC behavior vis-à-vis radiative transport at μm wavelengths in the infrared.
Reinforcing fibers within some uniform host matrix are most frequently encountered at two opposite extremes of dispersal regularity/irregularity, to wit, short tendrils randomly oriented, or else extended threads aligned along a single direction. Evidently it is only the FRC samples belonging to the second, aligned category which can be expected to exert any significant orientational influence upon radiative transport.
With this background in mind, the present note seeks to offer a very modest first insight into FRC radiative transfer by solving the albedo problem for an otherwise uniform half-space matrix 2 randomly seeded with long (idealized, to be sure, as infinitely long) fibers in sufficiently dilute 3 distribution parallel to the material interface.
One should hasten to acknowledge that, even though the present effort was initially motivated by relatively utilitarian considerations, radiant transport across randomized swarms of elongated particles, typified by ice needles 4 in terrestrial clouds [1], and by partially oriented interstellar dust grains [2], had long ago attracted the attention of astrophysicists seeking to analyze light reflected/emitted from a variety of cosmic environments, planetary atmospheres a e-mail: jan.grzesik@hotmail.com 1 Abbreviated henceforth as FRC, singular or plural as the local context may dictate. 2 Admittedly, in what follows, we allow this matrix to default to empty space. This is done so as to avoid having interface refractive complications obscure the details of our evolving radiative transport machinery. Needless to say, this physical/mathematical lacuna pleads to be filled by more realistic work in the future. 3 By dilute we mean here an average fiber separation of at least several wavelengths. 4 A thorough discussion of the extent to which finite cylinder scattering can be approximated by that due to its infinite idealization can be found in [1]. prominent among them. Indeed, a standard text [3], in wide circulation since first published in 1957, is devoted to light scattering by the types of particulates that astrophysicists expect to encounter. The global theory of radiative transport, on the other hand, has been dominated by the two magisterial monographs [4,5], with both acknowledging their debt to the invariant embedding viewpoint first advocated by Ambartsumian. Planar slab and/or bona fide halfspace albedo problems, akin to that presently considered, naturally made their way into the astrophysical research literature during attempts to understand light reflection from these atmospheres in their geometrically idealized limit. Early samples of such work, which draws on the methods of both Chandrasekhar and Sobolev, phrased in terms of Stokes four-element vectors so as to keep electromagnetic polarization 5 in full view, can be found in [6][7][8], and have doubtless spawned their own research progeny.
It is characteristic of an electromagnetic field impinging upon a uniaxial scatterer as now described to propagate along directions confined to cones having the fiber direction as axis and opened to the angle between that axis and the direction of incidence 6 , angle π/2 − φ in fig. 1 7 . Since conical confinement of this sort is indifferent to the incident field polarization, a welcome invitation presents itself to sidestep the Maxwell electromagnetic apparatus in favor of the much simpler Boltzmann equation governing radiative transfer. Such transfer need clearly be tracked only in its projection upon the plane perpendicular to fiber direction, the full conical flux being gotten therefrom under a simple multiplication by sec φ as indicated in fig. 1. Plainly put, the only positional variables with which the Boltzmann transport equation need concern itself are the angle θ about fiber direction ( fig. 1) and the depth of penetration, here taken as coördinate y, into the fiber-laden half-space. 5 Our more primitive method, let it be admitted at the outset, altogether fails to discriminate among electromagnetic polarization states, which, together with rays/wave normals, form a mutually orthogonal triad. 6 This geometrically engaging attribute is, regrettably, passed over in silence by the standard electromagnetics texts (such as those authored by Stratton, Jackson, Panofsky and Phillips, and so on). Sheltered though it may thus be, its persistent presence in the literature is securely anchored around [9] and in research papers of the type [10][11][12][13] to which it lends support. One finds there a derivation which relies upon the peculiarities of an essentially two-dimensional field propagation, complete with reference to the asymptotic behavior of the Hankel functions which figure therein. By contrast, an ab initio derivation from first principles, altogether liberated from any debt to special function features, can be found in [14]. This minor caveat notwithstanding, reference [9] is widely regarded as having at least rivaled, if not altogether eclipsed [3]. The manner of derivation adopted by [14] confirms moreover that the conical scattering attribute prevails regardless of cylinder cross-sectional shape, and is controlled by the (vector-valued) two-dimensional Fourier transform of the self-consistent ohmic/polarization currents radiating across that cross section. The currents themselves must of course be ascertained on the basis of an integral-equation self-consistency argument, or else by enforcement of somewhat more pedestrian boundary-value continuity requirements. One may note with some astonishment that [3] restricts itself to the case of a purely perpendicular wave incidence when setting out the specific details in its treatment of cylindrical scattering. 7 Figure 1 is freely available from the internet under a Google search "cylinder scattering". Its origin can be traced to [15], and it reappears on p. 264 of [9].
Our objectives here are very limited indeed, barely embryonic if one may be permitted to say so. On the one hand, we shall not venture into the prevailing technology of FRC manufacture, while, on the other, our analysis will merely seek to provide a methodological preamble to more refined, more sophisticated computations which may yet appear in the future from the pens of other hands. Indeed, on the methodological side we shall restrict ourselves, as already stated, to a uniform, uniaxially fibered half-space, and, moreover, to a scattering function within Boltzmann's equation which is invariant against azimuthal angle θ around fiber direction (coördinate axis z as indicated in fig. 1). Our primitive goal will be the interface albedo under plane wave irradiation at azimuthal angle θ 0 (and its cosine μ 0 = cos θ 0 ), reckoned from the positive y-axis and arbitrarily chosen in the impact range −π/2 < θ 0 < π/2 (see footnote 8 ).
Within these strict confines we shall pursue our albedo target along two routes, one quasi-numerical, based upon the Boltzmann equation being discretized along angular nodes conforming to a Gaussian quadrature net, the other more traditionally analytic and relying upon a Wiener-Hopf complex plane machinery. Essentially perfect agreement will be demonstrated between the albedo outcomes at the end of these distinct routes, and brought to the fore will also be a standard incoming/outgoing symmetry that harkens back to the Ambartsumian/Chandrasekhar/Case machinery of invariant embedding [4,5] and [16,17]. All of this, it is hoped, may yet stimulate more refined, more realistic analyses by other investigators engaged in radiant transport across fibered media.
The material thus described is organized as follows. After an obligatory setup of albedo geometry and the underlying Boltzmann equation, pride of place is given to the Wiener-Hopf machinery. Sketched next are the very simple details of transport equation discretization along the nodes of a concatenated Gaussian quadrature lattice extended across the full, −π < θ < π range of scattering azimuth θ. There follow polar albedo plots across the outgoing range π/2 < |θ| < π, ϑ = ±||θ| − π|, with special emphasis upon a virtually perfect agreement between calculations having an entirely dissimilar provenance. Our note then concludes with one or two suggestions as to possible extensions of this work.

Albedo setup
Right-handed Cartesian axes, their origin in the half-space boundary, are placed so that horizontal axis z is aligned with fibers as in fig. 1, axis y points into the scattering medium, while axis x is regarded as the vertical. As already indicated, azimuthal angle θ, −π < θ < π, is reckoned from axis y, positive or negative when conveying a right-hand rotation around axis ∓z.
Powerless to change, we betray an atavistic fidelity to the nomenclature of neutron transport by writing now ψ(θ, y) for the radiative intensity per unit increment of scattering angle θ. The total macroscopic cross section (absorption plus scattering, σ a + σ s ; dimension of an inverse length) is written as σ, with both of its constituents assumed to be spatially constant as to both position y and angle θ. Normalization by σ yields the dimensionless single-scattering albedo ω = σ s /σ, 0 ≤ ω < 1 (see footnote 9 ), and an optical path length τ = σy which is y when measured in units of one mean free path. Playing a key rôle is a probability {P (θ | θ )/(2π)}dθ that radiant intensity ψ(θ , τ), arriving along direction θ , is scattered into angular increment dθ around θ. One evidently must insist upon the normalization π −π P (θ | θ )dθ = 2π, regardless of incidence angle θ . In what follows we simply default to the primitive, isotropic case P (θ | θ ) = 1 regardless of depth τ . In standard notation, μ = cos θ.
With these obligatory definitions duly disposed of, we may simply set down as our basic energy bookkeeping (Boltzmann) equation 10 in the x-y plane transverse to the mandated fiber direction. 8 An accompanying fiber/ray tilt angle −π/2 < φ0 < π/2, as indicated by fig. 1, should likewise be retained in the recesses of one's mind, but need not otherwise enter the details of Boltzmann equation calculations. 9 We shall expressly assume here that the fibers absorb to a greater or lesser extent, so that indeed ω < 1. Calculation of the macroscopic parameters σ a,s is a moderately involved electromagnetic task in its own right, one that we do not propose to undertake here. How such computations must proceed can be found in [3] and [9], and, when scattering cylinders are radially stratified, in [10][11][12][13] and the references cited therein. 10 It has been a credo of radiative transport theory, since at least the publication of [4], that radiant flux requires for its complete description a four-element Stokes parameter vector. Our present flux ψ, merely in the interest of methodological simplicity, defaults to merely the top element of this vector, and its submission to Boltzmann's equation (RTE: radiative transfer/transport equation, in contemporary parlance) (1) betrays a primitive, phenomenological attitude to radiative transport, an attitude since refined on the basis of Maxwell's electromagnetic equations and reported, in chronological order, by [18,19], and [20]. Of these, [18] is the most succinct, [19] most detailed, whereas [20] threads a wonderful panorama of scientific history into its critique of the more antiquated phenomenology. Stated in so many words, our radiant flux ψ exempts itself from any obligation to account for electromagnetic field polarization around propagation cone generators. At the same time one must acknowledge that radiative polarimetry, which figures prominently in remote sensing of surface irregularities and in geologic prospecting, does make ample use of the matricial Stokes apparatus. Reference [18] provides the background theory for this terrestrial technology.

Wiener-Hopf solution for ψ d (θ, τ )
We displace attention from ψ d (θ, τ ) per se to its angular aggregate It turns out that nothing is lost thereby since, on the basis of (3), we can further write throughout the entire retrograde flight range π/2 < |θ| < π, wherein −1 < μ < 0. On setting τ = 0 we see that knowledge about ρ d (τ ) provides a direct stepping stone to the desired albedo. A path to solution for both ψ d (θ, τ ) and ρ d (τ ) is cleared by first subjecting these quantities to Laplace transformation, denoted by a circumflex, viz.,ψ with transform variable s initially required to have its real part positive subject to restrictions identified in the course of the ensuing argument. Indeed, we see at once, with τ set equal to zero, that (4) may now be simply rephrased as Laplace transformation of (3) next gives which can further be integrated over angle so as to yield The second integral on the right incorporates the null reëntrant flux condition at the half-space boundary τ = 0. Since it requires knowledge about the albedo, it is a priori unknown. Such information lapse notwithstanding, it is evident by inspection that that second integral is analytic in a left half-plane having s < 1.
We embark next on a standard Wiener-Hopf journey wherein the ingredients of (9) are suitably rearranged as to their left/right half-planes of analyticity, the end result being a recognition of (9) as a bridge uniting ostensibly disparate, left/right analytic functions into one global entire function bounded at infinity and therefore, by virtue of Liouville's theorem, a mere constant. Self-evident asymptotic estimates require that constant to be null, at which point all solution details fall into their rightful place.
We begin with the definite integral which is clearly analytic when | s| > 0, regardless of s, but, strictly speaking, is undefined without further qualification when instead s = 0 whereas | s| > 1. It can be rewritten in the usual way as a contour integral around the unit disk |ζ| = |e iθ | = 1 and thus evaluated by adding ±2πi times the sum of its residues at one, two, or none of the simple poles which its interior/exterior may contain.
In keeping with a precautionary remark above, we provide √ s + 1 with branch cuts radiating outward from ±1 to ±∞. Accordingly, i √ s 2 − 1 is real and negative when s = 0 while −1 < s < 1. Since also ζ + ζ − = 1, it follows that only one such pole is actually enclosed, and, because |ζ + | → ∞ as s → 0, whereas ζ − → 0, that enclosed pole lies at s = ζ ∓ in disk interior/exterior 11 . Hence Returning now to the main theme, and with the information from (13) explicitly displayed, we execute our first rearrangement by adding and subtracting μ 0 /(1 + sμ 0 ) on the right in (9) so as to arrive at √ wherein we further abbreviate by setting As is easily ascertained, the numerator on the left in (14) has two simple zeros at 12 whose presence dampens the prospect of its imminent, intended use as the argument of a logarithm. We correct for this by first dividing out (s − s + )(s − s − ) = s 2 − s 2 ± and then compensating for its adverse effect as s → ∞ by multiplying with s 2 − 1 (see footnote 13 ). And so we mold (14) into 11 We note from (12) that ζ± coincide at ∓1 when s = ±1. In general, But then, as s descends from s = 1 to s = 0 along its axis of reals, ζ− rises from ζ− = −1 to ζ− = 0, while simultaneously ζ+ migrates downward from ζ+ = −1 to ζ+ = −∞. Upward s movement from s = −1 to s = 0, with s = 0 throughout, drags ζ− downward from ζ− = 1 to ζ− = 0, and pushes ζ+ upward from ζ+ = 1 to ζ+ = ∞. Evaluation (13) is thus confirmed along the entire interval −1 < s < 1, s = 0, whence analytic continuation suffices to validate it also for complex arguments s. 12 In an appendix (sect. 3.2.1 below) we show that s± are the propagation constants associated with source-free solutions of (1) (restated there as (34)). In particular, s −, conveying an unbounded spatial growth exp(−s−τ ) for τ → +∞, has relevance for the kindred, right half-space Milne problem, wherein it represents some unspecified, deeply buried intensity source of unlimited strength. 13 Beyond this restricted purpose, multiplier s 2 − 1 evidently serves to suppress denominator infinities at branch points s = ±1.
One immediate, but only partial step toward rearranging (18) so as to achieve a left-hand side analytic in a right half-plane, and a right-hand side analytic in a left half-plane is to set The real difficulty lies in achieving a similar split for function κ(s), a task accomplished via the known recipe of passing first, seemingly redundantly, to the exponential of its logarithm, and then representing that latter as a Cauchy contour integral with up/down legs of infinite extent parallel to the imaginary axis of s. And so, if we set, for some real positive β, 0 < β < s + (see footnote 14 ), then we obtain at once the decomposition valid at least in the vertical strip −β < s < β. Moreover, functions κ ± (s) are readily seen to be analytic respectively across the overlapping left/right half-planes s < β or else s > −β, in that order, and indeed to vanish on approach to infinity, |s| → ∞. Equality (20), when rewritten now as represents thus yet another step in segregating quantities that are equal within a band of overlap but, otherwise, assert their analyticity throughout dissimilar half-planes.
Only one further adjustment remains to make that segregation complete, and that is subtraction from g(s) of its simple pole at s = −1/μ 0 . The subtraction must of course occur on both sides of (23), which thus yields as a bona fide bounded entire function, null at infinity, and thus null everywhere on the strength of Liouville's theorem. But that means that our solution is in hand, inasmuch as (15) gives, first and then the desired albedo follows from (7) as with μ < 0, μ 0 > 0, so that arguments of κ ± properly lie in their respective half-planes of analyticity. Somewhat greater symmetry in μ 0 and μ is attained by noting from (21) that since, by virtue of (19), κ(s) per se is symmetric, κ(s) = κ(−s), under argument reflection through the origin. On replacing μ with −μ 0 (27) also gives κ − (1/μ 0 ) = 1/κ + (−1/μ 0 ).

Fig. 2. Contour deformation for calculating κ−(s).
Hence Furthermore, if we set 15 then it follows at once that a symmetry under interchange of incoming/outgoing intensity directions at the half-space boundary, reciprocity traditionally associated with invariant embedding [4,5] and [16,17]. With a view to the first line in (29) it remains now only to deform the vertical contour whereby κ − (s) is defined in (21) into a form more convenient for numerical implementation.

Contour deformation
Of the two options among κ ± which eq. (29) provides, we choose the first and proceed to deform its contour around the branch cut extending from s = −1 to s = −∞, as shown schematically on fig. 2. Potential contributions from both the semi-circle at infinity and from the full, but vanishingly small circle around the branch point at s = −1 are dismissed in the usual way, and what remains is the difference of integrating immediately above and below the branch cut. We thus find in short order that the second line following under the natural variable substitution ω ζ 2 − 1 = tan ϑ.
In its latter guise the computation of κ − (s) becomes amenable to numerical integration. Entirely similar reasoning, but with the contour wrapped now around the complementary branch cut extending from s = 1 to s = ∞, yields And then, if we recall that (32) holds for s > −β (in fact, with even less restriction for s > −1), whereas (33) requires instead that s < β (which can be relaxed to the right even so far as s < 1), it becomes evident that eqs. (32)-(33) abide by the symmetry embodied in eqs. (27)- (28) and, indeed, can be gotten on the basis of that symmetry, one from the other, without any additional calculations 16 .

Free modes
In footnote 12 , accompanying eq. (17), we had already claimed that quantities s ± = ± √ 1 − ω 2 are associated with source-free transport modes. Here we verify this assertion in two ways, the second of which will anticipate a Milne-type integral equation for radiant density ρ d (τ ). Hence if we begin with (eq. (1) restated) and seek a solution having an exponential behavior ψ(θ, τ ) = e −λτ η(θ), with λ strictly real, we find in the usual way that (34) requires once e −λτ π −π η(θ)dθ has been cancelled from both sides. Reference to eqs. (10) and (13) and to the intervening discussion shows next that, in order to produce on the right in (35) a result that is manifestly real, we must trap λ between −1 and 1, obtaining or else λ = s ± as claimed.
The end goal of (17) can likewise be reached by converting (34) to an integral equation for ρ(τ ), at least when the scattering medium fills all of space, −∞ < τ < ∞. Such conversion is easily attained by use of integrating factors, which give when μ > 0, and 16 There is a considerable amount of hidden latitude in the choice of vertical contours in (21). Within obvious limits they need not be placed at strictly the same distance β, left and right, from the imaginary axis, and, again within limits, they could be allowed to undulate to a certain extent while maintaining overall a steady upward progress. But these are totally dispensable frivolities, having impact neither upon the transport phenomenology, nor upon the mathematics. Furthermore, each of κ ±(s) as it now stands defines two analytic functions in adjacent, open half-planes filling up the entire two dimensional s-space, save for a boundary along its vertical integration path. But once again this is a recherché aspect without physical relevance to our problem, and thus in no need of being belabored.
if instead μ < 0, the distinction embodying a self-evident physical causality. Integrating over all angles we thus obtain Setting ρ(τ ) = α e −λτ , with α some irrelevant but positive constant, and real λ less than 1 in magnitude, −1 < λ < 1, so as to assure quadrature convergence even when μ = 1, we find next, with interchange of integration order taken for and And then, with (40)-(41) brought to bear upon (39) we find which is nothing other than (36) once more 18 .

Milne-type integral equation
Equations (37)-(39) suggest that an integral equation can likewise be contrived for the somewhat more physically relevant diffuse density ρ d (τ ), associated with eq. (3) and having only the half-line τ ≥ 0 for its support. Analogues to (37)-(38) emerge now as when μ > 0, and when instead μ < 0. Integration over all angles then gives Of the three angular quadratures on the right in (45), only 17 That is to say, Fubini's theorem is assumed to be in force. 18 The middle line of (42) follows from a routine appeal to a half-angle trigonometric formula, whereas its last utilizes the evaluation of (10) as provided by (13). In this connection one readily verifies that the requisite inequality −1 < λ 2 /(2 − λ 2 ) < 1 is a consequence of having −1 < λ < 1.
submits to a closed-form evaluation 19 , the other two 20 being considerably more recondite. In particular, the kernel of integral equation (45) mimics its exponential integral analogue as normally encountered in connection with the Milne problem in a bona fide three-dimensional, isotropic scattering context [21]. Neither one of the kernels (47)-(48) seems to admit evaluation in closed form. This impediment notwithstanding, a Wiener-Hopf attack can be successfully mounted in the presence of kernel (48), and so presumably it could likewise be pursued with that of (47), essentially scrolling the analysis somewhat in reverse. But it is best at this point to ignore the siren call which beckons from deep within this detour.

Discrete ordinates
The Wiener-Hopf apparatus pivots around the global attribute ρ d (τ ), which no longer retains detailed flux memory ψ d (θ, τ ) along individual directions θ. Such details are restored, if only approximately, by discretizing the angular quadrature and then insisting that eq. (3) remain valid at the angular nodes 21 .
We thus imagine the discretized angular index l to run from 1 to some chosen N (see footnote 22 ) and organize the various discretized values ψ d (θ l , τ) into a column vector (49) 19 Result (46), confirmed by Mathematica, is obtained by emulating (13), but now in the form of an unlinked contour integral extended in the complex η-plane from η = −i to η = i, followed by appeal to a partial fraction decomposition of its integrand. Such decomposition exposes to view derivatives of logarithms whose ultimate remnant appears on the right in (46). Roots η± from (12) are now modified to read They both lie on the unit semicircle to the left, allowing the open contour integral above unfettered rein to advance from η = −i to η = i on the right. 20 About the third integral we can at least confidently assert, by virtue of its being symmetric in μ and μ0, that it is positive, so that it poses no danger whatsoever of pushing ρ d (τ ) into negative, physically forbidden territory. Its integrand, moreover, is clearly free from singularity at μ = μ 0 on the strength of L'Hôpital's rule. Of course, both first and second angular quadratures are positive by inspection. 21 A preferred, easily implemented discretization is obtained by partitioning the full, −π < θ < π angular range into any desired number of contiguous subintervals and, within each such, placing down a Gauss-Legendre (GL) net of some moderately high order. If 2a is the common subinterval length, then the canonical GL weights w appropriate to the interval (−1, 1) are correspondingly scaled down, w → aw. The global quadrature itself, evidently, is obtained as a simple sum over subinterval quadratures. Such concatenated GL discretization is easily deployed and has the merit of avoiding end-point bunching incident to use of just one master GL cycle of very high order (i.e., having many discretized nodes θ close to and above/below ∓π). Some interior bunching is of course inevitable, but in general there emerges a much more uniform composite GL net. 22 It is of considerable bookkeeping convenience here to have N divisible by four, with the second and third quartiles alluding to the forward moving intensity, while the first and fourth, when examined at the boundary τ = 0, underlie the angular albedo distribution.
Also required for the source term in (3) is the vector and an N × N matrix which accounts for the competition between absorption loss and scattering gain 23 .
It is a matter of empirical evidence that, in the Fortran codes which have been written around these ideas, all eigenvalues λ k of S, with 1 ≤ k ≤ N and Λ k being the corresponding eigenvector, have invariably turned out to be both real and distinct. We are admittedly unable to prove now that such a circumstance will prevail under all admissible parameter entries, but, should its truth be provisionally accepted as universal, then we are guaranteed by a well-known theorem that the eigenvector set {Λ k } N k=1 provides a basis for the N -dimensional linear vector space before us. Incidentally, the eigenvalues {λ k } N k=1 themselves have always been found to split evenly into positive and negative values, all negative, their complement {λ k } N k=N/2+1 all positive, symmetrically deployed in ± couplets, λ k = −λ N +1−k for 1 ≤ k ≤ N/2. The eigenvalue/eigenbasis framework thus revealed dominates all developments about to ensue, and, in order to impart a concrete symbolism to that framework, we set We thus duly arrive at as the discretized counterpart of (3), with real numbers {α k } N k=1 being the expansion coefficients of source vector (50) in our newly discovered eigenbasis, obtained by standard methods of linear algebra. Integrating upward from τ = 0 gives 24 If we next expand Ψ d (0) along the eigenbasis, 23 In (51), δm,n is the familiar Kronecker delta, equal to 1 when m = n and 0 otherwise. 24 As was similarly found in connection with the last term on the right in eq. (45), neither does there intrude into (55) any need to fear encroachment by some one λ k > 0 upon 1/μ0. Indeed, in both cases l'Hôpital's rule saves the day. But, should such l'Hôpital limit extraction ever be required, we would still be assured of an ultimate asymptotic decay, proceeding now as τ e −τ/μ 0 in that particular mode. with interface amplitudes {γ k } N k=1 , we find that Diverging exponentials with λ k < 0 are banished at once by setting, for 1 ≤ k ≤ N/2, with a complete impunity now as to any menace of having to cope with vanishing denominators, while the remaining amplitudes γ k at indices N/2 + 1 ≤ k ≤ N follow by suppressing any reëntrant flux, required now for all l in the range N/4 + 1 ≤ l ≤ 3N/4. Admittedly the N/2 × N/2 linear system (59) is still in need of its own solution, but that is a relatively standard matter vis-à-vis the wide availability of matrix inversion routines 25 .

Numerical examples
Since we have no wish to inundate the prospective reader with an avalanche of numerical data, we provide just two albedo examples based on the developments that have been sketched, showing a most welcome agreement between Wiener-Hopf versus discrete ordinates outcomes, an agreement which surely bolsters confidence in both. The above analysis has been implemented in Fortran code with both Wiener-Hopf and discrete ordinates viewpoints running in parallel 26 . Utility subroutines for eigenvalue/eigenvector computation, linear algebraic system solution, and numerical integration, the latter on behalf of the Wiener-Hopf albedo result (26) et seq., have all been drawn from the IMSLIB 27 libraries.
The sample computations of figs. 3 and 4 have been performed on an angular grid of 40 concatenated Gaussian level-10 quadratures, yielding 100 points per quadrant. As a way to ease visualizing a performance comparison, the discrete ordinates outcomes were then simply winnowed down by a factor of 5 so as to yield just 20 points per quadrant. It should come as no surprise that substantial absorption, as conveyed by a single scattering albedo ω = 0.60 in fig. 3, depresses the emergent flux across the board as compared to that from an essentially dissipation-free medium (ω = 0.99) in fig. 4 28 . The cumulative albedos perpendicular to the boundary and normalized by the similarly perpendicular in-

Comments
It is clear that the work now on view builds directly upon the standard formulation, as found in, say, [17] and [21], wherein the scattering phenomenon is fully three-dimensional, and the angular analysis is equipped with a weight appropriate to spherical quadrature. But, of course, as we all know, the devil is ever in the details, so that perhaps these efforts may yet justify their existence as a model for further work by others. The next logical step in the adaptation would be to emulate a flux development along the lines of Case-type singular eigenfunctions. In fact, a first step in this direction had already been taken quite a while ago, and was hastily assembled for presentation at the 1997 IEEE AP-S/URSI Symposium in Montréal, Canada [22] 29 . But this foray was only a hint of what can and should be done, and doubtless remains to this day in need of a patch up. Perhaps such a patch up may yet be undertaken by the undersigned, and/or at the hands of some other transport enthusiast(s).
And lastly, even though we seem to have shirked excessive contact with the existing transport literature, be its focus terrestrial or astrophysical, we most certainly have no wish to impart the patently erroneous impression that our pages furnish the one and only, the solitary foray into the realm of radiative transport across granular, inhomogeneous media. On the contrary! Our only claim is to attempt a contribution toward the understanding of radiant transport wherein a prevailing background symmetry, presently of the strictly linear sort, is sufficiently compelling to impress its stamp upon the propagating field structure. Antecedents in this latter direction, even if their overt concern may appear to be strictly electromagnetic, can be traced from [10][11][12][13], and the present work is a direct echo thereof, transplanted into its native, transport setting.
By contrast, most available treatments of radiant transport across granular media confine themselves, by virtue of both available analytic tools and the considerable domains of practical applicability, exemplified by thermal insulation blankets and polluted atmospheres, to particle and tendril swarms, the latter naturally randomized as to individual orientation. One such study, authored by Tien, a master of radiant heat transport, and Tong, one of his many students, can be found in [23], while a later complement, due to Lee, yet another former student, previously cited under [11][12][13], is reported as [24].
An early treatment of radiant transport in the presence of particulate dispersal is exemplified by the Purdue University dissertation of Love, Jr. [25], wherein one finds a vast amount of numerical work implementing a variety of scattering series solutions. The scope of this numerical effort surely deserves to be judged as monumental when viewed against the relatively impoverished computing resources available at the time of its creation. By contrast, in the purely theoretical work reported in [26], Whitaker provides for the temperature a diffusion equation whose highly modified heat flux seeks to describe porosity effects.
And, while radiant transport theory, seen from an engineering perspective, is well conveyed in the standard pillars [27,28], it is only the latter which devotes much attention to medium inhomogeneities.
Thanks are due to an anonymous referee for broadening our viewpoint so as to include reference to kindred astrophysical research. Even though the format at hand allows for little more than a passing mention of it, we, and hopefully also most potentials readers, will be sufficiently intrigued to explore the Stokes vector/matricial RTE refinements which it suggests. The newly encountered conical aspect of field propagation will doubtless render such exploration a non-trivial, intriguing task, a substantial reworking of the spherical coördinate material found under sources [4], and especially [18][19][20].
Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.