PHYSICAL REVIEW E
110
, 025306 (2024)
Boundary-layer structures arising in linear transport theory
E. L. Gaggioli
,
1
,
2
,
3
,
*
Laura C. Estrada
,
1
,
2
,
†
and Oscar P. Bruno
3
,
‡
1
Universidad de Buenos Aires
, Facultad de Ciencias Exactas y Naturales, Departamento de Física. Buenos Aires 1428, Argentina
2
CONICET–Universidad de Buenos Aires,
Instituto de Física de Buenos Aires (IFIBA). Buenos Aires
1428, Argentina
3
Department of Computing and Mathematical Sciences,
Caltech
, Pasadena, California 91125, USA
(Received 26 March 2024; accepted 17 July 2024; published 28 August 2024)
We consider boundary-layer structures that arise in connection with the transport of neutral particles (e.g.,
photons or neutrons) through a participating medium. Such boundary-layer structures were previously identified
by the authors in certain particular cases [
Phys. Rev. E
104
, L032801 (2021)
]. Extending the previous work to
anisotropic scattering and general Fresnel boundary conditions, this contribution presents computational algo-
rithms which (1) resolve the aforementioned layers as well as previously unreported boundary layers associated
with Fresnel boundary transmission and reflection, and (2) yield accurate simulations at fixed computational cost
for transport under phase functions with arbitrarily strong anisotropy. The present paper additionally includes (3)
Mathematical proofs which justify the numerical methods proposed for resolution of boundary-layer structures.
The impact of the new theory on algorithmic performance is demonstrated through a series of 1D computational
benchmarks that emulate typical photon- and neutron-transport applications such as, e.g., optical tomography,
and nuclear reactor analysis and design. Experimental results for transmission of photons through turbid media
are presented, exhibiting close agreement between simulated and experimental data. As illustrated by means
of a variety of numerical results, the proposed boundary-layer-based approach tackles transport problems with
unprecedented accuracy and efficiency.
DOI:
10.1103/PhysRevE.110.025306
I. INTRODUCTION
We are concerned with the problem of transport of neutral
particles such as neutrons [
1
,
2
] and photons [
3
,
4
], wherein
interactions between particles are negligible and only interac-
tions between the particles and a background medium need to
be taken into account. The transport processes considered, in
which the background medium may absorb, emit, and
/
or scat-
ter neutral particles, is governed by the linear single-species
version of the Boltzmann equation, namely, the linear trans-
port equation—which is alternatively known as the radiative
transfer equation (RTE) and the neutron-transport equation in
the contexts of photon and neutron transport [
5
,
6
], respec-
tively. The linear transport equation presents the challenge of
high dimensionality: in full three-dimensional (3D) space the
solution (particle density) depends on three spatial variables
as well as two angular propagation variables and time, giving
rise, in all, to a computationally demanding six-dimensional
problem. Important applications include radiative heat trans-
port [
7
,
8
], gas dynamics [
2
], radiation transport in stellar
and planetary atmospheres [
3
,
9
,
10
], cancer diagnosis [
11
–
13
],
radiation therapy dose planning [
14
], finger joint arthritis di-
agnosis [
15
,
16
], studies of the human brain function [
17
,
18
],
optical and fluorescence tomography [
19
–
21
], and neu-
tron transport for nuclear reactor design [
22
–
24
], among
others [
25
–
27
].
*
Contact author: egaggioli@df.uba.ar
†
Contact author: lestrada@df.uba.ar
‡
Contact author: obruno@caltech.edu
This paper presents numerical algorithms that enable ef-
fective treatment of certain boundary-layer structures arising
in the solutions of the linear transport equation which were
previously identified in Ref. [
28
] under the assumptions of
vacuum boundary conditions and isotropic scattering. As
noted in Ref. [
28
], if left unaddressed, the presence of un-
bounded derivatives in the boundary-layer regions poses a
significant challenge for the numerical solution of transport
equations. Expanding on the previous work, which, in partic-
ular, provides a means for accurate and inexpensive numerical
treatment of boundary layers, the methods introduced in the
present paper enable treatment of problems under strongly
anisotropic scattering (on the basis of a certain multiresolution
approach for the evaluation of the scattering integral), and they
can be applied under general Fresnel boundary conditions
(which, as shown in Sec.
III A
and throughout this paper,
themselves give rise to a different type of boundary-layer
structure). The proposed methods for regularization of both
types of boundary layers are based on the use of certain
regularizing changes of variables under which unbounded
derivatives are eliminated. These methods, which are pre-
sented here in a spatially one-dimensional (1D) context, are
directly applicable to fully 3D configurations. This paper
additionally includes (1) a mathematical analysis of the reg-
ularization process and associated algorithms, providing a
firm mathematical basis for the proposed methods and (2)
comparisons with experimental results for turbid media ex-
hibiting close agreement between experiments and numerical
simulations.
The linear-transport boundary layers identified in Ref. [
28
],
whose theoretical description and numerical resolution via
2470-0045/2024/110(2)/025306(23)
025306-1
©2024 American Physical Society
GAGGIOLI, ESTRADA, AND BRUNO
PHYSICAL REVIEW E
110
, 025306 (2024)
changes of variables are detailed and generalized in the
present paper, are inherent in the transport equation solu-
tions near boundaries, and by extension, to the solutions of
Boltzmann-type equations. These boundary layers are unre-
lated to the ones widely discussed in the literature over a
period of various decades [
29
–
33
], which concern the asymp-
totic diffusion limit of the transport equation, that is, equations
of the form of Eqs. (
1
)or(
3
) in the limit of large scatter-
ing coefficient
μ
s
proportional to
ε
−
1
and small absorption
μ
a
=
μ
t
−
μ
s
, proportional to
ε
, where
ε
is a small parame-
ter. The boundary layers that arise in the resulting diffusion
approximation are independent of the angular variable
ξ
.In
contrast, the boundary layers considered in this paper exist for
arbitrary scattering and absorption coefficients, and they are
only observed for small values of
ξ
(directions nearly parallel
to the boundary) and for directions near the critical direction
ξ
=
ξ
0
c
of total internal reflection.
The remainder of this paper is organized as follows.
Section
II
briefly sets up necessary notations and conven-
tions. Section
III
then develops the boundary-layer theory
for the linear transport equation and it presents mathematical
proofs showing that all associated boundary-layer structures
can be completely resolved by means of certain types of
changes of variables in the spatial and angular coordinates.
Section
IV
present numerical techniques based on the the-
oretical background provided by the previous sections as
well as a multiresolution algorithm tailored to handle effec-
tively the scattering integral under highly anisotropic phase
functions. Numerical and experimental results are then pre-
sented in Sec.
V
, highlighting the efficiency and accuracy
of the proposed algorithms, and including applications to
configurations involving collimated sources (such as those
necessary for modeling laser beams), as well as an ex-
perimental illustration demonstrating excellent agreement
between theory and experiment. Finally, Sec.
VI
presents
a few concluding remarks and it outlines future research
directions emerging from the present work, including, in
particular, the introduction of enhanced imaging techniques
based on joint inversion of fluorescence-microscopy and dif-
fuse optical-tomography data, as suggested in Ref. [
34
], with
the goal of overcoming some of the inherent limitations of the
fluorescence-microscopy technique in high-scattering media.
II. BACKGROUND
For simplicity and definiteness this paper focuses on 1D
configurations (Fig.
1
) governed by the scalar radiative trans-
fer equation. The radiative transfer model is applicable to
propagation of monoenergetic particles, including neutron
transport [
5
] as well as photon transport under unpolar-
ized radiation as is often found in the context of light
propagation in highly scattering media such as biological
tissue [
15
,
16
,
21
,
32
]; discussions concerning the possible exis-
tence of polarization effects can be found in Refs. [
35
,
36
]. The
more general vector radiative transfer equation that accounts
for multiple polarization states [
4
,
37
], however, can be treated
similarly (cf. Remark
1
). Both time dependent and steady state
configurations are considered in the results and discussion
section, but for definiteness, the boundary-layer analysis is
presented in the steady state context only. (A corresponding
FIG. 1. 1D finite “slab” geometry (from Ref. [
28
]):
ξ
=
cos(
θ
).
analysis for the time-dependent case, which can be treated
by means of similar techniques, is left for future work.) The
analysis given below can also be extended to the general
multidimensional case, for which the 1D configuration might
be considered as a limiting case near the domain boundary.
The 1D linear transport problem
ξ
∂
∂
x
u
(
x
,ξ
)
+
μ
t
(
x
)
u
(
x
,ξ
)
=
μ
s
(
x
)
(
x
,ξ
)
+
q
(
x
,ξ
)
,
u
(0
,ξ
)
=
0
(
ξ
)
+
R
0
(
ξ
)
u
(0
,ξ
R
)
,
for
ξ>
0
,
u
(1
,ξ
)
=
1
(
ξ
)
+
R
1
(
ξ
)
u
(1
,ξ
R
)
,
for
ξ<
0
,
(1)
where
(
x
,ξ
)
=
1
−
1
p
(
ξ,ξ
)
u
(
x
,ξ
)
d
ξ
,
(2)
describes the dynamics of the angular flux
u
(
x
,ξ
), where
x
,
ξ
=
cos(
θ
) and
ξ
=
cos(
θ
) denote the spatial variable and
the cosines of the relevant propagation angles, as depicted
in Fig.
1
, respectively. Here, calling
μ
a
(
x
) and
μ
s
(
x
)the
absorption and scattering coefficients,
μ
t
(
x
)
=
μ
s
(
x
)
+
μ
a
(
x
)
denotes the total transport coefficient. The phase function
p
(
ξ,ξ
) models the density of probability that particles inci-
dent at
x
in directions between
ξ
and
ξ
+
d
ξ
emerge from
x
in the direction
ξ
after a scattering event. The “volumetric
source”
q
models the emission of particles, such as photons or
neutrons, within the medium, as a result of, e.g., fluorescence
or nuclear fission. The terms
i
(
ξ
)(
i
=
0
,
1) can be used to
model a source of photons at the domain boundaries, as may
be given by a laser beam that injects radiation for optical
tomography applications, and finally,
R
0
,
1
denote the Fresnel
coefficients which, as discussed below in this section, them-
selves give rise to boundary-layer structures whose adequate
resolution impacts significantly on the accuracy of the numer-
ical solution of the transport equation.
The corresponding time-dependent problem,
1
c
∂
∂
t
+
ξ
∂
∂
x
+
μ
t
(
x
)
u
(
x
,ξ,
t
)
=
μ
s
(
x
)
(
x
,ξ,
t
)
+
q
(
x
,ξ,
t
)
,
u
(
x
,ξ,
t
=
0)
=
0
,
025306-2
BOUNDARY-LAYER STRUCTURES ARISING IN LINEAR ...
PHYSICAL REVIEW E
110
, 025306 (2024)
u
(0
,ξ,
t
)
=
0
(
ξ,
t
)
+
R
0
(
ξ
)
u
(0
,ξ
R
,
t
)
,
for
ξ>
0
,
u
(1
,ξ,
t
)
=
1
(
ξ,
t
)
+
R
1
(
ξ
)
u
(1
,ξ
R
,
t
)
,
for
ξ<
0
,
(3)
where
(
x
,ξ,
t
)
=
1
−
1
p
(
ξ,ξ
)
u
(
x
,ξ
,
t
)
d
ξ
(4)
and where
c
denotes the speed of the particles in the
participating media (which is taken to equal
c
=
1 through-
out this paper), incorporates a time derivative in the
differential equation as well as the prescription of the values
of the solution at the initial time and time-dependent bound-
ary conditions, in addition to various elements present in
Eq. (
1
).
The previous contribution [
28
] announced results concern-
ing boundary-layer structures for the transport equations in
the presence of so-called vacuum boundary conditions—that
is,
u
(
x
,ξ
)
=
0at
x
=
0 and
x
=
1 for the relevant ranges
of incoming directions, namely 0
<ξ
1for
x
=
0 and
−
1
ξ<
0for
x
=
1. [Clearly, the vacuum boundary con-
ditions coincide with Fresnel boundary conditions, displayed
in Eqs. (
1
) and (
3
), with
R
0
,
1
=
0 and
0
,
1
=
0.] Photon
transport theories (as, e.g., required for applications in optical
tomography) generally require the use of the Fresnel bound-
ary conditions with nonzero Fresnel coefficients. The Fresnel
coefficient, which depends on the indexes of refraction
n
and
n
s
of the participating medium and its surroundings (both
of which are assumed throughout this paper to be spatially
constant for simplicity), quantifies the fraction of the radiation
arriving at the boundary in a given direction, that is reflected in
the corresponding specular direction. For the 1D slab geome-
try considered in this paper, radiation incident on the boundary
in the direction
ξ
=
cos(
θ
) is specularly reflected into direc-
tion
ξ
R
=−
ξ
. Calling
α
0
=
π
−
θ
the incidence angle with
respect to the normal ˆ
ν
0
=
(1
,
0
,
0) at
x
=
0 and
α
1
=
θ
the
incidence angle with respect to the normal ˆ
ν
1
=
(
−
1
,
0
,
0) at
x
=
1, and letting
α
0
,
1
denote either
α
0
or
α
1
depending on
whether the normal at ˆ
ν
0
at
x
=
0or ˆ
ν
1
at
x
=
1 is considered,
the corresponding Fresnel coefficient
R
0
,
1
(similarly equal to
either
R
0
or
R
1
)isgivenby
R
0
,
1
(
ξ
)
=
⎧
⎪
⎪
⎪
⎪
⎨
⎪
⎪
⎪
⎪
⎩
n
−
n
s
n
+
n
s
2
if
α
0
,
1
=
0
,
1
2
sin
2
(
α
0
,
1
t
−
α
0
,
1
)
sin
2
(
α
0
,
1
t
+
α
0
,
1
)
+
tan
2
(
α
0
,
1
t
−
α
0
,
1
)
tan
2
(
α
0
,
1
t
+
α
0
,
1
)
if 0
<α
0
,
1
<α
0
,
1
c
,
1if
α
0
,
1
α
0
,
1
c
,
(5)
where the critical angle is given by sin(
α
0
,
1
c
)
=
n
s
n
, with asso-
ciated critical directions
ξ
0
c
=
cos
α
0
c
and
ξ
1
c
=
cos
α
1
c
(6)
(beyond which, for
n
>
n
s
, total internal reflection occurs
at the boundaries
x
=
0 and
x
=
1, respectively). Note that,
under the present assumptions of a constant index of refraction
we have
ξ
1
c
=−
ξ
0
c
.
(7)
The transmission angles
α
0
,
1
t
are obtained from Snell’s law of
refraction sin(
α
0
,
1
t
)
=
n
n
s
sin(
α
0
,
1
).
On the basis of the Fresnel coefficient we additionally
define the measurement operator
J
+
[
u
](
t
)
=
1
0
[1
−
R
1
(
ξ
)]
u
(1
,ξ,
t
)
d
ξ,
(8)
which quantifies the time-resolved outgoing flux of photons at
the boundary domain at
x
=
1.
III. BOUNDARY-LAYER THEORY FOR THE LINEAR
TRANSPORT EQUATION
A. Integral equation formulation: Fresnel and incidence
boundary layers
An integral equation equivalent to Eq. (
1
) can readily be
obtained by utilizing an integrating factor and substituting
the
x
=
0 and
x
=
1 boundary conditions (
1
) in the resulting
equation. For
ξ>
0 we thus obtain
u
(
x
,ξ
)
=
(
0
(
ξ
)
+
R
0
(
ξ
)
u
(0
,ξ
R
))
e
−
x
0
μ
t
(
y
)
dy
/ξ
+
e
−
x
0
μ
t
(
y
)
dy
/ξ
ξ
x
0
e
y
0
μ
t
(
z
)
dz
/ξ
q
(
y
,ξ
)
dy
+
e
−
x
0
μ
t
(
y
)
dy
/ξ
ξ
x
0
e
y
0
μ
t
(
z
)
dz
/ξ
μ
s
(
y
)
1
−
1
p
(
ξ,ξ
)
u
(
y
,ξ
)
d
ξ
dy
,ξ>
0
,
(9)
with a similar expression for
ξ<
0. In view of the ideas underlying boundary-layer theory [
38
], the particular case of Eq. (
9
)in
which the quantities
μ
s
,
μ
t
, and
q
are spatially constant captures the leading order boundary-layer behavior [
28
]ofthesolution
u
of the variable-coefficients problem—since, for a lowest-order asymptotic approximation
u
0
(
x
,ξ
), where
u
(
x
,ξ
)
∼
u
0
(
x
,ξ
)
as (
x
,ξ
)
→
(0
+
,
0
+
), the coefficients can indeed be assumed to be constant in a small neighborhood of the boundary points.
Assuming such constant coefficients and using Eqs. (
2
) and (
9
) leads to the near-boundary approximation
u
(
x
,ξ
)
∼
u
0
(
x
,ξ
)
=
(
0
(
ξ
)
+
R
0
(
ξ
)
u
(0
,ξ
R
))
e
−
μ
t
(0)
x
/ξ
+
e
−
μ
t
(0)
x
/ξ
ξ
x
0
e
μ
t
(0)
y
/ξ
q
(0
,ξ
)
dy
+
e
−
μ
t
(0)
x
/ξ
ξ
x
0
e
μ
t
(0)
y
/ξ
μ
s
(0)
(
y
,ξ
)
dy
,
as (
x
,ξ
)
→
(0
+
,
0
+
)
,
(10)
025306-3
GAGGIOLI, ESTRADA, AND BRUNO
PHYSICAL REVIEW E
110
, 025306 (2024)
with a similar result for (
x
,ξ
)
→
(1
−
,
0
−
). The exponen-
tial factor exp (
−
μ
t
(0)
x
/ξ
) represent fast transitions in the
density of particles traveling near the boundary point
x
=
0in
directions nearly parallel to the boundary (for which
ξ
is close
to zero). In particular, these exponential boundary-layer terms
entail unbounded derivatives in both, the
x
and
ξ
variables
as (
x
,ξ
)
→
(0
+
,
0
+
). Clearly, such boundary layers present a
challenge from a computational standpoint—since the numer-
ical solution of Eq. (
1
) requires, in particular,
x
differentiation
and integration respect to
ξ
across such structures.
The boundary-layer terms
0
(
ξ
)
e
−
x
0
μ
t
(
y
)
dy
/ξ
in Eq. (
9
)
and
0
(
ξ
)
e
−
μ
t
(0)
x
/ξ
in Eq. (
10
) represent the uncollided and
unabsorbed remainder of the incident flux entering through
the domain boundary
x
=
0; these terms are therefore called
“incidence boundary layers” (IBL). Similarly, the third and
fourth boundary-layer terms in Eqs. (
9
) and (
10
) incorporate
exponential increments and attenuations—arising from the
emission of particles resulting from the volumetric source
q
in the third terms, and from particle scattering
in the fourth
terms, and they are therefore called “volumetric source bound-
ary layers” (VSBL). The terms
R
0
(
ξ
)
u
(0
,ξ
R
)
e
−
x
0
μ
t
(
y
)
dy
/ξ
and
R
0
(
ξ
)
u
(0
,ξ
R
)
e
−
μ
t
(0)
x
/ξ
in these equations, finally, rep-
resent the uncollided and unabsorbed remainders of Fresnel
reflections of interior fields at the boundary
x
=
0, and they
incorporate boundary layers of a different kind, which we call
“Fresnel boundary layers” (FBL). Unlike the IBL and VSBL,
the Fresnel boundary layers around, e.g.,
x
=
0 arises from a
singularity of the Fresnel reflection coefficient
R
0
(
ξ
)atthe
boundary
ξ
0
c
between the angular regimes of partial and total
internal reflection, as discussed in what follows. Notably, un-
like the IBL and VSBL, the Fresnel boundary layers induce a
singular behavior throughout the spatial domain, and, strictly
speaking, they are not confined to the domain boundary. How-
ever, owing to the exponential factor that accompanies the
Fresnel coefficients, such singularities decay exponentially
fast with the distance to the boundary, they are therefore only
observed in a neighborhood of the boundary (see, e.g., Fig.
7
),
and, in that sense, they are boundary layers as well.
As mentioned above, the FBL arise as
ξ
approaches the to-
tal internal-reflection direction
ξ
0
c
from the right (respectively,
ξ
1
c
from the left). At such points all derivatives of
R
0
,
1
(
ξ
)
become infinite, which gives rise to fast transitions in the solu-
tion
u
around the points (
x
,ξ
)
=
(0
,ξ
0
c
) and (
x
,ξ
)
=
(1
,ξ
1
c
),
as illustrated in the aforementioned Fig.
7
. To see this for
R
0
(
ξ
) (the case
R
1
(
ξ
) is analogous) we note a singularity that
arises in this function from the term
α
0
t
=
arcsin
n
n
s
1
−
ξ
2
(11)
that appears repeatedly in Eq. (
5
). Expanding
n
n
s
√
1
−
ξ
2
around
ξ
=
ξ
0
c
(for which
n
n
s
1
−
(
ξ
0
c
)
2
=
1) we ob-
tain
n
n
s
√
1
−
ξ
2
∼
1
+
n
2
n
2
s
ξ
0
c
(
ξ
−
ξ
0
c
) which, together with
Eq. (
11
) and arcsin(
y
)
∼
π
2
−
√
2
√
1
−
y
for
y
1 yields
α
0
t
∼
π
2
−
n
n
s
2
ξ
0
c
ξ
−
ξ
0
c
,
(12)
or, more precisely,
α
0
t
=
F
ξ
−
ξ
0
c
,ξ>ξ
0
c
,
(13)
where
F
=
F
(
z
) is a smooth function of
z
around
z
=
0. These
expressions encapsulate the singular character
R
0
(
ξ
)
=
1
−
a
ξ
−
ξ
0
c
+
O
ξ
−
ξ
0
c
,
(14)
for some real coefficient
a
>
0, or, more generally,
R
0
,
1
(
ξ
)
=
S
0
,
1
ξ
−
ξ
0
,
1
c
1
2
,
(15)
of the complete Fresnel coefficient and its unbounded deriva-
tives for
ξ
ξ
0
c
and
ξ
ξ
1
c
, where
S
0
,
1
=
S
0
,
1
(
z
) is a smooth
function of
z
around
z
=
0. As indicated above, the exponen-
tially decaying factors that accompany the Fresnel coefficients
make the fast FBL transitions observable only in a small
region near the boundary points, as befits a boundary-layer
structure.
Remark 1.
As suggested in Sec.
II
, with exception of minor
variations required to account for the simultaneous presence
of multiple polarization states, the theoretical and computa-
tional methods presented in this paper are applicable in the
context of vector radiative transfer equations for polarized
light. In particular, since Snell’s law (
11
) is valid for arbitrary
polarization states, the singular character of the Fresnel re-
flection matrix [
37
, Eq. (10)] is once again, as in Eqs. (
13
)
and (
15
), given by smooth functions of
|
ξ
−
ξ
0
,
1
c
|
1
2
for
ξ
ξ
0
c
and
ξ
ξ
1
c
. It follows that the change of variables (
18
), which
was designed to regularize functions containing such types
of singularities, also produces the required regularization for
the Fresnel reflection matrices and associated vector-transport
solutions that arise under polarized radiation.
B. Boundary-layer resolving changes of variables
As discussed in what follows, certain angular and spatial
changes of variables can be utilized to fully resolve nu-
merically the near-singular IBL, VSBL, and FBL structures
described in the previous section; such a result was previously
announced for the VSBL structures in Ref. [
28
]. The theo-
retical results presented in Secs.
III C
and
III D
, which are
illustrated in Fig.
6
, show that, indeed, such transformations
“regularize” the problem: in the new variables the spatial
and angular derivatives remain uniformly bounded throughout
the spatial and angular domains, including points in space
arbitrarily close to the domain boundary and directions arbi-
trarily close to tangential to the boundary. As a result, use of
the new variables gives rise to rapidly convergent numerical
algorithms that, on the basis of relatively coarse discretiza-
tions, evaluate accurately the transport solution throughout the
angular and spatial domains, and, in particular, for arbitrarily
small values of
ξ
and for values of
x
arbitrarily close to the
boundaries
x
=
0 and
x
=
1.
As a simple example, let us consider the function
w
(
ξ
)
=
√
ξ
which, like the Fresnel coefficient, has infinite derivatives
at a point—in this case, the point
ξ
=
0. Clearly, use of the
change of variables
ξ
=
r
n
with an integer
n
>
1inanintegra-
tion problem results in the smoother integrand
n
w
(
r
n
)
r
n
−
1
=
nr
3
n
/
2
−
1
—which, of course, is infinitely differentiable for
n
025306-4
BOUNDARY-LAYER STRUCTURES ARISING IN LINEAR ...
PHYSICAL REVIEW E
110
, 025306 (2024)
even, and, in particular, for
n
=
2. Note that an equispaced
mesh in the
r
variable corresponds to a graded mesh in the
ξ
variable: the mesh grading transforms a curve with an infinite
slope in the
ξ
variable into a curve with a finite slope in
the
r
variable—a property that will be exploited not only in
connection with the Fresnel boundary layers, but also, with
values
n
>
2, to resolve the challenging nearly singular
ξ
dependence of the aforementioned IBL and VSBL.
As shown in Secs.
III C
and
III D
, changes of variables can
similarly be utilized to resolve the unbounded
x
derivatives
that occur in a family of functions such as
w
(
x
,ξ
)
=
1
−
e
−
x
ξ
.
Whereas none of these functions has an infinite derivative,
the derivatives of the family of functions are collectively un-
bounded, as, e.g., the
x
-derivatives of
w
ξ
(
x
) (e.g., the first
derivative
w
ξ
(
x
)
=−
1
ξ
e
−
x
ξ
) tend to infinity as
x
and
ξ
ap-
proach zero with fixed values of
x
/ξ
; similar unbounded
behavior arises for other derivatives of this function with
respect to
x
and
/
or
ξ
. As in the case of the square-root func-
tion considered above, changes of variables can be utilized
to resolve the
x
near-singularity of the complete family. As a
simple preliminary illustration we note that, under the change
of variables
x
=
x
(
v
)
=
e
v
e
v
+
1
,
(16)
the family
w
ξ
(
x
(
v
))
=
1
−
e
−
x
(
v
)
ξ
has uniformly bounded
derivatives with respect to
v
. Indeed, for the first derivative,
for example, we obtain
d
d
v
w
ξ
(
x
(
v
))
=−
(1
−
x
(
v
))
x
(
v
)
ξ
e
−
x
(
v
)
ξ
,
which is uniformly bounded, for all
v
and
ξ
, on account
of the boundedness of the functions (1
−
x
(
v
)) and
Xe
−
X
with
X
=
x
(
v
)
ξ
. A similar calculation shows that all of the
v
derivatives of
w
ξ
(
x
(
v
)) are uniformly bounded. Use of an
additional order-
n
algebraic change of variables
ξ
=
(
r
),
finally, ensures that the derivatives with respect to
v
of all
orders, and the derivatives with respect to
r
up to order (
n
−
1)
are uniformly bounded—for all
v
and
r
.
As mentioned above, the change of variables we use
to treat the Fresnel boundary-layer singularity expressed in
Eq. (
14
) is related to but different from the simple algebraic
transformation
ξ
=
r
n
discussed previously. Like the simple
transformation, the alternative algebraic transformation we
use induces a mesh grading via a power-
n
expression which
maintains an adequate number of discretization points away
from the finely graded region next to the singular point (which
corresponds to
ξ
=
r
=
0inthesimple
√
ξ
example consid-
ered above). In detail, instead of
ξ
=
r
n
, in this paper we
use rescaled versions [via linear functions
s
=
s
(
r
)] of the
algebraic Martensen-Kussmaul (MK) change of variables [
39
]
given by
h
N
(
s
)
=
2
π
[
v
(
s
)]
N
[
v
(
s
)]
N
+
[
v
(2
π
−
s
)]
N
,
0
s
2
π,
v
(
s
)
=
1
N
−
1
2
π
−
s
π
3
+
1
N
s
−
π
π
+
1
2
.
Roughly speaking, this transformation accumulates one half
of the grid points toward the endpoints 0 and 2
π
, while the
other half is distributed fairly uniformly within the interior of
the interval [0
,
2
π
]. In the context of angular boundary-layer
resolution we utilize this transformation with two different
values of
N
, in conjunction with the linear rescaling
s
=
s
(
r
)
=
⎧
⎨
⎩
π
r
ξ
0
c
for
0
<
r
ξ
0
c
,
π
r
−
ξ
0
c
1
−
ξ
0
c
for
ξ
0
c
<
r
1
.
(17)
In detail, in view of the relation (
7
), we introduce the angular
change of variables
ξ
=
(
r
)
=
ψ
(
r
)for0
r
1
,
−
ψ
(
−
r
)for
−
1
r
<
0
,
(18)
where, for a given value
n
(taken to equal
n
=
3
.
2 throughout
this paper),
ψ
(
r
)
=
h
n
(
s
(
r
))
π
ξ
0
c
,
0
<
r
ξ
0
c
,
ξ
0
c
+
1
−
ξ
0
c
π
h
2
(
s
(
r
))
,ξ
0
c
<
r
1
.
(19)
Note that this change of variables results in a numerical grid
in the
ξ
variable that is refined near
ξ
=
0 and
ξ
=±
ξ
0
c
;as
shown in the following section, such graded refinements pro-
vide the necessary resolution of boundary layers. As discussed
in the following section, the function
h
2
, that is to say,
h
N
with
N
=
2, which is used in the case
ξ
0
c
<
r
1 of definition (
19
),
exactly cancels the square-root singularity (
15
) in the Fresnel
coefficient. The function
h
N
with
N
=
n
used for the case
0
<
r
ξ
0
c
, in turn, is utilized to smoothen the singularity of
the solution
u
near
ξ
=
0. A different use of the MK function
h
N
(
s
) is made in Sec.
VA
, with a different linear rescaling, to
adequately discretize a collimated boundary source.
The formal analyses presented in the following section for
both the angular and spatial regularization processes is based
on use of the convergent Neumann series representation [
1
,
40
]
u
(
x
,ξ
)
=
∞
m
=
0
u
m
(
x
,ξ
)
,
(20)
of the solution
u
, where
u
m
(
x
,ξ
)
=
K
m
[
g
+
L
[
q
]](
x
,ξ
)isde-
fined in Appendix
A
as the result of the action of the
m
th
power
K
m
of a certain operator
K
on a “boundary-condition”
function
g
as well as the result
L
[
q
] of the action of an operator
L
on the source function
q
in Eq. (
1
). Physically,
u
m
(
x
,ξ
)
represents the density of photons at point
x
traveling with
direction given by
ξ
that have undergone a number
m
of
collision and scattering events.
As is often the case in practice [and as established in Ap-
pendix
B
for the 1D Henyey-Greenstein phase function given
in Eqs. (
44
) and (
45
)], the phase function
p
=
p
(
ξ,ξ
) and its
derivatives are assumed to be bounded: for each integer
j
0
there exists a constant
C
j
>
0 such that
∂
j
∂ξ
j
p
(
ξ,ξ
)
<
C
j
for
−
1
ξ,ξ
1
.
(21)
It follows that provided
u
=
u
(
x
,ξ
) is a bounded function of
x
and
ξ
(as it generally may be expected on physical grounds),
for each integer
j
0 there exists a constant
D
j
>
0 such
025306-5
GAGGIOLI, ESTRADA, AND BRUNO
PHYSICAL REVIEW E
110
, 025306 (2024)
that for
−
1
ξ
1 and 0
x
1 we have for the collision
integral (
2
)
∂
j
∂ξ
j
(
x
,ξ
)
<
D
j
,
(22)
a fact that will exploited in the theoretical study of angular
boundary layers presented in Sec.
III C
. Similarly, defining
m
(
x
,ξ
)
=
1
−
1
p
(
ξ,ξ
)
u
m
(
x
,ξ
)
d
ξ
,
(23)
we may generally assume
∂
j
∂ξ
j
m
(
x
,ξ
)
<
E
m
j
,
(24)
for certain constants
E
m
j
. Naturally, our formal study of so-
lution smoothness under the changes of variables introduced
above assumes that the boundary sources
0
(
ξ
) and
1
(
ξ
) and
the interior source
q
(
x
,ξ
) are smooth functions, with bounded
derivatives of all orders, with respect to
x
and
ξ
. For nota-
tional simplicity we present proofs for the case of spatially
constant absorption and scattering coefficients,
μ
a
(
x
)
=
μ
a
,
μ
s
(
x
)
=
μ
s
, and
μ
t
(
x
)
=
μ
t
=
μ
a
+
μ
s
.
C. Boundary-layer regularization I: Angular regularization
This section shows that the change of variables (
19
) regu-
larizes the exponential and algebraic boundary layers [namely,
the quantities (
A2
) and (
A3
) that, in view of Eqs. (
A7
)
through (
A11
), are a part of the terms (
A11
) of the Neumann
expansion (
A9
) of the solution, see also Fig.
7
] that occur in
the integrand
p
(
ξ,ξ
)
u
(
x
,ξ
) of the collision integral (
2
), with
respect to the integration variable
ξ
for
ξ
around
ξ
=
0.
(Without loss of generality we restrict our proof to the case
ξ>
0; the case
ξ<
0 follows analogously.) More precisely,
the results in this section show that the collision integrand that
results upon the change of integration variables
ξ
=
ψ
(
r
),
namely, [
p
(
ξ,ψ
(
r
))
u
(
x
,ψ
(
r
))
d
ψ
dr
(
r
)], has bounded deriva-
tives
∂
j
r
with respect to
r
, of orders
j
(
n
−
1) for 0
r
1
and for all
x
,0
x
1—and, thus, integration of this inte-
grand on the basis of, e.g., the Gauss-Legendre rule gives rise
to high-order convergence, of orders consistent with with clas-
sical error estimates [
41
] for the Gauss-Legendre quadrature
rule: the
-point Gauss-Legendre quadrature error decreases
like 32
V
/
[15
π
(
n
−
1)(2
−
n
+
2)
n
−
1
]for
n
−
1
2
, provided
the derivatives of the integrated function up to order (
n
−
1)
are bounded by the constant
V
>
0. Since, in the present case
ξ
0 the function
u
(
x
,ξ
) has bounded derivatives outside a
half-neighborhood of
ξ
=
0 to the right of
ξ
=
0 and a half
neighborhood of
ξ
=
ξ
0
c
to the right of
ξ
=
ξ
0
c
(see Fig.
7
), we
show that under the composition
u
(
x
,ψ
(
r
)) the resulting inte-
grand has bounded derivatives of order 0
j
(
n
−
1), for
0
r
1, in both, half-neighborhoods of
r
=
0 and
r
=
ξ
0
c
to the right of
r
=
0 and
r
=
ξ
0
c
, respectively.
We establish first the derivative boundedness near
r
=
0
and to the right of this point. To do this, noting that the
change of variables (
19
) behaves asymptotically like the
n
th
power function,
ψ
(
r
)
∼
r
n
,as
r
→
0
+
, it suffices to show that
the change of variables
ξ
=
ζ
1
(
r
)
=
r
n
results in the desired
bounded integrand derivatives for the region 0
r
<ξ
0
c
and
up to and including
r
=
ξ
0
c
from the left. Under the present hy-
pothesis Eq. (
21
) that the phase function
p
(
ξ,ξ
) has bounded
derivatives, and since
ζ
1
(
r
)
=
nr
n
−
1
, it suffices to show that
∂
j
∂
r
j
(
u
(
x
,
r
n
)
r
n
−
1
) is bounded for 0
r
<ξ
0
c
and up to and
including
r
=
ξ
0
c
from the left, and for 0
j
(
n
−
1). We
do this by showing, by induction in
m
, that the same is true for
each term
u
m
in the Neumann series (
20
)of
u
. That is to say
that for some positive constant
A
j
, we show inductively that
∂
j
∂
r
j
(
u
m
(
x
,
r
n
)
r
n
−
1
)
<
A
j
,
0
j
n
−
1
,
(25)
and thus, the
j
th derivative of the integrand in the scattering
integral is bounded for such values of
r
and
j
.
Proof of Eq. (
25
)
. We first consider the case
m
=
0for
which, in view of Eqs. (
A1
), (
A2
), and (
A11
) in Appendix
A
,
∂
j
/∂
r
j
(
u
0
(
x
,
r
n
)
r
n
−
1
) equals
∂
j
∂
r
j
(
0
(
r
n
)
e
−
μ
t
x
/
r
n
r
n
−
1
+
L
[
q
](
x
,
r
n
)
r
n
−
1
)
,
(26)
and we consider, in turn, each one of the two terms on the
right-hand side of Eq. (
26
). In regard to the first right-hand
term, given that, by assumption, the boundary source term
0
(
ξ
) has bounded derivatives of all orders, it suffices to
show that
e
−
μ
t
x
/
r
n
r
n
−
1
has bounded derivatives for the relevant
orders of differentiation. Using Eqs. (
A12
) and (
A13
), and
calling
X
=
μ
t
x
/
r
n
we obtain the expression
∂
j
∂
r
j
(
e
−
μ
t
x
/
r
n
r
n
−
1
)
=
r
n
−
j
−
1
j
=
0
j
−
−
1
s
=
0
a
,
s
X
j
−
−
s
e
−
X
,
(27)
with real coefficients
a
,
s
. Since, for any real constant
M
0,
X
M
e
−
X
is bounded for all
X
0
,
(28)
and since all of the exponents (
n
−
j
−
1) on the right-hand
side of the equation are nonnegative for 0
j
(
n
−
1), it
follows that the left-hand quantity in Eq. (
27
) is bounded for
this range of values of
j
, as desired. For the second term in
Eq. (
26
), in turn, using Eq. (
A2
) in the case
ξ>
0 considered
presently, together with Eq. (
A16
), we obtain
∂
j
∂
r
j
(
L
[
q
](
x
,
r
n
)
r
n
−
1
)
=
r
n
−
j
−
1
j
k
=
0
j
=
1
j
−
k
s
=
1
α
j
,
k
,,
s
r
sn
×
x
0
∂
s
q
∂ξ
s
(
y
,
r
n
)
μ
t
(
y
−
x
)
r
n
e
−
μ
t
(
y
−
x
)
/
r
n
dy
,
(29)
for certain coefficients
α
j
,
k
,,
s
. Then, using Eq. (
28
) with
X
=
μ
t
(
y
−
x
)
/
r
n
, and since, by hypothesis, the source term
q
(
y
,ξ
) has bounded derivatives, it follows that the integrand
on the right-hand side of Eq. (
29
), and, thus, the complete
right-hand side, is also bounded for all 0
j
n
−
1 and for
0
r
ξ
0
c
.The
m
=
0 case of the inductive proof has thus
been concluded.
To complete the inductive proof in the present case (in the
r
region to the left of
r
=
ξ
0
c
), we assume that for a given
m
the
derivatives
∂
j
r
(
r
n
−
1
u
m
) are bounded for 0
j
(
n
−
1), and
we show that the same is true for the derivatives
∂
j
r
(
r
n
−
1
u
m
+
1
),
025306-6
BOUNDARY-LAYER STRUCTURES ARISING IN LINEAR ...
PHYSICAL REVIEW E
110
, 025306 (2024)
or, equivalently, in view of Eq. (
A7
) through (
A11
), that the
terms
∂
j
∂
r
j
(
R
0
(
r
n
)
u
m
(0
,
−
r
n
)
e
−
μ
t
x
/
r
n
r
n
−
1
)
and
∂
j
∂
r
j
(
μ
s
L
[
S
[
u
m
]](
x
,
r
n
)
r
n
−
1
)
(30)
are bounded for 0
j
(
n
−
1). By the induction hypothesis
and in view of Eqs. (
27
) and (
28
) we see that the first term
derivatives are bounded for 0
j
(
n
−
1). For the second
term, in turn, using once again Eq. (
A2
), in conjunction with
Eqs. (
A4
), (
A16
), and (
23
), we obtain
∂
j
∂
r
j
(
L
[
S
[
u
m
]](
x
,
r
n
)
r
n
−
1
)
=
r
n
−
j
−
1
j
k
=
0
j
=
1
j
−
k
s
=
0
α
j
,
k
,,
s
r
sn
×
x
0
∂
s
m
∂ξ
s
(
y
,
r
n
)
μ
t
(
y
−
x
)
r
n
e
−
μ
t
(
y
−
x
)
/
r
n
dy
.
In view of Eq. (
24
) together with the bound (
28
) with
X
=
μ
t
(
y
−
x
)
/
r
n
we see that the integrand on the right-hand side
is bounded, and, thus that, for 0
j
(
n
−
1), so is the
complete second term in Eq. (
30
).
Having established the integrand derivative boundedness
under the change of variables (
19
) in the region 0
r
<
ξ
0
c
and up to and including
r
=
ξ
0
c
from the left, we now
proceed to establish the corresponding boundedness in the
remaining region
ξ
0
c
<
r
1 up to and including
r
=
ξ
0
c
from the right—in which the change of variables (
19
) be-
haves asymptotically like
ψ
(
r
)
∼
ξ
0
c
+
b
(
r
−
ξ
0
c
)
2
as
r
→
ξ
0
c
from the right for a certain real constant
b
.Inviewofthis
asymptotic character, it suffices to show that under the change
of variables
ξ
=
ζ
2
(
r
)
=
ξ
0
c
+
b
(
r
−
ξ
0
c
)
2
,the
j
th derivative
∂
j
∂
r
j
(
u
(
x
,ζ
2
(
r
))
d
ζ
2
(
r
)
dr
) is bounded for
ξ
0
c
r
1, which we
establish, as in the previous case 0
r
<ξ
0
c
, by showing,
inductively, that the same is true for each Neumann-series
term
u
m
. In other words, we show that
∂
j
∂
r
j
u
m
(
x
,ζ
2
(
r
))
d
ζ
2
(
r
)
dr
(31)
is bounded for all nonnegative integers
m
and for
ξ
0
c
<
r
1
up to and including
r
=
ξ
0
c
from the right.
Once again, we begin with the case
m
=
0, for which
Eq. (
31
) equals
∂
j
∂
r
j
0
(
ζ
2
(
r
))
e
−
μ
t
x
/ζ
2
(
r
)
d
ζ
2
(
r
)
dr
+
∂
j
∂
r
j
L
[
q
](
x
,ζ
2
(
r
))
d
ζ
2
(
r
)
dr
,
both of whose terms are clearly bounded since
ξ>ξ
0
c
>
0. To
complete the inductive proof we assume the derivatives (
31
)
are bounded for a certain integer
m
, and we show that the
same is true for
m
+
1. To do this we note that, under the
ζ
2
change of variables, the
j
th derivatives of the (
m
+
1)th term
of the Neumann expansion equal the sum of the following two
terms:
∂
j
∂
r
j
R
0
(
ζ
2
(
r
))
u
m
(0
,
−
ζ
2
(
r
))
e
−
μ
t
x
/ζ
2
(
r
)
d
ζ
2
(
r
)
dr
and
∂
j
∂
r
j
μ
s
L
[
S
[
u
m
]](
x
,ζ
2
(
r
))
d
ζ
2
(
r
)
dr
.
But, by the induction hypothesis
u
m
(
x
,ζ
2
(
r
)) has bounded
derivatives with respect to
r
for
ξ
c
<
r
1 up to and includ-
ing
r
=
ξ
0
c
from the right, and, thus, to show that the same is
true for
u
m
+
1
(
x
,ζ
2
(
r
)) it suffices to show that
∂
j
∂
r
j
R
0
(
ζ
2
(
r
))
is itself bounded. But, from Eq. (
15
), for a certain smooth
function
S
we have that
R
0
(
ζ
2
(
r
))
=
S
√
b
r
−
ξ
0
c
for
ξ
0
c
r
1. It follows that
R
0
(
ζ
2
(
r
)) is smooth for
ξ
0
c
<
r
1 and up to and including
r
=
ξ
0
c
from the right and,
thus the
j
th derivatives of
u
m
+
1
are bounded in that region,
as desired. The inductive proof is now complete, establishing
that
∂
j
∂
r
j
p
(
ξ,ψ
(
r
))
u
(
x
,ψ
(
r
))
d
ψ
dr
(
r
)
is bounded for all 0
r
1, 0
x
1 and 0
j
(
n
−
1).
D. Boundary-layer regularization II: Spatial regularization
The discussion in this section shows that the change of
variables (
16
) regularizes the solution
u
(
x
,ξ
) in the both the
left and right boundary-layer regions, namely
μ
t
x
ξ
<ε
(0
<
ξ
1) and
μ
t
(
x
−
1)
ξ
<ε
(
−
1
ξ<
0), respectively—which,
e.g., for small values of
|
ξ
|
, are small regions near
x
=
0
and
x
=
1 within the physical domain 0
x
1. (Here 0
<
ε
1 is an arbitrary number.) In the present context the so-
lution
u
(
x
,ξ
) is said to be regularized in the sense that the
v
-derivatives of the composition
u
(
x
(
v
)
,ξ
) of arbitrary order
are bounded for all values of
v
for which
x
(
v
) is contained
within the left and right boundary-layer regions. This fact is
established in this section (using
ε
=
1 for definiteness) by
showing that each term
u
m
in the Neumann series (
A9
) has
this property, that is to say, for each nonnegative integer
j
there exists a positive constant
B
j
such that
∂
j
∂
v
j
u
m
(
x
(
v
)
,ξ
)
<
B
j
(32)
for all
v
and
ξ
, with
−∞
<
v
<
∞
and
|
ξ
|
1, such
that (
x
(
v
)
,ξ
) lies in either of the two boundary-layer re-
gions. As in Sec.
III C
, without loss of generality we
restrict our proof to the case
ξ>
0; the case
ξ<
0 follows
similarly.
Remark 2.
Our spatial regularization proof relies on the
fact that all spatial derivatives of the solution
u
with respect
to
x
are bounded
outside
the boundary-layer regions—i.e.,
for
μ
t
x
ξ
>
1, 0
<ξ
1, and for
μ
t
(
x
−
1)
ξ
>
1,
−
1
ξ<
0—as
might be expected from standard asymptotic boundary-layer
theory with asymptotic matching [
38
]. Existing theoretical
results in suitable functional spaces [
42
] only provide limited
insights in this regard. But, certainly, computational investi-
gations available in the literature [
30
,
43
], as well as our own
025306-7
GAGGIOLI, ESTRADA, AND BRUNO
PHYSICAL REVIEW E
110
, 025306 (2024)
high-resolution simulations such as those presented in Fig.
6
,
computationally demonstrate the needed derivative bounded-
ness outside the boundary-layer regions. A theoretical proof
of the away-from-boundary derivative boundedness, which
would certainly be valuable, and which could be based on
consideration of asymptotics of integrals, is beyond the scope
of this paper and is left for future work.
Proof of Eq. (
32
)
. The proof is presented in what fol-
lows for all nonnegative integers
m
and
j
and for spatial
and angular points
x
(
v
) and
ξ
in the left boundary-layer
region
|
X
|=|
μ
t
x
(
v
)
/ξ
|
<
1,
ξ>
0 (the proof in the right
boundary layer
|
μ
t
(
x
(
v
)
−
1)
/ξ
|
<
1,
ξ<
0 is analogous).
Proceeding by induction in
m
we thus first consider the
j
th
derivative in the case
m
=
0 which, in view of Eq. (
A10
), is
given by
0
(
ξ
)
∂
j
∂
v
j
e
−
μ
t
x
(
v
)
/ξ
+
∂
j
∂
v
j
L
[
q
](
x
(
v
)
,ξ
)
.
(33)
Using
X
=
μ
t
x
(
v
)
ξ
and Eq. (
A18
), the first term in Eq. (
33
)is
seen to equal
0
(
ξ
)
e
−
X
j
k
=
1
j
=
0
b
k
,
X
k
x
(
v
)
(34)
for some real coefficients
b
k
,
; clearly, in view of Eq. (
28
)this
term is bounded for all relevant values of (
v
,ξ
). Integrating
by parts
j
+
1 times and utilizing Eq. (
A20
), for the second
term we obtain
∂
j
∂
v
j
L
[
q
](
x
(
v
)
,ξ
)
=
j
+
1
=
1
(
−
1)
−
1
ξ
−
1
μ
t
∂
j
∂
v
j
∂
−
1
∂
x
−
1
q
(
x
(
v
)
,ξ
)
−
e
−
μ
t
x
(
v
)
/ξ
∂
−
1
q
∂
x
−
1
(0
,ξ
)
+
(
−
1)
j
+
1
ξ
j
μ
j
+
1
t
e
−
μ
t
x
(
v
)
/ξ
×
j
m
=
1
μ
t
x
(
v
)
ξ
m
j
=
0
b
m
,
x
(
v
)
x
(
v
)
0
e
μ
t
y
/ξ
∂
j
+
1
∂
y
j
+
1
q
(
y
,ξ
)
dy
+
j
m
=
1
j
−
m
α
=
1
j
−
m
β
=
0
m
−
1
s
=
0
m
−
1
−
s
γ
=
1
j
−
m
δ
=
0
a
j
,
m
,α,β,γ ,δ
×
μ
t
x
(
v
)
ξ
α
+
γ
x
(
v
)
β
+
δ
∂
s
∂
v
s
∂
j
+
1
∂
x
j
+
1
q
(
x
(
v
)
,ξ
)
dx
(
v
)
d
v
.
Since, by assumption, the function
q
is smooth, with, say,
|
∂
∂
x
q
(
x
,ξ
)
|
<
C
for each integer
and all relevant values of
x
and
ξ
,
we have
x
0
e
μ
t
y
/ξ
∂
∂
y
q
(
y
,ξ
)
dy
C
ξ
μ
t
(
e
μ
t
x
/ξ
−
1)
,
which, upon substitution in the previous expression with
=
j
+
1, presents
∂
j
∂
v
j
L
[
q
](
x
(
v
)
,ξ
) as a sum of terms containing
nonnegative powers of the bounded quantities
ξ
,
x
(
v
), as well as derivatives of
x
(
v
) and derivatives of
q
, all of which are also
uniformly bounded.
To complete the inductive proof we assume that, for all integers
j
,the
j
derivative
∂
j
v
u
m
of the
m
th Neumann series term
u
m
is bounded in the boundary-layer regions, and we show that the same is true for the (
m
+
1)th term
u
m
+
1
.Todothis,wefirst
differentiate Eq. (
A11
) to obtain
∂
j
∂
v
j
u
m
+
1
(
x
(
v
)
,ξ
)
=
∂
j
∂
v
j
(
R
0
(
ξ
)
u
m
(0
,
−
ξ
)
e
−
μ
t
x
(
v
)
/ξ
)
+
∂
j
∂
v
j
(
μ
s
L
[
S
[
u
m
]](
x
(
v
)
,ξ
))
.
(35)
The derivatives in the first term on the right-hand side were already shown to be bounded as part of the
m
=
0 proof [cf. Eq. (
34
)
and associated text]. For the second right-hand term, in turn, using Eqs. (
23
) and (
A21
) with
f
=
m
we obtain
∂
j
∂
v
j
L
[
S
[
u
m
]](
x
(
v
)
,ξ
)
=
j
k
=
1
j
−
k
w
,α
=
1
k
−
1
s
=
0
k
−
s
−
1
β
=
0
s
=
0
s
−
δ,γ
=
1
d
j
,
k
,
w
,α,β,δ,γ ,
s
,
X
δ
+
w
+
1
μ
t
(
x
(
v
)
−
1)
k
−
s
−
β
x
(
v
)
α
+
β
+
γ
∂
∂
v
m
(
x
(
v
)
,ξ
)
+
e
−
μ
t
x
/ξ
ξ
x
(
v
)
0
e
μ
t
y
/ξ
m
(
y
,ξ
)
dy
j
k
=
1
X
k
j
=
0
b
k
,
.
In view of the induction hypothesis, Eq. (
23
) and Remark
2
we
see that the
th derivative term on the right-hand side of this
equation is bounded. Since
m
(
y
,ξ
) is also bounded [
j
=
0
in Eq. (
24
)], however, we see that
x
(
v
)
0
e
μ
t
y
/ξ
m
(
y
,ξ
)
dy
E
m
0
x
(
v
)
e
μ
t
x
(
v
)
/ξ
,
for some constant
E
m
0
which, upon substitution on the right-
hand produces a term bounded by a constant time
|
X
|
. Since
|
x
(
v
)
|
1, it follows that the right-hand side in this equation
is bounded by a linear combination of powers of
|
X
|
. Since
additionally
|
X
|
<
1 in the boundary-layer region, it follows
that the left-hand side in Eq. (
35
) is bounded in the boundary-
layer region as well, and the proof is complete.
Figure
6
illustrates this result by displaying the (clearly
bounded)
v
derivatives for a numerical solution of the full
transport problem in the complete spatioangular domain. A
generalization of this result to the case of spatially varying
025306-8
BOUNDARY-LAYER STRUCTURES ARISING IN LINEAR ...
PHYSICAL REVIEW E
110
, 025306 (2024)
parameters
μ
s
(
x
),
μ
a
(
x
), and
μ
t
(
x
) proceeds similarly, pro-
vided the derivatives of these functions are adequately
bounded.
IV. NUMERICAL METHODS
This section introduces numerical methods, based on the
theoretical results presented in Secs.
III C
and
III D
, for the nu-
merical solution of the time-independent and time-dependent
transport problems (
1
) and (
3
). Throughout this section, the
solution of the transport equation under the changes of vari-
ables (
16
) and (
18
) will be denoted either by
U
(
v
,
r
)
=
u
(
x
(
v
)
,
(
r
))
(36)
or by
U
(
v
,
r
,
t
)
=
u
(
x
(
v
)
,
(
r
)
,
t
)
,
(37)
depending on whether the time-independent or time-
dependent problem is considered.
A. Transport problem in the (
v,
r
)and(
v,
r
,
t
) variables
Using Eq. (
36
), upon application of the changes of
variables (
16
) and (
18
), the time independent transport prob-
lem (
1
) becomes
(2
+
2 cosh(
v
))
(
r
)
∂
∂
v
U
(
v
,
r
)
+
μ
t
(
x
(
v
))
U
(
v
,
r
)
=
μ
s
(
x
(
v
))
1
−
1
p
(
(
r
)
,
(
r
))
U
(
v
,
r
)
d
dr
(
r
)
dr
+
q
(
x
(
v
)
,
(
r
))
,
U
(
−∞
,
r
)
=
R
0
(
(
r
))
U
(
−∞
,
−
1
(
R
(
r
)))
+
0
(
(
r
))
,
0
r
1
,
U
(
∞
,
r
)
=
R
1
(
(
r
))
U
(
∞
,
−
1
(
R
(
r
)))
+
1
(
(
r
))
,
−
1
r
<
0
.
(38)
A similar expression results under such changes of variables for the time dependent transport problem (
3
):
1
c
∂
∂
t
+
(
r
)(2
+
2 cosh(
v
))
∂
∂
v
+
μ
t
(
x
(
v
))
U
(
v
,
r
,
t
)
=
μ
s
(
x
(
v
))
1
−
1
p
(
(
r
)
,
(
r
))
U
(
v
,
r
,
t
)
d
dr
(
r
)
dr
+
q
(
x
(
v
)
,
(
r
)
,
t
)
,
U
(
v
,
r
,
t
=
0)
=
0
,
U
(
−∞
,
r
,
t
)
=
R
0
(
ξ
)
U
(
−∞
,
−
1
(
R
(
r
))
,
t
)
+
0
(
(
r
)
,
t
)
,
0
r
1
,
U
(
∞
,
r
,
t
)
=
R
1
(
ξ
)
U
(
∞
,
−
1
(
R
(
r
))
,
t
)
+
1
(
(
r
)
,
t
)
,
−
1
r
<
0
.
(39)
Naturally, Eq. (
39
) is the relevant equation for time-
dependent problems, which, containing transient data in-
formation, provides in many cases the most useful model
for the solution of the inverse transport problem. Addition-
ally, this equation may be useful even for solution of the
time-independent problem—via time relaxation. In detail,
considering, e.g., boundary condition functions of the form
0
(
(
r
)
,
t
)
=
T
(
t
)
0
(
(
r
)) and
1
(
(
r
)
,
t
)
=
T
(
t
)
1
(
(
r
))
[where
T
(
t
) denotes a suitable time profile that smoothly
transitions from
T
(0)
=
0to
T
(
t
1
)
=
1forsome
t
1
>
0] and
evolving the system up to sufficiently large times
t
>
t
1
over
several hundred (respectively, several thousand) time steps for
low (respectively, large) values of the scattering coefficient
μ
s
,
results in convergence to the desired stationary solution. Alter-
natively, the time-independent problem can be solved directly
on the basis of direct discretization of the time-independent
equation (
38
). The time-relaxation approach, which does not
require inversion of the full spatioangular matrix system, does
depend on implicit solution of the time-dependent problem
for a sufficiently long time period, as described above. While
the time-independent-equation approach does not require time
evolution, in turn, the computational cost for the inversion of
the system matrix grows quickly as the discretization is re-
fined. In practice we have found excellent agreement between
the results provided by these approaches for the solution of
time-independent problems. Roughly speaking, further, we
have found that for high-accuracy and
/
or low-to-moderate
values of
μ
s
the approach based on the time-dependent
equation is preferable, while use of the time-independent
equation is advantageous for large values of the scattering
coefficient—which, in the time dependent-equation approach
demands long relaxation times.
The proposed discretized version of the time-independent
problem (
38
) can be expressed in the form
[
D
+
μ
t
I
−
μ
s
S
]
u
=
q
,
(40)
where
D
,
, and
S
denote a discrete differential operator in
the variable
v
; a matrix corresponding to the coefficient (2
+
2 cosh(
v
))
(
r
) multiplying the
v
derivative in the equation;
and the discretized scattering operator introduced in Eq. (
42
)
below, respectively. The discrete version
D
of
∂
∂
v
operator is
obtained by direct differentiation of Fourier series obtained
by means of the Fourier-continuation method (FC) [
44
–
47
],
which enables representation of general smooth nonperi-
odic functions by Fourier-series with high accuracy and
negligible numerical dispersion. The quantities
u
=
(
u
i
,
m
)
≈
(
u
(
x
(
v
i
)
,
(
r
m
))) and
q
=
(
q
i
,
m
)
=
(
q
(
x
(
v
i
)
,
(
r
m
))), in turn,
denote the numerical approximation of the solution and the
source function at the points (
x
(
v
i
)
,
(
r
m
)), where
v
i
and
r
m
denote the discretization points in the
v
and
r
variables,
respectively:
v
i
(respectively,
r
m
) provides a uniform dis-
cretization of the domain
−
v
min
v
i
v
max
(respectively,
multi-interval Gauss-Legendre discretizations, as detailed in
Sec.
IV B
).
The discrete scattering operator
S
is used to incorporate
in the discrete setting the scattering integral
displayed
in Eq. (
2
), which occurs on the right-hand side of Eq. (
1
).
In detail, upon application of the change of variables (
18
)
025306-9