Radiative transfer in half spaces of arbitrary dimension

TL;DR

This paper generalizes classical radiative transfer problems to arbitrary Euclidean dimensions, deriving solutions and properties that reveal universal behaviors and dimension-dependent differences in scattering and diffusion modes.

Contribution

It introduces a unified framework for radiative transfer in any dimension, extending classical 3D solutions to higher and lower dimensions with closed-form expressions.

Findings

Universal properties invariant to dimension d

Discrete diffusion mode not universal for d > 3 in absorbing media

Unexpected correspondences between different dimensions and scattering types

Abstract

We solve the classic albedo and Milne problems of plane-parallel illumination of an isotropically-scattering half-space when generalized to a Euclidean domain $\mathbb{R}^d$ for arbitrary $d \ge 1$. A continuous family of pseudo-problems and related $H$ functions arises and includes the classical 3D solutions, as well as 2D "Flatland" and rod-model solutions, as special cases. The Case-style eigenmode method is applied to the general problem and the internal scalar densities, emerging distributions, and their respective moments are expressed in closed-form. Universal properties invariant to dimension $d$ are highlighted and we find that a discrete diffusion mode is not universal for $d > 3$ in absorbing media. We also find unexpected correspondences between differing dimensions and between anisotropic 3D scattering and isotropic scattering in high dimension.

Tables

Table 1. Table 1: H 𝐻 H function moments of various orders α i subscript 𝛼 𝑖 \alpha_{i} and various dimensions d 𝑑 d , for c = 0.99 𝑐 0.99 c=0.99 .

	${\alpha }_1$	${\alpha }_2$	${\alpha }_3$	${\alpha }_4$	${\alpha }_5$	${\alpha }_6$	${\alpha }_7$	${\alpha }_8$
d=2	1.26655	1.02971	0.889917	0.794985	0.725119	0.670931	0.627316	0.591231
d=3	1.02718	0.721955	0.557658	0.45458	0.383773	0.332099	0.292712	0.261689
d=4	0.886147	0.556488	0.392232	0.295955	0.233692	0.190647	0.159405	0.135873
d=5	0.790625	0.452909	0.295462	0.208513	0.15522	0.120126	0.0957642	0.078151
d=6	0.720483	0.381903	0.232999	0.155012	0.109429	0.0809076	0.0616274	0.0487154
d=7	0.666181	0.330172	0.189885	0.119836	0.0806849	0.057	0.041798	0.03158

Table 2. Table 2: (Milne) extrapolation distance z 0 ​ ( d ) subscript 𝑧 0 𝑑 z_{0}(d) with ( c = 1 𝑐 1 c=1 ) in dimension d 𝑑 d .

$d$	$z_0\left(d\right)$
$1$	$1.0000000000$
$2$	$0.8183098862$
$3$	$0.7104460896$
$4$	$0.6366197724$
$5$	$0.5819457611$
$6$	$0.5393348406$
$7$	$0.5049086837$
$8$	$0.4763395251$

Table 3. Table 3: Two-media extrapolation distances z 0 ​ ( c 1 , c 2 , 2 ) subscript 𝑧 0 subscript 𝑐 1 subscript 𝑐 2 2 z_{0}(c_{1},c_{2},2) for 2D.

		$c_2$
$c_1$	0.6	0.8	0.9	0.95
0.8	1.5882436
0.9	1.3584815	2.1003544
0.99	1.2131854	1.6762401	2.3696245
1.01	1.1858938	1.6127114	2.2060217	5.5553
1.1	1.0791665	1.3925977	1.748799	2.7656641
1.2	0.9837306	1.222884	1.4652013	2.0342982
1.4	0.83958884	0.99755241	1.1399741	1.4286662

Table 4. Table 4: Two-media extrapolation distances z 0 ​ ( c 1 , c 2 , 3 ) subscript 𝑧 0 subscript 𝑐 1 subscript 𝑐 2 3 z_{0}(c_{1},c_{2},3) for 3D.

		$c_2$
$c_1$	0.6	0.8	0.9	0.95
0.8	1.4127134
0.9	1.175452	1.7835171
0.99	1.0265018	1.3985129	1.9595943
1.01	0.9987603	1.3402434	1.8181755	4.5360605
1.1	0.89132113	1.1377555	1.4192792	2.2264309
1.2	0.79703336	0.98183071	1.1697516	1.6118864
1.4	0.65879753	0.77729562	0.88422937	1.1004687

Table 5. Table 5: Two-media extrapolation distances z 0 ​ ( c 1 , c 2 , 4 ) subscript 𝑧 0 subscript 𝑐 1 subscript 𝑐 2 4 z_{0}(c_{1},c_{2},4) for 4D.

		$c_2$
$c_1$	0.6	0.8	0.9	0.95
0.8	1.3093552
0.9	1.0594888	1.5882436
0.99	0.90809624	1.2274169	1.7108654
1.01	0.8804248	1.1727963	1.5834373	3.9282009
1.1	0.77489064	0.9837306	1.222884	1.9098907
1.2	0.68446724	0.83958884	0.99755241	1.3693007
1.4	0.55567088	0.65367435	0.74209513	0.92048675

Equations

$${\Omega }_d=\frac{d{\pi }^{d/2}}{\Gamma \left(\frac{d}{2}+1\right)},$$

$$dE=I(\mathbf{x},\omega )\cos \theta d\sigma d\omega$$

$${\int }_{S^{d-1}}I(\mathbf{x},\omega )d\omega ={\Omega }_d$$

$$P({\omega }_i\to {\omega }_o)=1/{\Omega }_d.$$

$$dE=I(x,\mu )\mu G(\mu )d\mu d\sigma .$$

$$\frac{1}{2}{\int }_{-1}^{1}G(\mu )d\mu =1,$$

$$G(\mu )=\frac{2(1-{\mu }^2{)}^{\frac{d-3}{2}}\Gamma (\frac{d}{2})}{\sqrt{\pi }\Gamma \left(\frac{d-1}{2}\right)},d>1.$$

$$\frac{2}{\pi \sqrt{1-{\mu }^2}};1;\frac{4\sqrt{1-{\mu }^2}}{\pi };\frac{3}{2}(1-{\mu }^2);\frac{16(1-{\mu }^2{)}^{3/2}}{3\pi };\frac{15}{8}{\left(1-{\mu }^2\right)}^2.$$

$$\underset{d\to 1+}{\lim }G(\mu )=\delta (\mu -1)+\delta (\mu +1).$$

$$\frac{\sqrt{\pi }\Gamma \left(\frac{d+1}{2}\right)}{\Gamma \left(\frac{d}{2}\right)}{\int }_0^1\mu G(\mu )d\mu =1.$$

$$I(0,\mu )=\frac{\sqrt{\pi }\Gamma \left(\frac{d+1}{2}\right)}{\Gamma \left(\frac{d}{2}\right)},-1\le \mu \le 0$$

$${\pi }_d\equiv \frac{{\Omega }_d}{2}{\int }_0^1\mu G(\mu )d\mu =\frac{{\pi }^{\frac{d-1}{2}}}{\Gamma \left(\frac{d+1}{2}\right)},$$

$${\pi }_1=1,{\pi }_2=2,{\pi }_3=\pi ,{\pi }_4=\frac{4\pi }{3},...$$

$$f_L({\omega }_i,{\omega }_o)=\frac{k_d}{{\pi }_d}.$$

$$C(x)dx={\int }_{-1}^{1}I(x,{\mu }^{\prime })G({\mu }^{\prime }){\mu }^{\prime }\frac{dx}{{\mu }^{\prime }}d{\mu }^{\prime }.$$

$$\frac{c}{2}C(x)=\frac{c}{2}{\int }_{-1}^{1}I(x,{\mu }^{\prime })G({\mu }^{\prime })d{\mu }^{\prime }$$

$$\left(\mu \frac{\partial }{\partial x}+1\right)I(x,\mu )=\frac{c}{2}{\int }_{-1}^{1}I(x,{\mu }^{\prime })G({\mu }^{\prime })d{\mu }^{\prime },$$

$$\left(\mu \frac{\partial }{\partial x}+1\right)I(x,\mu )={\int }_{-1}^1\Psi ({\mu }^{\prime })I(x,{\mu }^{\prime })d{\mu }^{\prime }.$$

$$\Psi (\mu )=\frac{c}{2}G(\mu ).$$

$${\int }_0^1\Psi (\mu )d\mu =\frac{c}{2}\le \frac{1}{2},$$

