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 radiant 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 photon 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.


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 which belong 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.
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. 4 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 Figure 1. 5 Plainly put, the only positional variables with which the Boltzmann transport equation need concern itself are the angle θ about fiber direction ( Figure 1) and the depth of penetration, here taken as coördinate y, into the fiber laden half space.
Our objectives here are very limited indeed, barely embryonic if one 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 lay claim merely to 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 Figure 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. 6 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 quad-rature 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 Ambarzumian/Chandrasekhar/Case machinery of invariant embedding [7][8][9][10]. All of this, it is hoped, may yet stimulate more refined, more realistic analyses by other investigators engaged in radiant transport across fibered media.

Figure 1. Cylinder scattering cone C under oblique ray impact
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 axis z is aligned with fibers as in Figure 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 photon flux 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, 7 and an optical path length τ = σy gotten when y is measured in units of one mean free path. Playing a key rôle is a probability P (θ|θ )/(2π)dθ that a photon 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. In standard notation, µ = cos θ.
With these obligatory definitions duly disposed of, we may simply set down as our basic photon bookkeeping (Boltzmann) equation in the x − y plane transverse to the mandated fiber direction. Flux ψ(θ, τ ) is further decomposed as the sum of a Dirac delta function sheet e −τ /µ 0 δ(θ − θ 0 ) incident at angle θ 0 and properly attenuated thereafter, and a diffuse, multiply scattered component Decomposition (2) converts the homogeneous equation (1) into an inhomogeneous counterpart whose source (ω/2π)e −τ /µ 0 is now entirely overt. And then, as a boundary condition on (3) we further require that there be no reëtrant diffuse flux at half-space boundary τ = 0, viz., ψ d (θ, 0) = 0 when |θ| ≤ π/2. By contrast, it is the structure of the escaping, diffuse flux ψ d (θ, 0) for π/2 < |θ| < π that we identify with the half-space albedo.
3 Wiener-Hopf solution for ψ d (θ, τ ) We displace attention from ψ d (θ, τ ) per se to its angular aggregate which, following division by the speed of light c, simply measures the diffuse photon density. 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., and with transform variable s initially required to have its real part positive. Indeed, we see at once, with τ set equal to zero, theat (4) may now be simply rephrased as Laplace transformation of (3) accordingly 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| < 1, which we now assume, 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 residues at one, two, or none of its simple poles at In keeping with the precautionary remarks above, we provide ζ 2 − 1 = √ s − 1 √ s + 1 with branch cuts radiating outward from ±1 to ±∞. In particular, i √ s 2 − 1 is real and negative when s = 0 whereas −1 < s < 1. Since also ζ + ζ − = 1, it follows that only the one pole at ζ − is actually enclosed. And so Returning now to the main theme, and with the information from (13) explicitly displayed, we execute our first rearrangement by adding and subtracting 1/(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 8 whose presence dampens the prospect of its imminent 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. 9 And so we mold (14) into with 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 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 + , 10 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, and indeed to vanish on approach to infinity, |s| → ∞, s < β or else s > −β, in that order. 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 the simple pole at s = −1/µ 0 . The subtraction must of course occur on both sides of (20), 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 Somewhat greater symmetry in µ 0 and µ is attained by noting from (21) that since, in the present instance, µ is negative and, by virtue of (19), κ(s) per se is symmetric, κ(s) = κ(−s), under argument reflection through the origin. On replacing µ with −µ 0 (27) also gives Hence = ω 2π Furthermore, if we set 11 then it follows at once that a symmetry under interchange of incoming/outgoing photon directions at the half-space boundary, traditionally associated with invariant embedding [7][8][9][10]. 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 Figure 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. 12 We thus find in short order that the second line following under the natural variable substitution 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 > −β whereas (33) requires instead that s < β, 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. 13

Free modes
In Eq. (17) and its accompanying Footnote 8 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 photon density ρ d (τ ). Hence if we begin with and seek a solution having an exponential behavior ψ(θ, τ ) = e −λτ η(θ), with λ strictly real, we find in the usual way that (34) requires 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 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, we find next, with interchange of integration order taken for granted, 14 and And then, with (40)-(41) brought to bear upon (39) we find which is nothing other than (36) once more. 15

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 µ < 0. Integration over all angles then gives Of the three angular quadratures on the right in (45), only submits to a closed-form evaluation, 16 the other two 17 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 [12]. 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 photon flight 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. 18 We thus imagine the discretized angular index k to run from 1 to some chosen N 19 and organize the various discretized values ψ d (θ k , τ ) into a column vector . . .
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. 20 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 symmetrically deployed around the origin, {λ k } N/2 k=1 all negative, their complement {λ k } N k=N/2+1 all positive. 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 {m 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 21 If we next expand Ψ d (0) along the eigenbasis, wherein the various amplitudes {b k } N k=1 allude to "boundary," we find that Diverging exponentials with λ k < 0 are banished at once by setting, with complete impunity as to any menace of having to cope with vanishing denominators, for 1 ≤ k ≤ N/2, while the remaining amplitudes b k at indices N/2 + 1 ≤ k ≤ N follow by suppressing any reëntrant flux, required now ∀ 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. 22

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. 23 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 libraries.
The sample computations of Figures 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 point per quadrant. It should come as no surprise that substantial absorption, as conveyed by a single scattering albedo ω = 0.60 in Figure 3, depresses the the emergent flux across the board as compared to that from an essentially dissipation-free medium (ω = 0.99) in Figure 4. 24 The cu-mulative albedos, perpendicular to the boundary and normalized by the similarly perpendicular influx, when integrated across the indicated, 200−point angular grid, ring in at 0.29253 . . . in the strong absorption case (Figure 3), and 1.14332 . . . . . . when instead a pure scattering regime predominates. We are at a loss to provide any physically credible motivation for the horizon bulge in Figure 4, but of course would welcome all plausible suggestions.

Comments
It is clear that the work now on view builds directly upon the standard formulation, as found in, say, [10] and [12], 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 [13]. 25 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 much contact with the existing transport literature, 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 [4] and [11], 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 C. L. Tien, a master of radiant heat transport, and T. W. Tong, one of his many students, can be found in [14], while a later complement, due to S. C. Lee, yet another former student, is reported as [15].
An early treatment of radiant transport in the presence of particulate dispersal is exemplified by the Purdue University dissertation of Tom J. Love, Jr. [16], 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 [17], Stephen 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 [18]- [19], it is only the latter which devotes much attention to medium inhomogeneities. 19 It is of considerable bookkeeping convenience here to have N divisible by four, with the second and third quartiles alluding to forward photon flux, while the first and fourth, when examined at the boundary τ = 0, measure the angular albedo distribution. 20 In (51), δ kl is the familiar Kronecker delta, equal to 1 when k = l and 0 otherwise. 21 As was similarly found in connection with Eqs. (43) and (45), neither does there intrude into (54) any need to fear encroachment by some one λ k > 0 upon µ 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 .
22 This fresh inversion prospect should not provoke consternation. Indeed, it is understood that by this point we had already confronted an even more daunting, N × N matrix inversion on behalf of the m k amplitudes first declared in Eq. (54). 23 In point of fact our code, at least insofar as its discrete ordinates part is concerned, aims at somewhat greater generality wherein the fibered, scattering medium has finite thickness and is sandwiched, fore and aft, by uniform, nonscattering dielectric blankets. All such potential generalizations, however, are simply bypassed, through suitable parameter choice, when undertaking the calculations now reported in Figures  3 and 4. One may note in passing that computer run times underlying such calculations are, to all intents and purposes, entirely inconsequential, even on a personal laptop device of modest power. 24 In both figures the incidence angle θ 0 has been set at 45 o . But, even though we do not dwell on this aspect here, the albedo patterns are only weakly dependent upon µ 0 = cos θ 0 . 25 At first blush an IEEE/URSI conference may seem to be an unlikely, even incongruous venue for transport work of this sort. The level of incompatibility was lessened, however, by slightly stronger allusions to the obligatory electromagnetic background, which is barely mentioned in these pages. The present Figure  2 has, incidentally, been drawn from this presentation.