$$\Psi (\pm 1)=\{\begin{array}{l}\infty ,(1\le d<3)\\ c/2,(d=3)\\ 0,(d>3).\end{array}$$

$$C(x)=C_0(x)+c{\int }_0^{\infty }K(x-x^{\prime })C(x^{\prime })d{x}^{\prime },$$

$$K(x)=\frac{1}{2}{E}_1(|x|)=\frac{1}{2}{\int }_1^{\infty }\frac{{\text{e}}^{-|x|t}}{t}dt$$

$$K(x)=\frac{1}{2}{\text{e}}^{|x|}$$

$$K(x)=\frac{K_0(|x|)}{\pi }$$

$$K(x)=\frac{1}{2}\int_{0}^{1}\text{e}^{-|x|/\mu}\frac{1}{\mu}G(\mu)d\mu=\frac{1}{2\pi}\Gamma\left(\frac{d}{2}\right)G_{1,3}^{3,0}\left(\frac{x^{2}}{4}\left|\vphantom{\frac{1}{2}}\right.\begin{array}[]{c}\frac{d-1}{2}\\ 0,0,\frac{1}{2}\\ \end{array}\right).$$

$${\int }_{-\infty }^{\infty }K(x)dx=1$$

$$K(x)={\int }_0^{\infty }h(s){\text{e}}^{-|x|s}ds$$

$$h(s)=\frac{1}{2}\Theta (s-1)G(1/s)/s,$$

$$\underset{x\to 0}{\lim }K(x)=\infty .$$

$$\tilde{K}(t)\equiv {\int }_{-\infty }^{\infty }K(x){\text{e}}^{ixt}dx={}_2{F}_1\left(\frac{1}{2},1;\frac{d}{2};-t^2\right).$$

Full text

1 Introduction

Linear transport theory [1, 2, 3] has a long history of utility in many fields for modeling the motion of particles and waves in random media and is most often concerned with transport in a spatially three-dimensional (3D) domain (even when symmetries in the solution reduce the equations to a single spatial coordinate, prompting a “1D” label). It is, however, occasionally useful to consider Euclidean domains apart from $\mathbb{R}^{3}$. For infinite homogeneous media, monoenergetic problems in $\mathbb{R}^{d}$ have been mostly solved [4]. In this paper, we solve the classic problem of diffuse reflection from an isotropically-scattering half space in $\mathbb{R}^{d}$ for general $d\geq 1$. The familiar $H$ functions, rigorous diffusion lengths, extrapolation constants, and singular eigenfunctions of radiative transfer and neutron transport are given new mathematical context as special cases in a general family of solutions expressed in closed form using hypergeometric functions. We also derive a number of new analytic results for the “Flatland” ($d=2$) domain, which is used in a wide number of applied settings.

1.1 Motivation and related work

There are both practical and theoretical reasons to consider transport problems in general dimension. The 1D rod has long been a useful domain for transport research and education [5, 6, 7, 8]. The same can be said for the 2D “Flatland” domain [9], where it is possible to visualize the entire lightfield with 2D images [10].

Transport in both the rod and Flatland settings finds numerous application in practice. The rod model is equivalent to the two-stream approximation in plane-parallel atmospheric scattering [11, 12], which is still a common method of solution for radiation budgets [13]. Transport in Flatland also has many real-world applications, such as sea echo [14], seismology [15], animal migration [16, 17], and wave propagation and diffraction in plates and ice [18, 19, 20]. Also, planar waveguides comprised of dielectric plates with controlled or random patterns of holes lead to 2D transport and have proven useful for studying engineered disorder [21]. Similarly, bundles of aligned dielectric fibers, such as clumps of hair or fur, can also be treated with a Flatland approach [22, 23, 24, 25], where it is common to employ an approximate separable product of 1D and 2D solutions [26]. Reactor design also makes use of such 2D/1D decompositions [27].

Going beyond 3D, higher-order dimensions occasionally find application in practice. Exact time-dependent solutions in 2D and 4D have been combined to approximate the unknown 3D solution for the isotropic point source in infinite media [28], and later applied to a time-dependent searchlight problem using the method of images [29]. In the study of cosmic microwave background radiation it has even been considered to change dimension over the course of a single random flight [30].

The study of transport problem in the case of general dimension can reveal how dimension $d$ impacts various aspects of the solution. New insights about 3D transport have been found by identifying correspondences between problems of differing configurations and dimensionalities [31, 32, 33], which often appear unexpectedly. Most investigations of this nature have considered only infinite media. In this work, we identify some new exact correspondences regarding anisotropic scattering in a 3D half space.

The literature on stochastic processes includes a number of fully general studies in $\mathbb{R}^{d}$ (see [34, 4, 35, 36, 37, 38] and references therein) that reveal the influence of dimension. Such studies often permit $d\geq 1$ to be a free real parameter [39] by using the general surface area $\Omega_{d}$ of the unit sphere $S^{d-1}\equiv\{x\in\mathbb{R}^{d},||x||=1\}$ in $d$ dimensions,

[TABLE]

as expressed using the Gamma function [40], and taking the meaning of a non-integral dimension flight abstractly and “by definition” from the equations that follow. Where possible, we consider the same freedom on $d$ in this paper.

For infinite homogeneous media, Green’s functions for monoenergetic linear transport with isotropic [28, 41, 42, 38, 43, 44] and anisotropic [4, 45, 46] scattering are known in domains apart from 3D. The unidirectional point source was also considered in Flatland [47]. In these infinite domains, the non-universal role of diffusion as a “rigorous asymptote” of the full solution for general dimension with absorption was observed [48, 43]. We expand on these findings for the isotropic scattering case, and find a simple algebraic condition for diffusion asymptotics to arise.

For the case of bounded domains, isotropic scattering in a Flatland half space has been solved in a number of works [18, 19, 49, 25]. Slab geometry [23] and layered problems [50] in Flatland have also been solved. The study of inverse problems in plane-parallel domains of general dimension has been considered in a number of works (see, for example, [51]). We expand on these solutions by considering general dimension and by producing the singular eigenfunctions, whose orthogonality properties allow derivation of the moments of the internal scalar flux and angular distributions. These moment derivations complement the mean, variance and general moments previously produced for 3D [52, 53] and are useful for forming approximate searchlight approximations [54] and for guiding Monte Carlo estimators towards zero-variance [55].

1.2 Outline

Section 2 provides a discussion of general dimensional ($d$D) radiometry. The integral and integrodifferential transport equations are given in Section 3, together with the dispersion equation and related eigenvalues. The albedo problem is solved in Section 4 using Chandrasekhar $H$-functions. This section also discusses evaluation strategies for the $H$ functions, escape probabilities from within the half space, and a generalized enhancement factor for coherent backscattering. Case’s method is applied to the general problem in Section 5, and used to derive equations for spatial and directional moments. Extrapolation distances for the Milne problem for a half-space and two adjacent half spaces are derived in Section 6. We derive some universal properties of half-space transport in Section 7 before concluding in Section 8. A closed-form equation for the discrete eigenvalues is given in A and additional closed-form equations for the $H$ function are given in B. Monte Carlo sampling and connections between anisotropic scattering and nonclassical transport are given in C and D, respectively.

2 Radiometry in $\mathbb{R}^{d}$

Previous studies of the transport equation for $d$-space [56, 57] implicitly assumed the generalization of standard radiometric quantities, which we briefly review here, along with a derivation of the generalized Lambertian bidirectional reflectance distribution function (BRDF).

2.1 Radiance

Let us consider time-independent mono-energetic specific intensity (radiance) $I(\mathbf{x},\omega)$ at a position $\mathbf{x}\in\mathbb{R}^{d}$ such that the rate of energy $dE$ flowing across a surface element of area $d\sigma$ in directions comprising a solid angle $d\omega$ about direction $\omega\in S^{d-1}$ is

[TABLE]

where $\theta$ is the angle between $\omega$ and the outward surface normal of $d\sigma$. With $\Omega_{d}$ “steradians” in the unit sphere $S^{d-1}$, the integral of an angularly-uniform $(I(\mathbf{x},\omega)=1)$ field of unit radiance is

[TABLE]

and the isotropic phase function is

[TABLE]

2.2 Axial symmetry

In plane-parallel problems with axial symmetry the radiance is assumed to be symmetric in all but one axis and rotationally invariant about that axis, allowing us to express the transport equations over an integrated radiance that depends on a single position, the optical depth $x$, and a single direction parameter $\mu=\cos\theta$, the cosine with respect to the depth axis. The parameter $\mu$ indexes a pair of directions in Flatland, a cone of directions in 3D, and higher-dimensional hypercones for $d>3$.

We choose a definition of integrated radiance $I(x,\mu)$ that is constant in $\mu\in[-1,1]$ whenever $I(x,\omega)$ is constant in $\omega\in S^{d-1}$ by introducing the appropriate angular measure $G(\mu)d\mu$. The total rate of radiant energy flowing across $d\sigma$ at depth $x$ confined to directions with cosines in $d\mu$ about $\mu$ is then

[TABLE]

We normalize $G(\mu)$ such that

[TABLE]

a uniform integral over the sphere $S^{d-1}$. The function $G(\mu)$ that satisfies these conditions is ([4], Eq.(33))

[TABLE]

For integer dimensionalities $d=2$ to $7,$ $G(\mu)$ is

[TABLE]

Scattering in a 1D rod is included in this definition by noting [28]

[TABLE]

The probability that a single photon leaves an isotropic collision into $d\mu$ about $\mu$ is $(c/2)G(\mu)d\mu$, where $0<c\leq 1$ is the single-scattering albedo, and so the integrated radiance about $x$ arising from isotropic collisions happening at a unit rate is $I(x,\mu)=c/2$.

2.3 Uniform diffuse illumination

In additional to unidirectional illumination of the half space, we will also consider the case of uniform diffuse illumination, which we define to be uniform radiance in all directions arriving from the hemisphere to a given surface patch. Such a source produces flux across the patch with an intensity proportional to $I(\mu)\mu$. To produce a unit flux across a patch of unit area, we require a normalization constant such that a $\mu$-weighted integral over the hemisphere is $1$,

[TABLE]

The uniform diffuse boundary source condition is thus

[TABLE]

or $I(0,\omega)=1/\pi_{d}$ where

[TABLE]

2.4 Bidirectional reflectance distribution function

The BRDF gives the radiance $f_{L}(\omega_{i},\omega_{o})$ leaving surface area $d\sigma$ in direction $\omega_{o}$ due to a unit radiance arriving at $d\sigma$ from direction $\omega_{i}$. This form of expressing the diffuse reflection law is convenient for image synthesis [58] and for comparing the behaviour to other known BRDFs. Of particular interest is the Lambertian BRDF with total diffuse albedo $0<k_{d}\leq 1$, whose generalization to arbitrary dimension is

[TABLE]

3 Transport equations

We now review the plane-parallel transport equations for isotropic scattering in a halfspace in $\mathbb{R}^{d}$. Energy balance in an infinitesimal slab in plane geometry with axial symmetry yields a transport equation of one spatial and one angular variable. If the intensity distribution arriving at the slab of thickness $dx$ is given by $I(x,\mu^{\prime})$, the flux in direction $\mu^{\prime}$ crossing $dx$ is proportional to $\mu^{\prime}G(\mu^{\prime})$, and the track lengths extending through the slab are $dx/\mu^{\prime}$, so the rate of photons entering collisions within $dx$ is

[TABLE]

The inscattered contribution to $I(x,\mu)$ is thus

[TABLE]

and the full integrodifferential form of the transport equation is then

[TABLE]

which reduces to the familiar 3D form with $G(\mu)=1$ and the Flatland equation [22, 19, 25, 49] with $G(\mu)=2/(\pi\sqrt{1-\mu^{2}})$.

Equation (15) is a “pseudo problem” of the form studied by Chandrasekhar (Section 89 of [1]). In his notation,

[TABLE]

Chandrasekhar considered pseudo problems in relation to anisotropic scattering in a 3D half space. Multiple pseudo problems with polynomial characteristic functions $\Psi_{i}(\mu)$ arise in each case, and their related $H$ functions appear in the exact solution. No individual pseudo problem on its own corresponds to a complete transport problem, hence the label. For isotropic scattering in $d$-space, however, we see a single pseudo problem does describe the complete problem. Comparing Eqs.(15) and (16), we find the characteristic function for our problem to be

[TABLE]

For $d>1$, $\Psi(\mu)$ is an even, non-negative function satisfying

[TABLE]

provided $0\leq c\leq 1$. Further, for $d>1$, $\Psi(\mu)$ is also regular on $(-1,1)$. At the boundaries,

[TABLE]

Busbridge studied a very general class of pseudo problems, relaxing the assumption of polynomial characteristic function. With the above conditions satisfied, we can apply the findings of Chapter 2 of Busbridge [59].

Before solving Eq.(15), we consider the related integral equation for the collision-rate density $C(x)$ at optical depth $x$ in the half space. This can be formed by integrating the total attenuated intensity at $x$ arriving from collisions at each depth $x^{\prime}$ inside the half space,

[TABLE]

where $C_{0}(x)$ is the forcing function (the collision-rate density of first collisions from any external source, in this case). The kernels $K$ are symmetric and account for the total collision rate density at optical depth $x$ arising from energy that leaves a collision from a hyperplane at depth $x^{\prime}$.

Eq.(20) is an integral equation of the Wiener-Hopf (W-H) kind, named after the authors who first solved it for isotropic scattering in 3D. That original W-H equation was first posed by Chwolson [60], who considered the Schwarz-schild-Milne (exponential integral) displacement kernel

[TABLE]

in his study of the translucent appearance of milk glass111112 years later, the computer rendering of a glass of milk with multiple scattering was one of the iconic images in a seminal paper [61] that sparked the subsurface revolution in film rendering and earned the authors an Academy award..

For 1D and 2D, the kernels are also already known, and their W-H equations have been studied. The Picard-Lalesco kernel

[TABLE]

describes exponential flights in a rod. The MacDonald / Hankel kernel

[TABLE]

describes isotropic scattering in Flatland [49], and has also appeared in studies of wave diffraction problems [62, 63]. It is interesting that Fock, Case and Krein each considered the three kernels (21), (22), and (23) in papers [18, 64, 65] on general W-H methods, but did not explicitly identify them as pertaining to isotropic scattering in 3D, 1D and 2D, respectively. To the best of our knowledge, they have not previously been shown as members of the same unified family. We show this now, expressing the general kernel in terms of the plane-geometry measure for $d$-space,

[TABLE]

Here, $G_{1,3}^{3,0}$ is a Meijer $G$ function. Again, we have used the assumption of isotropic scattering (which can be lifted, at considerable complexity, and won’t be treated here) and also that the free-path distribution between collisions $p_{c}(x)$ is an exponential, $p_{c}(|x|/\mu)=\text{e}^{-|x|/\mu}$, which can be easily generalized for the case of complete-frequency redistribution in line formation [66] and non-classical media with non-exponential free paths [44] (see D.2).

The kernels in Eq.(24) are positive symmetric normalized displacement/convolution kernels

[TABLE]

of the Laplace type [66], expressible as

[TABLE]

where

[TABLE]

using the Heaviside theta function $\Theta(x)$. The kernels are singular at $x=0$ for $d>1$,

[TABLE]

The Fourier transform $\tilde{K}(t)$ of the kernels will play a central role in solving the albedo problem and can be expressed for the general case $d\geq 1$ using the hypergeometric function ${}_{2}F_{1}$ [40], by taking the Fourier transform of Eq.(24) and exchanging the order of integration,

[TABLE]

The common cases of $d\in\{1,2,3\}$ reduce to the familiar set

[TABLE]

of Fourier transforms of the Picard, MacDonald and Schwarz-schild-Milne kernels, respectively.

We see that a change in dimension $d$ amounts to a change of kernel $K(x)$ and related characteristic function $\Psi(\mu)$, and these two functions completely characterize the problem. The solutions that follow will have much in common with analogous variations of $K$ and $\Psi$ that arise when considering general phase function [1], reflectance conditions at the boundary [67], and line-formation and other energy-dependent problems [68, 69, 70]. As such, we will rely on general studies of W-H equations [18, 59, 65, 68, 71, 72, 70, 66].

3.1 Dispersion equation and eigenvalues

For infinite, half space or slab geometry problems, the solutions all depend on the eigenspectra of the transport kernel $K(x)$, which can include both a continuous and a discrete component. We review the eigenvalues now, and their conditions for existence, before solving the general albedo problem.

The discrete eigenvalues, when they exist, are real zeros $\nu_{0}$ of the dispersion function $\Lambda(z)$ [59], which is related to the Fourier transform of the kernel by

[TABLE]

or, equivalently, from the characteristic function $\Psi(\mu)$,

[TABLE]

From Eq.(29) we have the generalized dispersion equation in terms of a hypergeometric function

[TABLE]

in agreement with previous derivations in infinite spherical geometry [48, 28, 38]. By known properties of ${}_{2}F_{1}$ [40], $\Lambda(z)$ satisfies the differential equation

[TABLE]

Here we see $d=3$ as a special case, the unique dimension where the $z^{3}\Lambda^{\prime}(z)$ term vanishes.

The discrete eigenvalues fall into three cases [59], based on dimension $d$ and absorption $c$. For the case of conservative scattering $c=1$, double zeros at infinity arise for all $d\geq 1$. This follows immediately from the limit as $z\rightarrow\infty$,

[TABLE]

and from $\Lambda^{\prime}(\infty)=0$. For absorbing media, $0<c<1$, given the symmetry and non-negativity of $K$, there will be either [math] or $1$ pairs of real eigenvalues $\pm\nu_{0}$, satisfying $\Lambda(\pm\nu_{0})=0$ [65].

For $d\leq 3$, $\Psi(1)\neq 0$ and so $\Lambda(z)$ always has a real root $\nu_{0}>1$. For $d>3$, discrete eigenvalues will only exist when $\Lambda(1)<0$ [59]. After observing that

[TABLE]

we find that $\Lambda(z)$ admits a finite zero $\nu_{0}>1$ if and only if

[TABLE]

This condition simplifies several prior observations and conditions for diffusion modes disappearing in dimensions $d>3$ [48, 43]222It was incorrectly reported in [43] that the number of discrete eigenvalues increases past $d=4$. and is new, to the best of our knowledge. Closed form expressions are known for $d\in\{1,2,4,6\}$ [43],

[TABLE]

Mathematica is able to find the roots in 8D and 10D in closed form but we omit these bulky expressions for space reasons.

In 3D, the eigenvalues always exist, the dispersion equation being

[TABLE]

In odd dimensionalities $d\geq 5$, the discrete eigenvalues are (like in 3D) also solutions of transcendental equations of increasing complexity, such as in 5D,

[TABLE]

In A we derive a general closed-form expression for $\nu_{0}$ for any dimension $d\geq 1$.

A related function that plays an important role in the solution of half space problems is

[TABLE]

with $\mathcal{P}$ indicating the Cauchy principal value of any integral over $\nu$ or $\mu$ must be taken. For $\nu\in[-1,1]$ the principal value integral can be expressed in closed form,

[TABLE]

which simplifies to

[TABLE]

for $d\in\{1,2,3,4\}$, respectively. Equation (40) can also be written

[TABLE]

which shows how $c=(d-3)/(d-2)$ is a special value. Again, using known properties of ${}_{2}F_{1}$, it is straightforward to find the differential equation satisfied by $\lambda(z)$,

[TABLE]

with conditions

[TABLE]

From the Plemelj-Sohotski formula[73] the boundary values of $\Lambda(z)$ are

[TABLE]

from which it follows that

[TABLE]

Also, $\Lambda(0)=\Lambda^{\pm}(0)=1$. It is straightforward to verify that Eq.(45) satisfies the differential equation in Eq.(34).

Both $\Lambda(z)$ and $\lambda(z)$ are shown in Figure 1 for a variety of dimensions $d$. For the illustrated case of $c=0.63$ we see that there are no discrete roots $\nu_{0}>1$ of $\Lambda(\nu_{0})$ for $d>4$ and also that $\lambda(\nu)$ admits either [math], $1$, or $2$ roots $0<\nu<1$.

4 The albedo problem

Let us now consider a unit current azimuthally-symmetric plane-parallel illumination arriving along direction cosine $\mu_{\ell}$, given by

[TABLE]

for $\mu$ and $\mu_{\ell}\in[0,1]$. The density of initial collisions inside the half space at optical depth $x\geq 0$ is then,

[TABLE]

which is independent of dimension. Here, $C_{0}(x)dx$ is the rate of particles entering their first collision within $dx$ about $x$.

The inhomogeneous integral equation (20) with the forcing function $C_{0}(x)$ in Eq.(48) defines the $dD$ albedo problem with a unidirectional source. Its solution can be found in a number of ways, and directly yields the internal distribution and, indirectly, the emerging distribution.

Once the collision rate density $C(x)$ is found, the scalar flux of particles in flight inside the medium is known immediately because we are assuming a classical medium with no correlation or memory between scattering events. Given our parametrization of optical depth $x$, the scalar flux is proportional to $C(x)$ by a unit constant, the mean-free path, which is simply a change of units.

The integrated radiance at the boundary is

[TABLE]

and the rate of photons leaving the medium along directions in $d\mu$ about $\mu$ is $I(0,\mu;-\mu_{\ell})\mu G(\mu)d\mu$. Eq.(49) is simply a reduced-intensity calculation of radiance leaving collisions at depth $x$, using the exponential Beer-Lambert law. The quantity $(c/2)C(x)G(\mu)d\mu$ is the integrated radiance leaving collisions at $x$ into $d\mu$ about $\mu$, $\exp(-x/\mu)$ is the Beer-Lambert calculation, and $(1/\mu)dx$ is the source measure at depth $x$; i.e. the length of the line segment tilted by $\mu$ inside the slab of thickness $dx$ from which the source of collided photons arises.

Equation (49) gives the law of diffuse reflection for the half space and shows that the emerging distribution is related to the internal distribution by a Laplace transform. The Laplace transform of the internal distribution is given in terms of $H$-functions, which we consider next.

4.1 $d$D Chandrasekhar $H$-functions

Given the characteristic function $\Psi(\mu)=(c/2)G(\mu)$ for isotropic scattering in $d$D (Section 3), the related $H$ functions satisfy the integral equation [59],

[TABLE]

for $0\leq\mu\leq 1$ or, more generally

[TABLE]

Regardless of dimension $d\geq 1$ or absorption $0<c\leq 1$, the solution of Eq.(50) is known in closed form by the Fock/Chandrasekhar equation [18, 59]

[TABLE]

Figure 2 illustrates values of $1/H(z)$.

Let us pause for a moment to consider the significance of the appearance of $\Lambda(i/k)$ in Eq.(52). By Eq.(31),

[TABLE]

which is the Fourier-space Neumann-series Green’s function for the isotropic plane source in infinite geometry. So we see an exact infinite-space solution inside the $H$ function expression for the half space,

[TABLE]

Ivanov [66] noted a similar relationship for general displacement kernels $K$. From his analysis we also have, for all $d\geq 1$,

[TABLE]

4.1.1 $H$-function moments

The moments of the $H$ functions can form an additional integral equation for $H$, and also arise in later expressions for the internal and emerging distributions of the albedo problem (and related extrapolation distances), and so we look at these now.

If Eq.(51) is expanded around infinity, one obtains

[TABLE]

where $\alpha_{n}$ are moments of the $H$-function defined by

[TABLE]

In her general study of pseudo problems, Busbridge [59, 74] found that the $\alpha_{n}$ moments of Eq.(57) satisfy the recurrence equations

[TABLE]

where moments related to the characteristic function are given by

[TABLE]

For our present study of isotropic scattering in $d$-space, we find

[TABLE]

which, for $n\in\{0,1,2,3\}$, is

[TABLE]

In Flatland, $g_{2n}$ reduces to

[TABLE]

differing from the familiar 3D case, $g_{2n}=(2n+1)^{-1}$.

For nonconservative scattering $0<c<1$, the odd moments must be evaluated numerically, and the even moments are given from the recurrence relations. For example,

[TABLE]

For conservative ($c=1$) scattering, it is the even moments which must be numerically evaluated, and the odd moments deduced,

[TABLE]

Figure 3 and Table 1 illustrate a peculiar, but not exact, pairing of moments for different dimensions that we cannot explain and leave as an interesting area for future investigation.

Equation (62) can be used to rewrite Eq.(50) as

[TABLE]

which is the form that should be used if evaluating $H$ directly from an integral equation [1].

4.1.2 $H(z)$ calculation methods

Given their central role in the solutions that follow, we now consider important details regarding the uniqueness and evaluation strategies for the $H$ functions.

Equation (52) is the unique solution of Eq.(51) unless there is a finite discrete eigenvalue $\nu_{0}$ of the dispersion Eq.(33). In the latter case, there is one other non-physical solution [59], which is not relevant to our problem. Equation (52) is the unique solution of Eq.(66) in all cases [59].

Fox [75] considered the expression of $H$ with general characteristic $\Psi$ as the solution to a Riemann-Hilbert problem involving the function

[TABLE]

where $\theta(0)=0$ and $0\leq\theta(t)\leq\pi$. When $(d-3)/(d-2)<c<1$ and one real eigenvalue $\nu_{0}>1$ exists, a useful extension of Fox’s approach is [71, 76]

[TABLE]

For $c\leq(d-3)/(d-2)$, with no eigenvalues $\nu_{0}>1$, the simpler result is [71]

[TABLE]

We noted improved numerical stability over Eq.(52) for $d>4$ when using these last two forms of $H(z)$. In practice, when $\lambda(\nu)$ has zeros, negative $\tan^{-1}$ results need to be detected manually and remapped to ensure $0<\theta(t)<\pi$ when performing numerical evaluation of Eqs.(68) and (69).

Carlstedt and Mullikin [71] declare that when there are no discrete roots $\nu_{0}>1$ of $\Lambda(z)$ that there are then no roots of $\lambda(\nu)$ for $\nu\in[-1,1]$, which we found to not hold in general (Figure 1 shows that $\lambda(\nu)$ will admit $0,1$ or $2$ roots, depending on $d$ and $c$), but we noted no issues in applying Eqs.(68) and (69), provided $\theta(t)$ was strictly non-negative.

Additional closed-form expressions for computing $H(z)$ are given in B.

4.2 Law of diffuse reflection

We can derive the law of diffuse reflection for the half space by solving for the internal collision rate density $C(x)$ due to external unidirectional illumination along cosine $\mu_{\ell}$, the solution of the W-H equation (20) with source term (48). This collision rate can be found using the Green’s function for an internal isotropic plane source at depth $x_{0}$ (with $C_{0}(x)=\delta(x-x_{0})$).

The Green’s function $\mathbb{G}(x,x_{0})$ is a source function such that $\mathbb{G}(x,x_{0})dx$ is the rate of photons leaving collisions (or the source directly) from depths $dx$ about $x$. If we define the Laplace transform

[TABLE]

then we have, from Ivanov ([66], Eqs. (19) and (21)), that the double Laplace transform of the Green’s function is

[TABLE]

Before considering the external source, we first note a number of exact properties that relate to the life of a photon in the half space. From $\mathbb{G}(x,x_{0})$, we can find the collision rate density due to isotropic emission (or leaving a collision) at depth $x_{0}$. We convert the source function $\mathbb{G}(x,x_{0})$ to a collision rate density by removing the Dirac delta for direct emission (since this is not a collision), which is

[TABLE]

and then apply a factor $1/c$, to convert densities for leaving collisions into densities for entering collisions,

[TABLE]

The double Laplace inversion of $\bar{\bar{C}}(s;s_{0})$ gives the collision rate density $C(x;x_{0})$ at any position $x\geq 0$ inside the half space, due to an isotropic plane source at depth $x_{0}$. We also easily have the mean number of collisions, by taking the Laplace inversion of $\bar{C}(0;s_{0})$, which is, by Eq.(55),

[TABLE]

Of those collisions, $1-c$ will absorb the photon, so the escape probability $p(x_{0})$ after isotropic emission (or leaving a collision) at $x_{0}$ is then

[TABLE]

For unidirectional illumination along cosine $\mu_{\ell}$, the collision-rate density inside the half space will be

[TABLE]

After combination of Eqs.(49) and (76), we find

[TABLE]

which is the generalized law of diffuse reflection for a half space of general dimension. The probability that a photon arriving along cosine $\mu_{\ell}$ escapes the half space along a direction within $d\mu$ of $\mu$ is $\mu I(0,\mu;-\mu_{\ell})G(\mu)d\mu$, and so the total albedo of the half space is

[TABLE]

where we have used Eq.(66). Equations (77) and (78) are in agreement with previous derivations for Flatland [49] and show how the familiar expressions for 3D are universal over dimension $d\geq 1$, with all variation due to the $H$ function and the measure $G(\mu)d\mu$. The variation of the emergent distribution and albedo with respect to dimension $d$ is illustrated in Figures 4 and 5.

4.2.1 Low-order scattering

The once- and twice-scattered portions of the reflection law and albedo shine further light on the structure of the solutions and can also provide accurate approximations for high absorption. These are found via Taylor expansions about $c=0$ of Equations (77) and (78). The once-scattered reflection law reduces to

[TABLE]

In the Taylor expansion of the reflection law for the twice-scattered emergent distribution, we encounter an integral that we observe is equal to the Laplace transform of the kernel,

[TABLE]

The twice-scattered component of the reflection is then,

[TABLE]

We found the general form of the Laplace transform to reduce to

[TABLE]

Special cases include the known results for Flatland [49]

[TABLE]

for 3D [77]

[TABLE]

and our new result for 4D,

[TABLE]

Figure 6 shows Monte Carlo validation of the emergent distributions for a 4D half space.

The escape probability after exactly one scattering also involves the Laplace transform of the kernel,

[TABLE]

An expression for $R(\mu_{\ell}|1)$ in terms of the characteristic function is also known [72]

[TABLE]

4.2.2 Uniform diffuse illumination

We now consider the case of uniform diffuse illumination of the half space, with boundary condition given by Eq.(10). Integrating the diffuse reflection law (77) against the source function

[TABLE]

and applying Eq.(66), we find the emergent distribution

[TABLE]

with albedo

[TABLE]

The emergent distribution $I(0,\mu)$ is proportional to Eq.(78), as it must be, by the optical reciprocity theorem. For the conservative case, $c=1$, $R$ reduces to Eq.(9), which is $1$ by construction.

We find a number of analytic results for the single- and double-scattered components of the albedo for diffuse illumination, from a Taylor expansion about $c=0$ of Eq.(91). For the probability of escape after one collision, in Flatland we find

[TABLE]

for 3D[78]

[TABLE]

and for $d\geq 4$

[TABLE]

It appears that for even dimensionalities $d\geq 4$, the escape probabilities are rational multiplies of $c$, including the remarkably simple single-scattering albedo for 4D, $R_{1}=c/5$. The double-scattering albedo for diffuse illumination in 4D is also a rational factor of $c^{2}$, $R_{2}=4c^{2}/35$.

4.3 Coherent backscattering

In applications such as remote sensing it can be important to consider the effects of coherent backscattering. Under the conditions of weak localization, the reflectance will be amplified in the backward direction by an enhancement factor $\eta$ that reaches a maximum for normal incidence and $c=1$. The maximal enhancement factors for 2D and 3D half spaces are [22, 24]

[TABLE]

which generalize to arbitrary dimension $d$ using Eq.(52), producing a monotonically increasing factor as dimension $d$ increases (Figure 7), beginning at $\eta=7/4$ for a 1D rod and approaching $\eta\rightarrow 2$.

5 The albedo problem by Case’s method

5.1 Introduction

In this section we consider Case’s approach [79] to solving the $d$D albedo problem to illustrate its relationship to the Wiener-Hopf and resolvent approaches. To solve Eq.(15) with the classic Case-style eigenmode method[79], we separate variables with the ansatz

[TABLE]

for eigenvalues in the spectrum $\nu\in\sigma$, to obtain

[TABLE]

If we impose the normalization condition for the eigenmodes to be

[TABLE]

then $\sigma=\{\nu\in[\,-1,1]\cup\pm\nu_{0}\}$ is the eigenvalue spectrum and the eigenmodes $\phi(\nu,\mu)$ satisfy the equation

[TABLE]

identical to the form in a 3D domain (relative to our angular measure $G(\mu)$). From Eq.(99) the discrete eigenmodes, when they occur, satisfy

[TABLE]

with roots $\nu_{0}$ obtained from the dispersion function (33), which we have summarized in Section 3.1. The eigenmodes for the continuum are

[TABLE]

From Eqs.(7) and (98) it follows that the continuum eigenmodes also satisfy

[TABLE]

In this section, the collimated incident and diffuse illuminations for the albedo problem are selected to be

[TABLE]

with $I(x,\mu;\mu_{\ell})$ and $I_{k}(x,\mu)$ tending to [math] as $x\rightarrow\infty$. The special case $k=-1$ corresponds to a unit current illumination. The notation $\psi(\mu)$ will be used henceforth to denote either surface illumination.

The constraint $I(x,\mu)\rightarrow 0$ as $x\rightarrow\infty$ forces the eigenmode expansion for the albedo problem to be written as

[TABLE]

where $\sigma+=\{\nu\in[\,0,1]\,\cup\,\nu_{0}\}$ defines the half-spectrum subset of the full spectrum of eigenvalues $\sigma$. From Eqs.(103), (104), and (105) it follows that the expansion coefficients $A(\nu_{0})$ and $A(\nu)$, $\nu\in[0,1]$, are to be determined from

[TABLE]

We follow the approach in [80] and [53]333Tables 1, 3 and 4 of reference [53] have numerical errors. The values for the column labeled $j_{ratio,1}(\mu_{0})$ in Table 1 are all incorrect. The correct values for Table 3 with $c=0.7$ are $<\!\!x^{2}(\mu_{0})\!\!>\,\,=3.36647$ and $<\!\!x^{3}(\mu_{0})\!\!>\,\,=12.45400$ for $\mu_{0}=0.9$ and

$<\!\!x^{2}(\mu_{0})\!\!>\,\,=3.83320$ and $<\!\!x^{3}(\mu_{0})\!\!>\,\,=14.86410$ for $\mu_{0}=1.0$. The top right value of Table 4 should be $1.91257$. to construct a Chandrasekhar $H(\mu)$ function by forcing the eigenfunctions $\phi(\nu,\mu)$ to obey a half-range-in-$\mu$ transport equation analogous to Eq.(97),

[TABLE]

Equation (97) then forces the constraint

[TABLE]

With the substitution of $\phi(\nu,\mu)$ from Eqs.(100) and (101) into Eq.(108), the following equations are valid,

[TABLE]

and

[TABLE]

This suggests we construct $H(\mu)$ by considering in the complex plane the equation

[TABLE]

and examining the analyticity properties of $H(\mu)$, $0\leq\mu\leq 1$. The factors $\Lambda(z)H(z)$ and $H(-z)$ are both continuous across $(-1,0)$; similarly, $\Lambda(z)H(-z)$ and $H(z)$ are continuous across $(0,1)$ so $\Lambda(z)H(z)H(-z)$ is analytic along $(-1,1)$. The points $z=\pm 1$ can be included so $\Lambda(z)H(z)H(-z)$ is analytic in the entire plane and from Liouville’s theorem approaches a constant. With $H(0)=\Lambda(0)=1$ the Wiener-Hopf identity results,

[TABLE]

with $H(z)$ satisfying Eq.(51) subject to the constraint $1/H(-\nu_{0})=0$ imposed by Eq.(109). Thus, for $H(\mu)$, $0\leq\mu\leq 1$, Eq.(50) is recovered. Equation (66) is needed, for example, to show that multiplication of Eq.(107) by $G(\mu)H(\mu)$ and integration over $\mu\in[0,1]$ gives

[TABLE]

after use of Eq.(62).

The results of Eqs.(108) and (113) can be generalized, following a partial fraction analysis along with use of Eq.(57), to show that

[TABLE]

if the convention $\sum_{j=1}^{0}X_{j}\equiv 0$ is understood here and elsewhere. For $k=-1$, Eq.(108) gives $U_{0}(\nu)=1$, $\nu\in\sigma+$.

5.2 Orthogonality relations

Multiply Eq.(107) by $\nu^{-1}\phi(\nu^{\prime},\mu)$ and, in a second equation for $\nu^{\prime}$, multiply by $\nu^{\prime-1}\phi(\nu,\mu)$ and then integrate both results over $\mu\in[0,1]$ and subtract to obtain

[TABLE]

The corresponding normalization equations when $(d-3)/(d-2)<c<1$ are

[TABLE]

for $\nu,\nu^{\prime}\in(0,1)$.

The discrete normalization (116) can be derived by first differentiating Eq.(66), written for $-z$, and using Eq.(112); then for $(d-3)/(d-2)<c<1$ and with Eq.(100) and $\Lambda(\nu_{0})=0$, it follows that

[TABLE]

For $d\in\{1,2,3,4\}$ we find, respectively,

[TABLE]

Another convenient expression for $N(\nu_{0})$, generalizing the derivation of Case and Zwiefel [79] (p. 68), is found from writing the dispersion equation (33)

[TABLE]

where

[TABLE]

Differentiating (121) with respect to $c$, we find

[TABLE]

From Eq.(121), we also have $F(\nu_{0})=1/c$. Combining with Eq.(118), we find

[TABLE]

The continuum normalization of

[TABLE]

can be derived with the Poincaré-Bertrand formula and $\Lambda^{\pm}(\nu)$ from Eq.(45) yielding

[TABLE]

where $\lambda(\nu)$ is given by Eq.(40). Use of that formula enables an interchange of the order of integration from $\int_{0}^{1}d\mu\int_{\sigma+}d\nu$ to $\int_{\sigma+}d\nu\int_{0}^{1}d\mu$ when using the orthogonality relations.

With the orthogonality relations we can determine the half-range expansion coefficients $A(\nu)$ in Eq.(106) from

[TABLE]

Completeness of the eigenfunctions is assured by virtue of the closure relation

[TABLE]

for $\mu,\,\,\mu_{\ell}\in(0,1)$. This last equation can be derived as in [81, 82], again with the Poincaré-Bertrand formula, or easily confirmed by solving Eq.(127) for $\psi(\mu)=\delta(\mu-\mu_{\ell})$ and substituting the result into Eq.(106) to verify closure.

Three identities directly follow from Eq.(128). For the first, multiply by $1/\mu_{\ell}$ and integrate over $\mu_{\ell}\in[0,1]$; after interchanging orders of integration and using Eq.(108) it follows that

[TABLE]

Yet another identity can be derived if Eq.(128) is multiplied by $\mu_{\ell}^{n}d\mu_{\ell}$ and $G(\mu)H(\mu)d\mu$ and the result integrated for both variables over $[0,1]$ to obtain

[TABLE]

For $n=0$ the last result is consistent with Eq.(113).

Other equations can be derived by rewriting Eq.(66) for $-z$ and $-z^{\prime}$ and subtracting the resulting equations to obtain

[TABLE]

before specializing $z$ and $z^{\prime}$ to variables $\nu$ and/or $\mu$. After taking the appropriate limits as $z$ and $z^{\prime}$ approach eigenvalues, all the orthogonality results can be condensed into the formula

[TABLE]

where $\Xi(\nu)=0$ for $0\leq\nu\leq 1$ and 1 otherwise. Thus, we find

[TABLE]

Equation (133) leads to the reflection relation (or “albedo operator”) that allows us to conveniently express the outgoing radiation from the surface in terms of the ingoing radiation. To derive that equation, observe from Eq.(106) that

[TABLE]

so insert Eq.(133) into Eq.(134), interchange the order of integrations, and use Eq.(106) to obtain

[TABLE]

5.3 Albedo problem spatial moments

Substitution of Eqs.(103) and (104) into Eq.(127), followed by the use of Eq.(105), yields general equations for the angular intensities for the collimated and diffuse illuminations, respectively, as

[TABLE]

The objective here is to compute the $n$th spatial moments of the flux within the half-space, as given for $n=0,\!1,\!2,$ etc., by

[TABLE]

Also of interest are the ratios of fluxes that give the $n$th mean distances of travel in $d$D before absorption or escape from the half space,

[TABLE]

Substitution of Eq.(105) into either Eq.(137) or (138) gives a spatial moment

[TABLE]

After substitution of $A(\nu)$ from Eq.(127), followed by integration over $x$ and $\mu$ and use of Eq.(98), the result is

[TABLE]

that is the $d$D equivalent to the 3D Eq.(I42)444Eq.(XX) from [53] will be cited as Eq.(IXX) that has no factor $G(\mu)$ in the integral.

Use of the Dirac delta $\psi(\mu)$ of Eq.(103) in Eq.(141) gives

[TABLE]

where

[TABLE]

Multiplication of closure relation (128) by $\mu^{n+1}H(\mu)G(\mu)d\mu$ and integration over $\mu$ gives

[TABLE]

After substitution of Eq.(114) into (144) and a rearrangement of terms, it follows from Eq.(142) that

[TABLE]

with the understanding that $\sum_{j=1}^{0}X_{j}\equiv 0$.

Equations (139a) and (145) then give the recursion relation

[TABLE]

with $\left<x^{0}(\mu_{\ell})\right>\equiv 1.$ The first two equations, for example, are

[TABLE]

The results for $\left<x^{2}(\mu_{\ell})\right>$ and $\left<x(\mu_{\ell})\right>$ are useful because they can be combined to give the variance $V(x(\mu_{\ell}))$ of the spatial distribution,

[TABLE]

For a diffuse illumination, the spatial moments of Eq.(141) for a diffuse illumination can be derived for $n\geq 1$ with Eq.(145) used as a Green’s function to obtain(I55)

[TABLE]

with the starting condition

[TABLE]

Equations (139b) and (150) then lead to

[TABLE]

for $k=-1,\ 0,\ 1$, etc. We denote the variance of the optical depth of the photon under diffuse illumination by

[TABLE]

The $d$D half-space results in Eqs.(145) to (153) are the same as for the 3D half-space except the integrals in $H(\mu)$ and $\alpha_{n}$ have the additional factor $G(\mu)$. The equations for the mean and variance of the optical depth of a photon under unidirectional and diffuse illuminations are illustrated and compared to Monte Carlo simulation in Figures 8 and 9. Monte Carlo sampling methods for $\mathbb{R}^{d}$ are described in C.

5.4 Albedo problem directional moments

Moments of the direction of photons emerging from the surface of the half-space also can be derived, in analogy to the emerging directional moments for a 3D half-space[53]. Equation (135) for the collimated incident beam boundary condition of Eq.(103) immediately gives

[TABLE]

in agreement with Eq.(77) once the differences between Eqs.(47) and (103) are accounted for. We choose to weight $I(0,-\mu;\mu_{\ell})$ with the powers $\mu^{n}G(\mu)$, $n=0,\ 1,$ etc., so for the outward and inward moments of a collimated incident illumination we have

[TABLE]

The ratios $R_{n}$ for $n=0,\ 1,$ etc.​​​ ,

[TABLE]

are the fractions of the incident current propagated in the $n$th outward directional moment, with $R_{0}$ the probability that an entering photon will escape the half space.

With the partial fraction analysis of

[TABLE]

and Eqs.(66) and (156), the directional moments for a collimated surface illumination equal the 3D moments of Eq.(I37a) except for the extra factor $G(\mu_{\ell})$,

[TABLE]

so with $j_{\,in,n}(\mu_{\ell})=\mu_{\ell}^{n+1}G(\mu_{\ell})$ the ratios $R_{n}(\mu_{\ell})$ can be computed for $n=0,\ 1,$ etc. Three examples are

[TABLE]

with $R_{0}(\mu_{\ell})$ giving the albedo of the halfspace, in agreement with Eq.(78). The ratios in Eqs.(161) and (162) are identical to the corresponding 3D ratios of Eqs.(I39a) and (I39b) except for the different definitions for $H(\mu_{\ell})$ and $\alpha_{n}$. We verify these in Figure 10.

For a diffuse illumination, Eqs.(104) and (135) initially give, for $k\geq 0$,

[TABLE]

With the help of Eqs.(66) and (159), this becomes

[TABLE]

With the definitions

[TABLE]

it follows for $k\geq 0$ that

[TABLE]

We define the directional moments

[TABLE]

where again $R_{0,k}$ gives the albedo of the half-space. For uniform diffuse illumination $k=0$, we find

[TABLE]

which we verify in Figure 11.

Turning now to the special case of $k=-1$ that corresponds to a unit incident current, Eqs.(50) and (135) give

[TABLE]

The albedo of the half-space in this instance is

[TABLE]

after use of Eq.(62).

For other values of $k$, using Eqs.(166) and (167), we find

[TABLE]

which differs from $j_{\,in,n,-1}$ for 3D.

6 The Milne problem

6.1 Extrapolated endpoint for $c=1$

In this section we extend the derivation of Hopf ([83] p. 85, [84] p.225) to find the distance $z_{0}$ such that the rigorous asymptotic portion of the scalar flux in the Milne problem with absorption $c=1$, when extrapolated outside of the medium, reaches $\phi_{as}(-z_{0})=0$. With Hopf’s $\kappa(u)=\Lambda_{c=1}(1/u)$ in our notation, we find

[TABLE]

which we find to simplify further

[TABLE]

The distance $z_{0}$ is also known to be related to the second moment of the $H$ function ([59] p. 56)

[TABLE]

From the characteristic function we find

[TABLE]

From the last 2 equations, with $c=1$, we find

[TABLE]

For the 1D rod we find $z_{0}(1)=1$. For Flatland we find

[TABLE]

in agreement with [49]. For 3D, Eq.(176) matches Hopf’s equation (forgiving the typesetting error) and is best evaluated using the procedure proposed by Plazcek and Seidel [85], which we extend to higher odd dimensions. For 4D we find another closed form solution

[TABLE]

In 5D, the integral

[TABLE]

can be evaluating using Plazcek and Seidel’s approach giving

[TABLE]

With d = 6,

[TABLE]

Mathematica is able to determine a bulky but closed form solution for $z_{0}(8)$, which contains an isolated $4/\pi$, which is interesting to compare to the other even-dimension solutions. Numerical constants for the first $8$ values of $d$ are summarized in Table 2.

6.2 Extrapolated endpoint for $c<1$

With absorption, we consider Case’s approach to the Milne problem with boundary conditions

[TABLE]

where the asymptotic portion of the intensity can be written as

[TABLE]

The extrapolation distance $z_{0}(c,d)$ now depends on absorption and dimension, and the asymptotic flux satisfies $\int_{-1}^{1}I_{as}(-z_{0}(c,d),\mu)d\mu=0$, leading to

[TABLE]

The expansion coefficient $A(\nu_{0})$ is determined from Eq.(127)

[TABLE]

with $\psi(\mu)=-\phi(-\nu_{0},\mu)$. But from Eq.(132), it follows for $\nu^{\prime}\in\sigma_{+}$ that

[TABLE]

and so the combination of Eqs.(187) and (189) gives

[TABLE]

in agreement with Busbridge ([59], Eq.(23.13)).

We found numerical agreement with Eq.(190) and a previous derivation for Flatland and 3D [49]

[TABLE]

which generalizes to general dimension when $0\leq\theta(t)\leq\pi$ is determined from Eq.(67).

For monoenergetic isotropic scattering and Fresnel-

matched boundaries in $d\leq 3$, the extrapolation distance always exists, but an extrapolation distance is not universal for $d>3$, and disappears when $c<(d-3)/(d-2)$, as seen in Figure 12.

For the Milne problem emerging intensity, we find from Eqs.(135) and (189) that

[TABLE]

Because $I(0,\mu)=0$ and $G(-\mu)=G(\mu)$ for $\mu\in[0,1]$, the directional moments

[TABLE]

are, from Eqs.(108) and (113),

[TABLE]

6.3 Extrapolated endpoint for adjacent half-spaces

Problems of two adjacent half-spaces with isotropic scattering also can be analyzed[86, 87, 88, 89, 80, 90], with different absorption and scattering properties denoted by the single scattering albedo values $c_{k}=c_{1}$ for $x>0$ and $c_{k}=c_{2}$ for $x<0$. The characteristic functions now are $\Psi_{k}(\mu)=(c_{k}/2)G(\mu)$, for $k=1,\ 2$, where $G(\mu)$ again is given for $d$D in Eq.(7). Completeness of the discrete and continuum eigenfunctions,

[TABLE]

leads to the expansions

[TABLE]

with

[TABLE]

For continuity at the interface at $x=0$,

[TABLE]

where

[TABLE]

If we denote

[TABLE]

then Eq.(197) can be written as

[TABLE]

The weight function $\mu W(\mu)$ for two-media orthogonality relations can be expressed in terms of the $H$-functions for each medium[86] with

[TABLE]

The adjacent half-space orthogonality relation is

[TABLE]

and the associated normalization integrals are

[TABLE]

In a similar manner[80], other adjacent-media equations follow, as in Eq.(132).

With the exponentially growing source in $x>0$, $\mathcal{I}(\mu)=-\phi_{1}(-\nu_{01},\mu)$ and the adjacent-medium extrapolation distance $z_{0}(c_{1},c_{2},d)$ becomes

[TABLE]

Tabulated results for 3D are available[88, 80] and Figure 13 and Tables 3, 4, and 5 illustrate results for different values of $c$ and $d$.

7 Half-space and adjacent half-space properties invariant to $d$

A well known property of sums of independent identically distributed (i.i.d.) symmetric random numbers is that the mean first-passage time (number of terms in the sum) for the sum to go negative is universal and does not depend on the distribution of the individual terms. Early proofs of this theorem were given by Sparre-Andersen and Pollaczek in the 1950s [91, 92].

The depth coordinate of a photon undergoing isotropic scattering in a half space corresponds precisely to such a stochastic process, with depth displacements between collisions drawn from the transport kernel $K(x)$, which varies with $d$. Escape occurs the first time the sum of displacements goes negative. For internal sources, unidirectional illumination, and uniform diffuse illumination, we have seen that the mean number of collisions and escape probabilities for a half space depend on $d$. However, when the first displacement into the medium is also drawn from $K(x)$, by placing an isotropic point or plane source at the boundary interface, then the i.i.d. property is satisfied, and we find mean collision rates and escape probabilities that are universal and that are given by very simple equations.

7.1 Half-space invariants

Frankel and Nelson [93] considered a half space in 3D with an isotropic plane source at the boundary, and gave a derivation of the albedo, including the generating function for escape over number of collisions. Their derivation does not rely on the exponential free-path distribution between collisions, nor that the medium is three dimensional, and so, in some sense, constitutes one of the earliest proofs of this universal property of sums of random variables. Their derivation finds the universal albedo for the one-sided isotropic surface source (only emitting into the medium) to be

[TABLE]

This result agrees with our prior derivation for a one-sided isotropic plane source at the boundary, Eq.(173). Equivalently, if the first displacement is [math], which is the case under grazing illumination $\mu_{\ell}=0$, then we also find a universal albedo, directly from Eq.(78), given that $H(0)=1$,

[TABLE]

A closely-related universal property is the collision rate density at the boundary $C(0)$ due to a unit uniform isotropic source inside the medium. Davison [2] noted the relationship of this problem, by reciprocity, to that of the total collision rate inside the medium due to an isotropic plane source at the boundary. This deviates slightly from the earlier example by now considering a two-sided source at the boundary, with half of the photons escaping directly outward causing no collisions. The mean number of collisions inside the medium (and by reciprocity, $C(0)$ for the constant source problem), is now half of the number created by the one-sided isotropic source,

[TABLE]

agreeing with prior derivations for 3D [2] and 2D [49].

7.2 Adjacent half-space invariants

Let us now consider another related scenario: two adjacent half spaces with different single-scattering albedos, $c_{1}\neq c_{2}$, with an isotropic point or plane source at the boundary. With emission and scattering being isotropic, it suffices to consider only a 1D rod, with displacements given by a normalized continuous symmetric distribution $K(x)$. We can compute the mean number of collisions using Stokes’ “glass plates” analysis of interreflections between the two half rods. Let us define a reflection operator $\mathcal{R}(c)$ for a half rod for light arriving at the interface, using Eq.(204) and a related absorption operator $\mathcal{A}(c)=1-\mathcal{R}(c)$. The total absorption $\mathcal{A}_{1}$ in rod $1$ can be found by considering that the initial intensity entering rod $1$ is $1/2$ from the photons moving towards rod $1$ initially, and $(1/2)R(c_{2})$ from the photons that reflect off of rod $2$ before arriving at rod $1$. Combining these intensities, $((1/2)+(1/2)R(c_{2}))\mathcal{A}(c_{1})$ will be absorbed in rod $1$ with no further collisions. But some portion of this intensity will reflect off of rod $1$, and off of rod $2$, arriving back at rod $1$, etc. Completing the geometric series, we find,

[TABLE]

The total number of collisions in rod $1$ is then

[TABLE]

This agrees with known results for the collision rate density $C(0)$ at the boundary between two half spaces due to a uniform unit source in the half space with absorption $c_{1}$ [87], which again follows from the reciprocity theorem.

The total number of collisions $C_{2}$ in rod 2 is $C_{1}$ with $c_{1}$ and $c_{2}$ swapped, and the total number of collisions in the entire system given by

[TABLE]

To the best of our knowledge, Eq.(209) is new and contains within it two classical universal properties of linear transport. When both absorption albedos match, $c=c_{1}=c_{2}$, we form a homogeneous infinite medium with the desired mean number of collisions $C_{1}+C_{2}=1/(1-c)$. Also, when $c_{2}=0$, we recover the $\sqrt{\epsilon}$ law [66] for the half space.

We found Monte Carlo simulation to validate Eq.(209) for a variety of dimensions and free-path distributions between collisions.

8 Conclusion

We have solved the half-space albedo and Milne problems of linear transport theory in a generalized setting of Euclidean $\mathbb{R}^{d}$ space. The familiar solution of the 3D problem is given new mathematical context as a member of a new family of pseudo problems whose general solutions are given exactly in terms of hypergeometric functions. The 3D problem is a unique member of this family, the only dimension where the characteristic $\Psi(1)=c/2$ is not [math] and not $\infty$ at $\mu=1$, and the largest dimension that includes a discrete diffusion component for all absorption levels $c>0$.

In Section 89 of his book[1], Chandrasekhar concludes “$\dots$ a full discussion of the relation of the ‘pseudo-problems’ $\dots$ will throw some light on the basic structure of the theory of radiative transfer” and we have found this to be the case. The new family of pseudo problems includes a two dimensional wave propagation problem that corresponds exactly to isotropic scattering in Flatland half spaces, and an old solution of that problem has led to an exact analytic $H$ function for Flatland, which should be useful in practical settings.

That this family also includes $H$ functions arising for 3D anisotropic scattering, as given in D, was unexpected. These, and the strange, near correspondences between $H$ function moments in various dimension, seem worthy of future attention.

9 Acknowledgments

Thanks to B. D. Ganapol for pointing out the numerical errors in [53] and for tracking down a twice-scattered BRDF solution for 3D [77].

Appendix A Analytic expression for discrete eigenvalues

From Eqs.(68) and (112) with $z=0$ we follow Siewert [94] to find a closed-form expression for the discrete eigenvalues in any real dimensionality $d>1$,

[TABLE]

provided $(d-3)/(d-2)<c<1$. The term $\theta(t)$ is known analytically for all $d\geq 1$ by combining Eqs.(67) and (40) and must be restricted to be positive.

Appendix B Exact analytic forms of $H$

In our initial study of the family of $H$ functions for isotropic scattering in $\mathbb{R}^{d}$ we found analytic forms for several special cases, as well as simple relationships between different dimensionalities at different levels of absorption.

B.1 The $H$ function for the 1D rod

For the 1D rod, the $H$ function reduces to a simple expression,

[TABLE]

B.2 The Flatland $H$ function

In 2D Flatland, analytic forms have been reported [49] for normal incidence ($\mu=1$)

[TABLE]

involving a hypergeometric function ${}_{3}F_{2}$ with derivative

[TABLE]

The same ${}_{3}F_{2}$ function appears in the case of conservative scattering in Flatland [49],

[TABLE]

We now show that Flatland $H$ can be expressed analytically for all $0<\mu<1,0<c<1$ by exploiting an integral form presented by Fock [18] for a diffraction problem with the MacDonald kernel of Eq.(23). Fock found (see also [65]),

[TABLE]

If we define

[TABLE]

where $\text{Li}_{2}$ is the dilogarithm function,

[TABLE]

then we can write the analytic form of Flatland $H$,

[TABLE]

We also observe that

[TABLE]

which coincides with Eqs.(212) and (214) by noting that

[TABLE]

Similar expressions involving $\text{Li}_{2}(z)$ and ${}_{3}F_{2}$ also arise in the case of a Cauchy random flight in a 1D rod [95] and path integrals over half spaces [96]. These analytic forms have advantages relative to the numerical evaluation of integrals using quadrature and related methods. Software such as Mathematica can easily produce trusted benchmark values of any desired precision when given such analytic inputs. Using Eq.(218), we found complete agreement with the benchmark values for Flatland $H$ reported in a previous paper [49]. We caution that the above expressions are restricted to $0<\mu<1$, and do not apply to terms such as $H(\nu_{0})$, which appear in the Milne extrapolation distance equations.

The new analytic form of $H$ also permits additional exact simplified benchmark values involving Catalan’s constant $C=0.9159655942...$. To add to the previously reported result [24, 49]

[TABLE]

we observe that

[TABLE]

B.3 The 4D $H$-Function

With a similar approach to the derivation in Flatland [49] we also found a closely-related analytic form of $H$ in 4D with conservative scattering,

[TABLE]

B.4 Interdimensional correspondences

We also noticed an interesting relationship between dimensionalities two apart for specific values of absorption. Let us write Eq.(52) for $H$ in a new notation, where $f(t)$ is a general input and the argument to the logarithm,

[TABLE]

We have, by linearity of the integral and properties of the log and exp,

[TABLE]

If we now write Eq.(33) as a function with three arguments,

[TABLE]

we conjecture

[TABLE]

for $d\geq 3$ (which we verified for $d$ up to $30$). Given that

[TABLE]

we therefore conjecture that

[TABLE]

By this notation we mean

[TABLE]

This relates the $H$ function for conservative scattering at some dimensionality $d$ to the $H$ function in two dimensions higher where absorption is set to $c=\frac{d-3}{d-2}$, which is the absorption level where the discrete eigenspectra disappears. For example,

[TABLE]

and so on. These are, to the best of our knowledge, the only known analytic expressions for (or relationships between) $H$-functions that correspond to non-trivial multiple-scattering processes.

Similar unexpected correspondences have been noted in infinite $dD$ transport processes, such as Brownian motion, where a particle leaving the origin spends the same total time in a sphere in dimension $d+2$ as the first time to escape the sphere of the same radius in dimension $d$, for any $d>1$ and any radius [31].

Appendix C Monte-Carlo sampling in $\mathbb{R}^{d}$

C.1 Isotropic direction sampling

Uniform random sampling of directions on the $d$D sphere $S^{d-1}$ are generated by normalizing a vector of $d$ independent normally-distributed random numbers with a mean of [math], which we generate using Box-Muller transforms.

C.2 Unit uniform diffuse illumination sampling

For uniform diffuse illumination, we can invert the cumulative distribution function

[TABLE]

and solve for $\mu$ to find

[TABLE]

This gives a Monte-Carlo sampling procedure for uniform diffuse illumination of the half space by selecting incoming cosines $\mu$ from Eq.(234) where $0<\xi<1$ is a uniform random number.

Appendix D Connections to anisotropic scattering and nonclassical transport

In studying the variation of the transport kernel and related eigenspectra as a function of the dimensionality $d$ of the half space with isotropic scattering, we have come across transport kernels and related $H$ functions that exactly describe transport in 3D domains under anisotropic scattering or with correlation between the scattering centers in the medium. We briefly identify these correspondences here.

D.1 Anisotropic scattering

In the exact law of diffuse reflection for a 3D half space with linearly-anisotropic scattering ([1] p. 140) the solution involves an $H$ function with characteristic

[TABLE]

where the phase function $(1+b\cos\theta)$ has anisotropy parameter $-1<b<1$. This maps precisely to the $H$ function for isotropic scattering in 5D

[TABLE]

with a reduced absorption level

[TABLE]

Similarly, in the case of Rayleigh scattering in 3D, one of the three $H$ functions has a characteristic

[TABLE]

which corresponds to isotropic scattering in 7D with single-scattering albedo $c_{7}=1/10$. Results such as Eqs.(68) and (69) may offer previously unconsidered ways to compute transport solutions for these 3D problems.

D.2 Nonclassical transport

In the case of a random medium with spatially-correlated scatterers [97, 98], the exponential distribution of path-lengths between collisions in Eq.(24) becomes a general, non-negative, normalized distribution $p_{c}(s)$ and the plane-parallel generalized Schwarzschild-Milne kernel for the collision rate density becomes [44]

[TABLE]

In infinite spherical geometries it was shown [43] that the diffusion asymptotics of isotropic scattering in a classical medium of dimensionality $d$ can be exhibited by a non-exponential transport process in dimensionality $d^{\prime}$ by appropriate choice of $p_{c}(s)$. The same procedure can translate plane-parallel collision kernels $K(x)$ from one dimensionality to another. Of particular interest are the Picard and MacDonald kernels, whose related $H$ functions have analytic expressions. Either kernel can be exhibited in $\mathbb{R}^{d}$ under isotropic scattering by selecting the free-paths between collisions from a family of distributions involving Bessel-$K$ functions,

[TABLE]

whose Fourier-transformed propagator in $\mathbb{R}^{d}$ is a diffusion mode to a power $p$ [44],

[TABLE]

The Picard kernel in Eq.(22) arises when $p=1$, which happens in Flatland [43] via

[TABLE]

and in 3D via [43]

[TABLE]

The MacDonald kernel in Eq.(23) arises when $p=1/2$, which happens in 3D via

[TABLE]

The $H$ function for $p=2$ can also be written analytically by combining the 1D $H$ Function and Eq.(226) to yield

[TABLE]

While the same kernel can describe two different scattering processes in half spaces of differing dimensionalities, the solutions of the respective albedo problems are different because the rate density of initial collisions $C_{0}(z)$ is non-exponential in the case of correlation, and the relationship between the scalar collision rate and emerging distribution is less directly expressed via Laplace transforms. Complete details of solving half space albedo problems with correlation will be considered in a future paper [99].