{psinputs}

Transport models for wave propagation in scattering media with nonlinear absorption

Joseph Kraisler Department of Applied Physics and Applied Mathematics, Columbia University, New York, NY, 10027; jek2199@columbia.edu Wei Li Department of Mathematical Sciences, Depaul University, Chicago, IL 60604; wei.li@depaul.edu Kui Ren Department of Applied Physics and Applied Mathematics, Columbia University, New York, NY 10027;kr2002@columbia.edu John C. Schotland Department of Mathematics and Department of Physics, Yale University, New Haven, CT 06511; john.schotland@yale.edu Yimin Zhong Department of Mathematics and Statistics, Auburn University, Auburn, AL 36830; yimin.zhong@auburn.edu

Abstract

This work considers the propagation of high-frequency waves in highly-scattering media where physical absorption of a nonlinear nature occurs. Using the classical tools of the Wigner transform and multiscale analysis, we derive semilinear radiative transport models for the phase-space intensity and the diffusive limits of such transport models. As an application, we consider an inverse problem for the semilinear transport equation, where we reconstruct the absorption coefficients of the equation from a functional of its solution. We obtain a uniqueness result on the inverse problem.

Key words. wave propagation, nonlinear media, nonlinear absorption, semilinear radiative transport equation, semilinear diffusion equation, inverse problem

1 Introduction

The derivation of kinetic models for wave propagation in highly-scattering media is a classical subject [32, 15] that has received significant attention in the past two decades due to its importance in many emerging applications [11, 13, 25, 39, 41, 35]. A significant amount of progress has been made on both the mathematical justification of the derivation (such as those based on multiscale analysis of the Wigner transform) [6, 21, 24, 30, 37, 49] and computational validation of the derived kinetic models [7, 9, 33, 47, 53]; see [1, 2, 4, 5, 12, 23, 25, 28] and references therein for additional investigations in this field. The obtained models for imaging in complex media have also been utilized in many different settings [8, 9, 13, 16].

In this work, we are interested in developing kinetic models when nonlinear absorption occurs during wave propagation [10, 38, 48, 57]. We are mainly motivated by applications where reconstructing the absorption of the underlying medium from internal or boundary observations is of practical interest. Such applications include, for instance, the case of reconstructing the two-photon absorption coefficient of biological tissues with optical or photoacoustic measurements [20, 36, 42, 44, 45, 51, 54, 55, 56].

Our derivation will be carried out in the frequency domain, where wave propagation is described by the standard Helmholtz equation. The nonlinear absorption mechanism we are interested we study is modeled as a zeroth-order perturbation to the second-order Helmholtz differential operator. This is the essential factor that makes it possible for us to perform the standard calculations in this field for our nonlinear problem. Even though this is a mathematically less attractive nonlinearity to study, the derivation does provide us a formal justification of the semilinear radiative transport model, see, for instance (20), used in many applications. We refer interested readers to [22] for the derivation of transport models for a different type of nonlinearity that makes use a mean-field approximation.

The rest of this paper is organized as follows. In Section 2, we derive the radiative transport model for media with quadratic and higher-order absorption. We then discuss the diffusive limit of the derived transport models in Section 3. As an application, we study in Section 4 an inverse medium problem for our semilinear radiative transport equation. Concluding remarks are offered in Section 5.

2 Derivation of the transport equation

For simplicity, we first consider the case of quadratic nonlinear absorption. This could serve as a model of light propagation in media with two-photon absorption [51, 54]. We will then generalize the result to the case of a general polynomial nonlinearity. Let the wave field $p(z,\mathbf{x})$ be the solution to the scalar wave equation in the time-harmonic form, that is,

\Delta_{\mathbf{x}}p+\frac{\partial^{2}p}{\partial z^{2}}+k^{2}n^{2}p=0,

(1)

where $\Delta_{\mathbf{x}}$ is the transverse Laplacian in $\mathbf{x}\in\mathbb{R}^{d}$ ( $d\geq 1$ ), $k$ is the wave number, and $n=n(z,\mathbf{x},p)$ is the refractive index. We assume that the refractive index takes the form

n^{2}=1-2\sigma V(\frac{z}{\ell_{z}},\frac{\mathbf{x}}{\ell_{\mathbf{x}}})+ik^{-1}\mu\widetilde{K}(z,\mathbf{x})|p|^{2}\,,

(2)

where $V$ is a real bounded stationary random field with zero mean, with $\ell_{\mathbf{x}}$ and $\ell_{z}$ are the transverse and longitudinal correlation lengths of the random field, respectively. The deterministic function $\widetilde{K}$ is non-negative and measures the strength of the second-order absorption. The parameters $\sigma$ and $\mu$ are the scaling factors quantifying the amplitudes of the fluctuation and the second-order absorption, respectively. Assuming that the field $p$ possesses a beam-like structure propagating in $z$ direction and neglecting the back scattering, we may write $p(z,\mathbf{x})=e^{ikz}\psi(z,\mathbf{x})$ with complex amplitude $\psi(z,\mathbf{x})$ satisfing the following equation

\frac{\partial^{2}\psi}{\partial z^{2}}+2ik\frac{\partial\psi}{\partial z}+\Delta_{\mathbf{x}}\psi+k^{2}(n^{2}-1)\psi=0\,.

(3)

Let $L_{\mathbf{x}}$ and $L_{z}$ be the characteristic lengths of propagation in $\mathbf{x}$ and $z$ directions, respectively. We rescale the variables $\mathbf{x}\mapsto L_{\mathbf{x}}\mathbf{x}$ and $z\mapsto L_{z}z$ and $\mathbf{x}$ , $z$ are now $\mathcal{O}(1)$ . Then with the newly-defined variables, we may write $(n^{2}-1)$ as

n^{2}-1=-2\sigma V(\frac{L_{z}z}{\ell_{z}},\frac{L_{\mathbf{x}}\mathbf{x}}{\ell_{\mathbf{x}}})+ik^{-1}\mu\widetilde{K}(L_{z}\mathbf{z},L_{\mathbf{x}}\mathbf{x})|\psi|^{2}\,.

(4)

The equation (3) now becomes

\frac{1}{L_{z}^{2}}\frac{\partial^{2}\psi}{\partial z^{2}}+2\frac{ik}{L_{z}}\frac{\partial\psi}{\partial z}+\frac{1}{L_{\mathbf{x}}^{2}}\Delta_{\mathbf{x}}\psi-2k^{2}\sigma V(\frac{L_{z}z}{\ell_{z}},\frac{L_{\mathbf{x}}\mathbf{x}}{\ell_{\mathbf{x}}})\psi+ik\mu\widetilde{K}(L_{z}z,L_{\mathbf{x}}\mathbf{x})|\psi|^{2}\psi=0\,.

(5)

With the small aperture assumption that $L_{\mathbf{x}}\ll L_{z}$ , we can formally approximate the above equation by the paraxial wave equation

i\frac{\partial\psi}{\partial z}+\frac{L_{z}}{2kL_{\mathbf{x}}^{2}}\Delta_{\mathbf{x}}\psi-kL_{z}\sigma V(\frac{L_{z}z}{\ell_{z}},\frac{L_{\mathbf{x}}\mathbf{x}}{\ell_{\mathbf{x}}})\psi+i\frac{L_{z}}{2}\mu\widetilde{K}(L_{z}z,L_{\mathbf{x}}\mathbf{x})|\psi|^{2}\psi=0\,.

(6)

Our derivation works in the regime where the longitudinal propagation distance $L_{z}$ is much larger than the correlation length $\ell_{z}$ and the correlation length is much larger than the wavelength, that is $L_{z}\gg\ell_{z}$ and $\ell_{z}\gg\lambda:=\frac{2\pi}{k}$ . We, therefore, introduce the small parameter $\varepsilon$ , and assume the scaling relations in the weak-coupling regime:

\frac{\ell_{\mathbf{x}}}{L_{\mathbf{x}}}=\frac{\ell_{z}}{L_{z}}=\varepsilon\ll 1,\quad k\ell_{\mathbf{x}}^{2}=\ell_{z},\quad\sigma=\frac{1}{k\ell_{z}}\sqrt{\varepsilon},\quad\mu\propto\frac{1}{L_{z}}.

(7)

Let us denote the rescaled wave field by $\psi_{\varepsilon}(z,\mathbf{x})$ and take $K(z,\mathbf{x})=L_{z}\mu\widetilde{K}(L_{z}z,L_{\mathbf{x}}\mathbf{x})$ . Then the paraxial wave equation turns into

i\frac{\partial\psi_{\varepsilon}}{\partial z}+\frac{\varepsilon}{2}\Delta_{\mathbf{x}}\psi_{\varepsilon}-\frac{1}{\sqrt{\varepsilon}}V(\frac{z}{\varepsilon},\frac{\mathbf{x}}{\varepsilon})\psi_{\varepsilon}+\frac{i}{2}K(z,\mathbf{x})|\psi_{\varepsilon}|^{2}\psi_{\varepsilon}=0.

(8)

We then take the Wigner transform of $\psi_{\varepsilon}$ :

W_{\varepsilon}(z,\mathbf{x},\mathbf{k})=\int_{\mathbb{R}^{d}}e^{i\mathbf{k}\cdot\mathbf{y}}\psi_{\varepsilon}(\mathbf{x}-\frac{\varepsilon\mathbf{y}}{2},z)\overline{\psi_{\varepsilon}(\mathbf{x}+\frac{\varepsilon\mathbf{y}}{2},z)}\frac{d\mathbf{y}}{(2\pi)^{d}},

(9)

It is then standard to check $W_{\varepsilon}(z,\mathbf{x},\mathbf{k})$ satisfies the Liouville equation:

\frac{\partial W_{\varepsilon}}{\partial z}+\mathbf{k}\cdot\nabla_{\mathbf{x}}W_{\varepsilon}+\mathcal{L}_{V,\varepsilon}W_{\varepsilon}+\mathcal{L}_{K,\varepsilon}W_{\varepsilon}=0,

(10)

where

	$\displaystyle\mathcal{L}_{V,\varepsilon}W_{\varepsilon}$	$\displaystyle=\frac{i}{\sqrt{\varepsilon}}\int_{\mathbb{R}^{d}}e^{-i\mathbf{p}\cdot\mathbf{x}/\varepsilon}\left(W_{\varepsilon}(z,\mathbf{x},\mathbf{k}+\frac{\mathbf{p}}{2})-W_{\varepsilon}(z,\mathbf{x},\mathbf{k}-\frac{\mathbf{p}}{2})\right)\frac{\widehat{V}(\frac{z}{\varepsilon},\mathbf{p})d\mathbf{p}}{(2\pi)^{d}},$		(11)
	$\displaystyle\mathcal{L}_{K,\varepsilon}W_{\varepsilon}$	$\displaystyle=\frac{1}{2}\int_{\mathbb{R}^{d}}e^{-i\mathbf{p}\cdot x}\left(W_{\varepsilon}(z,\mathbf{x},\mathbf{k}+\frac{\varepsilon\mathbf{p}}{2})+W_{\varepsilon}(z,\mathbf{x},\mathbf{k}-\frac{\varepsilon\mathbf{p}}{2})\right)\frac{\widehat{S}_{\varepsilon}(z,\mathbf{p})d\mathbf{p}}{(2\pi)^{d}}\,.$		(11)

Here $S_{\varepsilon}(z,\mathbf{x})$ is defined as

\displaystyle S_{\varepsilon}(z,\mathbf{x}):=K(z,\mathbf{x})|\psi_{\varepsilon}(z,\mathbf{x})|^{2}=K(z,\mathbf{x})\int_{\mathbb{R}^{d}}W_{\varepsilon}(z,\mathbf{x},\mathbf{k})d\mathbf{k}\,,

(12)

while $\widehat{V}$ and $\widehat{S}_{\varepsilon}$ denote the Fourier transform ( $\mathbf{x}\to\mathbf{p}$ ) of $V$ and $S_{\varepsilon}$ respectively. We use the standard Fourier transform definition

\widehat{f}(\mathbf{p})=\int_{\mathbb{R}^{d}}e^{i\mathbf{p}\cdot\mathbf{x}}f(\mathbf{x})\frac{d\mathbf{x}}{(2\pi)^{d}}\,.

2.1 Multiscale expansion

In order to find the asymptotic limit as $\varepsilon\to 0$ , we introduce $\mathbf{y}=\mathbf{x}/\varepsilon$ as the fast variable and denote $W_{\varepsilon}(z,\mathbf{x},\mathbf{k})=W_{\varepsilon}(z,\mathbf{x},\mathbf{y},\mathbf{k})$ . Formally we write $W_{\varepsilon}$ in its asymptotic expansion in $\varepsilon$ :

W_{\varepsilon}(z,\mathbf{x},\mathbf{y},\mathbf{k})=W_{0}(z,\mathbf{x},\mathbf{y},\mathbf{k})+\sqrt{\varepsilon}W_{1}(z,\mathbf{x},\mathbf{y},\mathbf{k})+\varepsilon W_{2}(z,\mathbf{x},\mathbf{y},\mathbf{k})+\cdots\,.

Using (12), we may also expand $S_{\varepsilon}(z,\mathbf{x})=S_{\varepsilon}(z,\mathbf{x},\mathbf{y})$ accordingly as

S_{\varepsilon}(z,\mathbf{x},\mathbf{y})=S_{0}(z,\mathbf{x},\mathbf{y})+\sqrt{\varepsilon}S_{1}(z,\mathbf{x},\mathbf{y})+\varepsilon S_{2}(z,\mathbf{x},\mathbf{y})+\cdots\,.

where

S_{0}(z,\mathbf{x},\mathbf{y})=K(z,\mathbf{x})\int_{\mathbb{R}^{d}}W_{0}(z,\mathbf{x},\mathbf{y},\mathbf{k})d\mathbf{k}\,.

We can now plug in the transform $\nabla_{\mathbf{x}}\to\nabla_{\mathbf{x}}+\frac{1}{\varepsilon}\nabla_{\mathbf{y}}$ into (10) to conclude that the leading order equation at $\mathcal{O}(\varepsilon^{-1})$ implies $\mathbf{k}\cdot\nabla_{\mathbf{y}}W_{0}=0$ . This is equivalent to

W_{0}(z,\mathbf{x},\mathbf{y},\mathbf{k})=W_{0}(z,\mathbf{x},\mathbf{k}),\ \ \ \mbox{and}\ \ \ S_{0}(z,\mathbf{x},\mathbf{y})=S_{0}(z,\mathbf{x})\,.

For the order of $\mathcal{O}(\varepsilon^{-1/2})$ , we have

\mathbf{k}\cdot\nabla_{\mathbf{y}}W_{1}+\alpha W_{1}+i\int_{\mathbb{R}^{d}}e^{-i\mathbf{p}\cdot\mathbf{y}}(W_{0}(z,\mathbf{x},\mathbf{k}+\frac{\mathbf{p}}{2})-W_{0}(z,\mathbf{x},\mathbf{k}-\frac{\mathbf{p}}{2}))\widehat{V}(\frac{z}{\varepsilon},\mathbf{p})\frac{d\mathbf{p}}{(2\pi)^{d}}=0\,,

(13)

where $\alpha\to 0^{+}$ . This gives that the Fourier transform of $W_{1}$ , $\widehat{W}_{1}$ , is

\widehat{W}_{1}(z,\mathbf{x},\mathbf{p},\mathbf{k})=\frac{(W_{0}(z,\mathbf{x},\mathbf{k}+\frac{\mathbf{p}}{2})-W_{0}(z,\mathbf{x},\mathbf{k}-\frac{\mathbf{p}}{2}))\widehat{V}(\frac{z}{\varepsilon},\mathbf{p})}{\mathbf{p}\cdot\mathbf{k}+i\alpha}\,.

(14)

Finally, we derive the equation for $\mathcal{O}(1)$ terms. To handle the nonlinearity, we let $s>d+2$ and assume exist positive constants $C_{1}$ , $C_{2}$ and $C_{3}$ such that

	$\displaystyle\\|W_{0}(z,\mathbf{x},\cdot)\\|_{C^{1}(\mathbb{R}^{d})}+\\|W_{0}(z,\mathbf{x},\cdot)\\|_{L^{1}(\mathbb{R}^{d})}<C_{1},$	$\displaystyle\forall(z,\mathbf{x})\in\mathbb{R}^{d+1},$
	$\displaystyle\left\\|\int_{\mathbb{R}^{d}}W_{0}(z,\cdot,\mathbf{k})d\mathbf{k}\right\\|_{H^{s}(\mathbb{R}^{d})}<C_{2},$	$\displaystyle\forall z\in\mathbb{R},$
	$\displaystyle\\|K(z,\cdot)\\|_{H^{s}(\mathbb{R}^{d})}<C_{3},$	$\displaystyle\forall z\in\mathbb{R}\,.$

Then we have that

W_{0}(z,\mathbf{x},\mathbf{k}-\frac{\varepsilon}{2}\mathbf{p})+W_{0}(z,\mathbf{x},\mathbf{k}+\frac{\varepsilon}{2}\mathbf{p})=2W_{0}(z,\mathbf{x},\mathbf{k})+\mathcal{O}(\varepsilon|\mathbf{p}|)\,,

(15)

where $\mathcal{O}(\varepsilon|\mathbf{p}|)$ has a uniform constant. Moreover, we have that

$\displaystyle\left\|\int_{\mathbb{R}^{d}}\frac{\|\mathbf{p}\|\|\widehat{S}_{0}(z,\mathbf{p})\|d\mathbf{p}}{(2\pi)^{d}}\right\|^{2}$	$\displaystyle\leq\frac{1}{(2\pi)^{2d}}\int_{\mathbb{R}^{d}}\frac{\|\mathbf{p}\|^{2}}{(1+\|\mathbf{p}\|^{2})^{s}}d\mathbf{p}\int_{\mathbb{R}^{d}}(1+\|\mathbf{p}\|^{2})^{\frac{s}{2}}\|\widehat{S}_{0}(z,\mathbf{p})\|^{2}d\mathbf{p}$	(16)
	$\displaystyle=\frac{1}{(2\pi)^{2d}}\left(\int_{\mathbb{R}^{d}}\frac{\|\mathbf{p}\|^{2}}{(1+\|\mathbf{p}\|^{2})^{s}}d\mathbf{p}\right)\\|S_{0}(z,\cdot)\\|_{H^{s}(\mathbb{R}^{d})}^{2}$
	$\displaystyle\leq C_{s}\\|K(z,\cdot)\\|_{H^{s}(\mathbb{R}^{d})}^{2}\left\\|\int_{\mathbb{R}^{d}}W_{0}(z,\cdot,\mathbf{k})d\mathbf{k}\right\\|_{H^{s}(\mathbb{R}^{d})}^{2}$

is also uniformly bounded, where $C_{s}>0$ is a constant depending only on $s$ . The last inequality holds since $s>d+2$ implies $s>d/2$ . Therefore the $\mathcal{O}(1)$ term is

		$\displaystyle\frac{\partial W_{0}}{\partial z}+\mathbf{k}\cdot\nabla_{\mathbf{x}}W_{0}+\mathbf{k}\cdot\nabla_{\mathbf{y}}W_{2}$		(17)
		$\displaystyle+i\int_{\mathbb{R}^{d}}e^{-i\mathbf{p}\cdot\mathbf{y}}(W_{1}(z,\mathbf{x},\mathbf{y},\mathbf{k}+\frac{\mathbf{p}}{2})-W_{1}(z,\mathbf{x},\mathbf{y},\mathbf{k}-\frac{\mathbf{p}}{2}))\widehat{V}(\frac{z}{\varepsilon},\mathbf{p})\frac{d\mathbf{p}}{(2\pi)^{d}}$
		$\displaystyle+W_{0}(z,\mathbf{x},\mathbf{k})S_{0}(z,\mathbf{x})=0\,.$

In order to close the equation, we still need to add the orthogonal relation between $W_{0}$ and $W_{2}$ , that is, $\mathbb{E}[\mathbf{k}\cdot\nabla_{\mathbf{y}}W_{2}]=0$ . Hence, we have

		$\displaystyle\frac{\partial W_{0}}{\partial z}+\mathbf{k}\cdot\nabla_{\mathbf{x}}W_{0}$		(18)
		$\displaystyle+\mathbb{E}\left[i\int_{\mathbb{R}^{d}}e^{-i\mathbf{p}\cdot\mathbf{y}}(W_{1}(z,\mathbf{x},\mathbf{y},\mathbf{k}+\frac{\mathbf{p}}{2})-W_{1}(z,\mathbf{x},\mathbf{y},\mathbf{k}-\frac{\mathbf{p}}{2}))\widehat{V}(\frac{z}{\varepsilon},\mathbf{p})\frac{d\mathbf{p}}{(2\pi)^{d}}\right]$
		$\displaystyle+W_{0}(z,\mathbf{x},\mathbf{k})S_{0}(z,\mathbf{x})=0\,.$

Let $R$ be the correlation function of $V$ , and assume that the power spectrum satisfies

\mathbb{E}[\widehat{V}(z,\mathbf{p})\widehat{V}(z,\mathbf{q})]=(2\pi)^{d}\widehat{R}(\mathbf{p})\delta(\mathbf{p}+\mathbf{q})\,.

(19)

Then the expectation term in (18) converges weakly to

		$\displaystyle\mathbb{E}\left[i\int_{\mathbb{R}^{d}}e^{-i\mathbf{p}\cdot\mathbf{x}/\varepsilon}(W_{1}(z,\mathbf{x},\frac{\mathbf{x}}{\varepsilon},\mathbf{k}+\frac{\mathbf{p}}{2})-W_{1}(z,\mathbf{x},\frac{\mathbf{x}}{\varepsilon},\mathbf{k}-\frac{\mathbf{p}}{2}))\widehat{V}(\frac{z}{\varepsilon},\mathbf{p})\frac{d\mathbf{p}}{(2\pi)^{d}}\right]$
		$\displaystyle\to 4\pi\int_{\mathbb{R}^{d}}\widehat{R}(\mathbf{p}-\mathbf{k})[W_{0}(z,\mathbf{x},\mathbf{k})-W_{0}(z,\mathbf{x},\mathbf{p})]\delta(\|\mathbf{k}\|^{2}-\|\mathbf{p}\|^{2})d\mathbf{p}\,.$

Therefore the final radiative transport equation of $W_{0}$ is

		$\displaystyle\frac{\partial W_{0}}{\partial z}+\mathbf{k}\cdot\nabla_{\mathbf{x}}W_{0}+4\pi\int_{\mathbb{R}^{d}}\widehat{R}(\mathbf{p}-\mathbf{k})[W_{0}(z,\mathbf{x},\mathbf{k})-W_{0}(z,\mathbf{x},\mathbf{p})]\delta(\|\mathbf{k}\|^{2}-\|\mathbf{p}\|^{2})d\mathbf{p}$		(20)
		$\displaystyle+\left(K(z,\mathbf{x})\int_{\mathbb{R}^{d}}W_{0}(z,\mathbf{x},\mathbf{k}^{\prime})d\mathbf{k}^{\prime}\right)W_{0}(z,\mathbf{x},\mathbf{k})=0\,.$		(20)

The last term on the left side of (20) stands for the quadratic absorption, which shows that the absorption coefficient linearly depends on the angular average of $W_{0}(\mathbf{x},z,\mathbf{k})$ .

2.2 Extension to higher order absorption

We now extend the previous result to the case of general polynomial absorption. This could be a physical model for the propagation of light in media with multi-photon absorption [14]. The absorption term in the refractive index is modeled by

n^{2}=1-2\sigma V(\frac{z}{\ell_{z}},\frac{\mathbf{x}}{\ell_{\mathbf{x}}})+ik^{-1}\mu\sum_{l=0}^{m}\widetilde{K}_{l}(z,\mathbf{x})|\psi|^{2l},

(21)

where $\widetilde{K}_{l}$ stands for the absorption strength of $(l+1)$ -th order. Following the same derivation of the paraxial wave equation, we have a new Liouville equation in the form of (10) where the only modification is the term

S_{\varepsilon}(\mathbf{x},z)=\sum_{l=0}^{L}K_{l}(z,\mathbf{x})|\psi_{\varepsilon}(z,\mathbf{x})|^{2l},\quad K_{l}(z,\mathbf{x}):=\mu L_{z}\widetilde{K}_{l}(L_{z}z,L_{\mathbf{x}}\mathbf{x})\,.

(22)

Assuming that $K_{l}$ and $W_{0}$ are sufficiently regular, we can follow the same procedure and obtain the radiative transport equation for $W_{0}$ in the setting of the polynomial absorption:

		$\displaystyle\frac{\partial W_{0}}{\partial z}+\mathbf{k}\cdot\nabla_{\mathbf{x}}W_{0}+4\pi\int_{\mathbb{R}^{d}}\widehat{R}(\mathbf{p}-\mathbf{k})[W_{0}(z,\mathbf{x},\mathbf{k})-W_{0}(z,\mathbf{x},\mathbf{p})]\delta(\|\mathbf{k}\|^{2}-\|\mathbf{p}\|^{2})d\mathbf{p}$		(23)
		$\displaystyle+\left(\sum_{l=0}^{L}K_{l}(z,\mathbf{x})\left[\int_{\mathbb{R}^{d}}W_{0}(z,\mathbf{x},\mathbf{k}^{\prime})d\mathbf{k}^{\prime}\right]^{l}\right)W_{0}(z,\mathbf{x},\mathbf{k})=0\,.$		(23)

2.3 The nonlinear RTE in the standard form

For a monochromatic solution $W_{0}(z,\mathbf{x},\mathbf{k})$ which is supported on $\|\mathbf{k}\|=1$ , the above radiative transport equation becomes

		$\displaystyle\frac{\partial W_{0}}{\partial z}+\mathbf{k}\cdot\nabla_{\mathbf{x}}W_{0}+4\pi\int_{\mathbb{S}^{d-1}}\widehat{R}(\mathbf{p}-\mathbf{k})[W_{0}(z,\mathbf{x},\mathbf{k})-W_{0}(z,\mathbf{x},\mathbf{p})]d\mathbf{p}$		(24)
		$\displaystyle+\left(\sum_{l=0}^{L}K_{l}(z,\mathbf{x})\left[\int_{\mathbb{S}^{d-1}}W_{0}(z,\mathbf{x},\mathbf{k}^{\prime})d\mathbf{k}^{\prime}\right]^{l}\right)W_{0}(z,\mathbf{x},\mathbf{k})=0\,,$		(24)

where $\mathbb{S}^{d-1}$ is the unit sphere in $\mathbb{R}^{d}$ . It can further be put in the standard form

\displaystyle\frac{\partial W_{0}}{\partial z}+\mathbf{k}\cdot\nabla_{\mathbf{x}}W_{0}+\Sigma_{a}(\langle{W_{0}}\rangle)W_{0}+\Sigma_{s}(W_{0}-\mathcal{K}W_{0})

\displaystyle=0.

(25)

Here the total energy of the field at $\mathbf{x}$ is denoted by

\langle{W_{0}}\rangle:=\int_{\mathbb{S}^{d-1}}W_{0}(z,\mathbf{x},\mathbf{k})d\mathbf{k}

and

\Sigma_{a}(\langle{W_{0}}\rangle)=\sum_{l=0}^{L}\Sigma_{a,l}\langle{W_{0}}\rangle^{l}

(26)

is the effective absorption coefficient with $\Sigma_{a,l}(z,\mathbf{x})=K_{l}(z,\mathbf{x})$ is the absorption coefficient of $(l+1)$ -th order. The scattering coefficient is $\Sigma_{s}(\mathbf{k}):=4\pi\displaystyle\int_{\mathbb{S}^{d-1}}\widehat{R}(\mathbf{p}-\mathbf{k})d\mathbf{p}$ and the scattering operator $\mathcal{K}$ is defined as

\mathcal{K}u(z,\mathbf{x},\mathbf{k}):=\int_{\mathbb{S}^{d-1}}p(\mathbf{k},\mathbf{k}^{\prime})u(z,\mathbf{x},\mathbf{k}^{\prime})d\mathbf{k}^{\prime}

(27)

where the scattering phase function

p(\mathbf{k},\mathbf{k}^{\prime})=\frac{4\pi\widehat{R}(\mathbf{k}-\mathbf{k}^{\prime})}{\Sigma_{s}(\mathbf{k})}.

For the problem to be physical, the absorption strengths $K_{l}$ ( $0\leq l\leq L$ ) should be all non-negative.

3 The diffusion limit

We now study the diffusion limit of the transport equation for $W_{0}(z,\mathbf{x},\mathbf{k})$ with general polynomial nonlinear absorption. We assume that the physical domain $\Omega$ is convex with smooth boundary $\partial\Omega$ . We focus on the following nonlinear radiative transport equation:

$\displaystyle\frac{\partial W_{0}}{\partial z}+\mathbf{k}\cdot\nabla_{\mathbf{x}}W_{0}+\Sigma_{a}(\langle{W_{0}}\rangle)W_{0}+\Sigma_{s}(W_{0}-\mathcal{K}W_{0})$	$\displaystyle=0\quad$	in	$\displaystyle X,$	(28)
$\displaystyle W_{0}(z,\mathbf{x},\mathbf{k})$	$\displaystyle=0\quad$	on	$\displaystyle\Gamma_{-},$
$\displaystyle W_{0}(0,\mathbf{x},\mathbf{k})$	$\displaystyle=f(\mathbf{x})\quad$	on	$\displaystyle X_{0},$

where $X=(0,T)\times\Omega\times\mathbb{S}^{d-1}$ , $\Gamma_{-}:=\{(z,\mathbf{x},\mathbf{k})\in(0,T)\times\partial\Omega\times\mathbb{S}^{d-1}\mid\mathbf{n}_{\mathbf{x}}\cdot\mathbf{k}<0\}$ , and $\mathbf{n}_{\mathbf{x}}$ is the unit outward normal vector at $\mathbf{x}\in\partial\Omega$ , $X_{0}=\Omega\times\mathbb{S}^{d-1}$ .

We will focus on power spectra of the form $\widehat{R}(\mathbf{p}-\mathbf{k})=\widehat{R}(\mathbf{p}\cdot\mathbf{k})$ , in which case the scattering coefficient

\Sigma_{s}(\mathbf{k})\equiv\Sigma_{s}

is a constant, and the scattering phase function is $p(\mathbf{k},\mathbf{k}^{\prime})=4\pi\widehat{R}(\mathbf{k}\cdot\mathbf{k}^{\prime})/\Sigma_{s}$ . Assume the scattering phase function $p(\mathbf{k}\cdot\mathbf{k}^{\prime})$ is bounded below and above by positive constants $\underline{\theta},\overline{\theta}\geq 0$ such that

\underline{\theta}\leq p(\mathbf{k},\mathbf{k}^{\prime})\leq\overline{\theta}\,.

(29)

Note that $\underline{\theta}$ satisfies the condition $\nu_{d-1}\underline{\theta}\leq 1$ , where $\nu_{d-1}$ is the measure of the ( $d-1$ ) dimensional unit sphere. For simplicity, we require that the initial condition $f(\mathbf{x})\in L^{\infty}(\Omega\times\mathbb{S}^{d-1})$ is $\mathbf{x}$ -dependent and non-negative. The absorption coefficients obey the condition

0<\underline{\Sigma_{a,l}}\leq\Sigma_{a,l}(z,\mathbf{x})\leq\overline{\Sigma_{a,l}}<\infty,

(30)

for some constants $\underline{\Sigma_{a,l}}$ and $\overline{\Sigma_{a,l}}$ .

We need the following lemma to have a well-posed problem.

Lemma 3.1 (Kellogg [34]).

Let $\mathcal{M}$ be a bounded convex open subset of a real Banach space, and $F:\overline{\mathcal{M}}\to\overline{\mathcal{M}}$ is a compact continuous map which is continuously Fréchet differentiable on T. If (i) for each $m\in\mathcal{M}$ , $1$ is not an eigenvalue of $F^{\prime}(m)$ , and (ii) for each $m\in\partial\mathcal{M}$ , $m\neq F(m)$ , then $F$ has a unique fixed point in $\mathcal{M}$ .

We define the space, using $dS(\mathbf{x})$ to denote the surface measure of $\partial\Omega$ ,

	$\displaystyle\mathcal{W}_{p}=\{u\in L^{p}((0,T)\times\Omega\times\mathbb{S}^{d-1});\frac{\partial u}{\partial z}+\mathbf{k}\cdot\nabla u\in L^{p}((0,T)\times\Omega\times\mathbb{S}^{d-1});$
	$\displaystyle u(0,\cdot,\cdot)\in L^{p}(\Omega\times\mathbb{S}^{d-1});u\|_{\Gamma_{-}}\in L^{p}((0,T)\times\Gamma_{-},\|\mathbf{n}_{\mathbf{x}}\cdot\mathbf{k}\|dz\,d\mathbf{k}dS(\mathbf{x}))\}$

and the equipped norm

\|u\|_{\mathcal{W}_{p}}^{p}=\|u\|_{L^{p}((0,T)\times\Omega\times\mathbb{R}^{d-1})}^{p}+\|\partial_{z}u+\mathbf{k}\cdot\nabla u\|_{L^{p}((0,T)\times\Omega\times\mathbb{R}^{d-1})}^{p}.

Theorem 3.2.

Let the initial condition $f(\mathbf{x})$ satisfy:

(i)\qquad\|f\|_{\infty}\leq\frac{1}{\nu_{d-1}}\inf_{(z,\mathbf{x})\in(0,T)\times\Omega}\frac{\Sigma_{a,k-1}}{k\Sigma_{a,k}},\quad\forall 1\leq k\leq L

(ii)\qquad\nu_{d-1}\Sigma_{a}^{\prime}(\|f\|_{\infty})\|f\|_{\infty}\leq 2\Sigma_{s}\underline{\theta},\quad\forall(z,\mathbf{x})\in(0,T)\times\Omega,

where $\Sigma_{a}^{\prime}$ is the Fréchet derivative of $\Sigma_{a}$ , that is,

\Sigma_{a}^{\prime}(m)=\sum_{s=1}^{L}s\Sigma_{a,l}(z,\mathbf{x})m^{s-1}\,.

Then equation (28) admits a unique non-negative solution in $\mathcal{W}_{\infty}((0,T)\times\Omega\times\mathbb{S}^{d-1})$ .

Proof.

Let $F:m\mapsto\phi$ be the map defined through the relation $\langle{\phi}\rangle=Fm$ , where $\phi$ solves the transport equation

$\displaystyle\frac{\partial\phi}{\partial z}+\mathbf{k}\cdot\nabla\phi+\Sigma_{a}(m)\phi+\Sigma_{s}\left(\phi-\mathcal{K}\phi\right)$	$\displaystyle=0\quad$	in	$\displaystyle X,$	(31)
$\displaystyle\phi(z,\mathbf{x},\mathbf{k})$	$\displaystyle=0\quad$	on	$\displaystyle\Gamma_{-},$
$\displaystyle\phi(0,\mathbf{x},\mathbf{k})$	$\displaystyle=f(\mathbf{x})\quad$	on	$\displaystyle X_{0}\,.$

We denote by $\mathcal{M}$ the set

\mathcal{M}=\{m\in L^{\infty}((0,T)\times\Omega)\cap L^{2}((0,T)\times\Omega)\mid 0\leq m\leq\|f\|_{\infty}+\delta\}

with $\delta>0$ being arbitrary. It is straightforward to check that $\mathcal{M}$ is convex, bounded, and closed under the usual $L^{2}$ topology. For any $m\in\mathcal{M}$ , we have that $\Sigma_{a}(m)$ is non-negative. Therefore, by the maximum principle for the linear transport equation (31) (see for instance [19]), $\langle{\phi}\rangle\in\mathcal{M}$ . This shows that $F:\mathcal{M}\to\mathcal{M}$ . To show that $F$ is a continuous operator, we denote by $\phi_{1}$ and $\phi_{2}$ the solutions to (31) corresponding to $m_{1}\in\mathcal{M}$ and $m_{2}\in\mathcal{M}$ respectively. We then introduce $w=\phi_{1}-\phi_{2}$ . It is then clear that $w$ solves the linear transport equation:

$\displaystyle\frac{\partial w}{\partial z}+\mathbf{k}\cdot\nabla w+\Sigma_{a}(m_{1})w+\Sigma_{s}\left(w-\mathcal{K}w\right)$	$\displaystyle=(\Sigma_{a}(m_{2})-\Sigma_{a}(m_{1}))\phi_{2}\quad$	in	$\displaystyle X,$	(32)
$\displaystyle w(z,\mathbf{x},\mathbf{k})$	$\displaystyle=0\quad$	on	$\displaystyle\Gamma_{-},$
$\displaystyle w(0,\mathbf{x},\mathbf{k})$	$\displaystyle=0\quad$	on	$\displaystyle X_{0}\,.$

By standard linear theory [19], this equation admit a unique $w\in\mathcal{W}_{2}$ such that

	$\displaystyle\\|w\\|_{\mathcal{W}_{2}((0,T)\times\Omega\times\mathbb{S}^{d-1})}$	$\displaystyle\leq C\\|(\Sigma_{a}(m_{2})-\Sigma_{a}(m_{1}))\phi_{2}\\|_{L^{2}((0,T)\times\Omega\times\mathbb{S}^{d-1})}$		(33)
		$\displaystyle\leq C\Sigma_{a}^{\prime}(\\|f\\|_{\infty}+\delta)\\|f\\|_{\infty}\\|m_{2}-m_{1}\\|_{L^{2}((0,T)\times\Omega)}$		(33)

for some constant $C>0$ . Using the averaging lemma [27], we obtain that there exists a constant $C^{\prime}>0$ that

	$\displaystyle\\|Fm_{1}-Fm_{2}\\|_{H^{1/2}((0,T)\times\Omega)}$	$\displaystyle=\\|\langle{w}\rangle\\|_{H^{1/2}((0,T)\times\Omega)}\leq C^{\prime}\\|w\\|_{\mathcal{W}_{2}((0,T)\times\Omega\times\mathbb{S}^{d-1})}$		(34)
		$\displaystyle\leq CC^{\prime}\Sigma_{a}^{\prime}(\\|f\\|_{\infty}+\delta)\\|f\\|_{\infty}\\|m_{2}-m_{1}\\|_{L^{2}((0,T)\times\Omega)}.$		(34)

Combining this with the Kondrachov embedding theorem, we have shown that $F:\mathcal{M}\to\mathcal{M}$ is a continuous compact operator. By the Schauder fixed point theorem, there exists a fixed point for $\mathcal{M}$ , and hence (28) has a non-negative solution.

To prove the uniqueness of the solution, we use Lemma 3.1. We first observe that for any fixed point $m^{\ast}$ of $F$ , it must satisfy the conditions: $m^{\ast}\leq\|f\|_{\infty}$ and $m^{\ast}>0$ . This is due to the fact that $f$ is strictly positive. Hence there are no fixed points on the boundary $\partial\mathcal{M}$ . Next, we show that $F^{\prime}(m)$ cannot have $1$ as its eigenvalue. Let $\phi$ be the solution to (31) with $m\in\mathcal{M}$ and $\delta m$ be a perturbation that $m+\delta m\in\mathcal{M}$ , we denote by $w$ the solution to the following equation:

$\displaystyle\frac{\partial w}{\partial z}+\mathbf{k}\cdot\nabla w+\Sigma_{a}(m)w+\Sigma_{s}\left(w-\mathcal{K}w\right)$	$\displaystyle=-\Sigma_{a}^{\prime}(m)\delta m\phi\quad$	in	$\displaystyle X,$	(35)
$\displaystyle w(z,\mathbf{x},\mathbf{k})$	$\displaystyle=0\quad$	on	$\displaystyle\Gamma_{-},$
$\displaystyle w(0,\mathbf{x},\mathbf{k})$	$\displaystyle=0\quad$	on	$\displaystyle X_{0}\,.$

Then the Fréchet derivative $F^{\prime}(m)[\delta m]=\langle{w}\rangle$ . Suppose $F^{\prime}(m)$ indeed has $1$ as its eigenvalue and $\langle{w}\rangle$ as the corresponding nonzero eigenfunction. Then $F^{\prime}(\langle{w}\rangle)=\langle{w}\rangle$ and

$\displaystyle\frac{\partial w}{\partial z}+\mathbf{k}\cdot\nabla w+\Sigma_{a}(m)w+\Sigma_{s}\left(w-\mathcal{K}w\right)$	$\displaystyle=-\Sigma_{a}^{\prime}(m)\langle{w}\rangle\phi\quad$	in	$\displaystyle X,$	(36)
$\displaystyle w(z,\mathbf{x},\mathbf{k})$	$\displaystyle=0\quad$	on	$\displaystyle\Gamma_{-},$
$\displaystyle w(0,\mathbf{x},\mathbf{k})$	$\displaystyle=0\quad$	on	$\displaystyle X_{0}\,.$

Using the notations $\Sigma_{t}=\Sigma_{a}(m)+\Sigma_{s}$ and $R=\Sigma_{s}\mathcal{K}w-\Sigma_{a}^{\prime}(m)\langle{w}\rangle\phi$ , we can write the solution to (36) as

	$\displaystyle w(z,\mathbf{x},\mathbf{k})$	$\displaystyle=\int_{0}^{z\wedge\tau_{-}(\mathbf{x},\mathbf{k})}\exp^{-\int_{0}^{s}\Sigma_{t}(z-l,\mathbf{x}-l\mathbf{k})dl}R(z-s,\mathbf{x}-s\mathbf{k},\mathbf{k})ds$
		$\displaystyle=\int_{0}^{z\wedge\tau_{-}(\mathbf{x},\mathbf{k})}\Sigma_{t}(z-s,\mathbf{x}-s\mathbf{k})e^{-\int_{0}^{s}\Sigma_{t}(z-l,\mathbf{x}-l\mathbf{k})dl}\frac{R(z-s,\mathbf{x}-s\mathbf{k},\mathbf{k})}{\Sigma_{t}(z-s,\mathbf{x}-s\mathbf{k})}ds\,,$

which then gives the bound

$\displaystyle\|w(z,\mathbf{x},\mathbf{k})\|$	$\displaystyle\leq\int_{0}^{z\wedge\tau_{-}(\mathbf{x},\mathbf{k})}\Sigma_{t}(z-s,\mathbf{x}-s\mathbf{k})e^{-\int_{0}^{s}\Sigma_{t}(z-l,\mathbf{x}-l\mathbf{k})dl}\left\|\frac{R(z-s,\mathbf{x}-s\mathbf{k},\mathbf{k})}{\Sigma_{t}(z-s,\mathbf{x}-s\mathbf{k})}\right\|ds$	(37)
	$\displaystyle\leq\left(1-e^{-\int_{0}^{z\wedge\tau_{-}(\mathbf{x},\mathbf{k})}\Sigma_{t}(z-l,\mathbf{x}-l\mathbf{k})dl}\right)\sup_{0\leq s\leq z\wedge\tau_{-}(\mathbf{x},\mathbf{k})}\left\|\frac{R(z-s,\mathbf{x}-s\mathbf{k},\mathbf{k})}{\Sigma_{t}(z-s,\mathbf{x}-s\mathbf{k})}\right\|$
	$\displaystyle\leq\gamma\sup_{(0,T)\times\Omega\times\mathbb{S}^{d-1}}\left\|\frac{R(z,\mathbf{x},\mathbf{k})}{\Sigma_{t}(z,\mathbf{x})}\right\|$

for some $0<\gamma<1$ . Here $a\wedge b=\min(a,b)$ and $\tau_{-}(\mathbf{x},\mathbf{k})$ is the distance from $\mathbf{x}$ to $\Gamma_{-}$ in the direction of $-\mathbf{k}$ . The next step is to show that $\sup_{(0,T)\times\Omega\times\mathbb{S}^{d-1}}\left|\frac{R(z,\mathbf{x},\mathbf{k})}{\Sigma_{t}(z,\mathbf{x})}\right|\leq\|w\|_{\infty}$ , which, when combined with the bound in (37), leads to the bound $|w(z,\mathbf{x},\mathbf{k})|\leq\gamma\|w\|_{\infty}$ (and hence $w=0$ ). This contradicts the assumption that $\langle{w}\rangle$ is the eigenfunction corresponding to the eigenvalue $1$ of $F^{\prime}$ . We derive the bound under the two assumptions in the theorem.

(a)

When condition (i) is satisfied, we deduce from it that

\nu_{d-1}\Sigma_{a}^{\prime}(m)\|f\|_{\infty}\leq\Sigma_{a}(m)\,.

Meanwhile, we also have that

|R(z,\mathbf{x},\mathbf{k})|\leq\Sigma_{s}\|w\|_{\infty}+\nu_{d-1}\Sigma_{a}^{\prime}(m)\|f\|_{\infty}\|w\|_{\infty},

Combining the above two bounds gives us

|R(z,\mathbf{x},\mathbf{k})|\leq(\Sigma_{s}+\Sigma_{a}(m))\|w\|_{\infty}=\Sigma_{t}\|w\|_{\infty}.

Therefore, we have $\sup_{(0,T)\times\Omega\times\mathbb{S}^{d-1}}\left|\frac{R(z,\mathbf{x},\mathbf{k})}{\Sigma_{t}(z,\mathbf{x})}\right|\leq\|w\|_{\infty}$ .

(b)

For the case when condition (ii) is satisfied, we first observe that

R(z,\mathbf{x},\mathbf{k})=\Sigma_{s}\left(\mathcal{K}w-\underline{\theta}\langle{w}\rangle\right)+\left(\Sigma_{s}\underline{\theta}-\Sigma_{a}^{\prime}(m)\phi\right)\langle{w}\rangle

implies

\left|\frac{R(z,\mathbf{x},\mathbf{k})}{\Sigma_{t}(z,\mathbf{x})}\right|\leq\frac{\Sigma_{s}(1-\nu_{d-1}\underline{\theta})+\nu_{d-1}|\Sigma_{s}\underline{\theta}-\Sigma_{a}^{\prime}(m)\phi|}{\Sigma_{t}(z,\mathbf{x})}\|w\|_{\infty}.

(38)

When $\|f\|_{\infty}$ satisfies

\|f\|_{\infty}\leq\frac{1}{\nu_{d-1}}\frac{2\Sigma_{s}\underline{\theta}}{\Sigma_{a}^{\prime}(\|f\|_{\infty}+\delta)},\quad\forall(z,\mathbf{x})\in(0,T)\times\Omega,

(39)

we obtain that $\Sigma_{s}(1-\nu_{d-1}\underline{\theta})+\nu_{d-1}|\Sigma_{s}\underline{\theta}-\Sigma_{a}^{\prime}(m)\phi|\leq\Sigma_{s}$ which, when combined with (38), implies that

\left|\frac{R(z,\mathbf{x},\mathbf{k})}{\Sigma_{t}(z,\mathbf{x})}\right|\leq\|w\|_{\infty}.

Since $\delta>0$ is arbitrary, taking $\delta\to 0$ , the condition (39) becomes

\|f\|_{\infty}\leq\frac{1}{\nu_{d-1}}\frac{2\Sigma_{s}\underline{\theta}}{\Sigma_{a}^{\prime}(\|f\|_{\infty})},\quad\forall(z,\mathbf{x})\in(0,T)\times\Omega\,,

which is simply (ii).

The proof is now complete. ∎

To study the diffusion limit, we introduce a small parameter $\epsilon$ and denote by $W_{\epsilon}$ the solution to the scaled radiative transport equation:

$\displaystyle\frac{\partial W_{\epsilon}}{\partial z}+\frac{1}{\epsilon}\mathbf{k}\cdot\nabla W_{\epsilon}+\Sigma_{a}(\langle{W_{\epsilon}}\rangle)W_{\epsilon}+\frac{1}{\epsilon^{2}}\Sigma_{s}\left(W_{\epsilon}-\mathcal{K}W_{\epsilon}\right)$	$\displaystyle=0\quad$	in	$\displaystyle X,$	(40)
$\displaystyle W_{\epsilon}(z,\mathbf{x},\mathbf{k})$	$\displaystyle=0\quad$	on	$\displaystyle\Gamma_{-},$
$\displaystyle W_{\epsilon}(0,\mathbf{x},\mathbf{k})$	$\displaystyle=f(\mathbf{x})\quad$	on	$\displaystyle X_{0}\,.$

We have the following corollary using condition (ii) in Theorem 3.2.

Corollary 3.3.

If $\epsilon$ is sufficiently small such that

\Sigma_{a}^{\prime}(\|f\|_{\infty})\|f\|_{\infty}\leq\frac{1}{\nu_{d-1}}\frac{2\underline{\theta}\Sigma_{s}}{\epsilon^{2}},

then (40) admits a unique non-negative solution in $\mathcal{W}_{\infty}((0,T)\times\Omega\times\mathbb{S}^{d-1})$ . Moreover, the solution satisfies $\|W_{\epsilon}\|_{\infty}\leq\|f\|_{\infty}$ .

3.1 Asymptotic expansion

Let $\epsilon$ be sufficiently small such that (40) admits a unique solution. We formally expand the solution in powers of $\epsilon$ :

W_{\epsilon}(z,\mathbf{x},\mathbf{v})=W_{0}(z,\mathbf{x},\mathbf{v})+\epsilon W_{1}(z,\mathbf{x},\mathbf{v})+\epsilon^{2}W_{2}(z,\mathbf{x},\mathbf{v})+\phi_{\epsilon}(z,\mathbf{x},\mathbf{v})\,.

(41)

Substituting this into (40) and matching the equations at orders $\epsilon^{-2}$ , $\epsilon^{-1}$ and $\epsilon^{0}$ gives the system:

$\displaystyle(\mathcal{I}-\mathcal{K})W_{0}$	$\displaystyle=0,$	(42)
$\displaystyle\Sigma_{s}(\mathcal{I}-\mathcal{K})W_{1}+\mathbf{k}\cdot\nabla W_{0}$	$\displaystyle=0,$
$\displaystyle\frac{\partial W_{0}}{\partial z}+\Sigma_{s}(\mathcal{I}-\mathcal{K})W_{2}+\mathbf{k}\cdot\nabla W_{1}+\Sigma_{a}({W_{0}})W_{0}$	$\displaystyle=0\,.$

Following standard procedures [15, 19], we obtain from the first two equations that

W_{0}(z,\mathbf{x},\mathbf{k})=W_{0}(z,\mathbf{x}),\quad W_{1}(z,\mathbf{x},\mathbf{k})=-\sum_{i=1}^{d}\frac{D_{i}(\mathbf{k})}{\Sigma_{s}}\partial_{x_{i}}W_{0}(z,\mathbf{x})\,,

(43)

where $D_{i}(v)\in L^{\infty}(\mathbb{S}^{d-1})$ are the unique solutions to

(\mathcal{I}-\mathcal{K})D_{i}(\mathbf{k})=\mathbf{k}\cdot\mathbf{e}_{i},\quad\int_{\mathbb{S}^{d-1}}D_{i}(\mathbf{k})d\mathbf{k}=0\,,\quad i=1,2,\dots,d.

(44)

Next, we integrate the third equation in (42) over $\mathbb{S}^{d-1}$ , and utilize the fact that

\langle{(\mathcal{I}-\mathcal{K})W_{2}}\rangle=0

to get the equation for $W_{0}$ . This leads to

\int_{\mathbb{S}^{d-1}}\left(\frac{\partial W_{0}}{\partial z}+\left(\mathbf{k}\cdot\nabla W_{1}+\Sigma_{a}(W_{0})W_{0}\right)\right)d\mathbf{k}=0\,.

Since $W_{0}$ is independent of $\mathbf{k}$ , we have that $W_{0}$ solves the diffusion equation:

$\displaystyle\frac{\partial W_{0}}{\partial z}-\nabla\cdot\left(\frac{A}{\Sigma_{s}}\nabla W_{0}\right)+\Sigma_{a}(W_{0})W_{0}$	$\displaystyle=0,\quad$	in	$\displaystyle(0,T)\times\Omega,$	(45)
$\displaystyle W_{0}(z,\mathbf{x})$	$\displaystyle=0,\quad$	on	$\displaystyle(0,T)\times\partial\Omega,$
$\displaystyle W_{0}(0,\mathbf{x})$	$\displaystyle=f(\mathbf{x})\quad$	on	$\displaystyle\Omega,$

where the matrix $A$ is positive definite and is defined as

A_{ij}=\int_{\mathbb{S}^{d-1}}\mathbf{k}\cdot\mathbf{e}_{i}D_{j}(\mathbf{k})d\mathbf{k}\,.

Under our assumption that the scattering phase function $p(\mathbf{k}\cdot\mathbf{k}^{\prime})$ is rotation invariant, we have that $A_{ij}=\frac{1}{d(1-g)}\delta_{ij}$ , where $g\in(-1,1)$ is the anisotropy parameter.

Let $\phi_{\epsilon}$ be the remainder in the expansion. Then we can check that $\phi_{\epsilon}$ satisfies the equation:

$\displaystyle\partial_{z}\phi_{\epsilon}+\frac{1}{\epsilon}\mathbf{k}\cdot\nabla\phi_{\epsilon}+\frac{\Sigma_{s}}{\epsilon^{2}}\left(\mathcal{I}-\mathcal{K}\right)\phi_{\epsilon}+\Sigma_{a}(\langle{W_{\epsilon}}\rangle)\phi_{\epsilon}$	$\displaystyle=\epsilon h_{1},\;$	in	$\displaystyle X$	(46)
$\displaystyle\phi_{\epsilon}(z,\mathbf{x},\mathbf{k})$	$\displaystyle=\epsilon h_{2},\quad$	on	$\displaystyle\Gamma_{-},$
$\displaystyle\phi_{\epsilon}(0,\mathbf{x},\mathbf{k})$	$\displaystyle=\epsilon h_{3},\quad$	on	$\displaystyle X_{0}\,,$

where the functions $h_{1}$ , $h_{2}$ , and $h_{3}$ are respectively of the forms:

	$\displaystyle h_{1}(z,\mathbf{x},\mathbf{k})$	$\displaystyle=-\partial_{z}W_{1}-\mathbf{k}\cdot\nabla W_{2}-\Sigma_{a}(\langle{W_{\epsilon}}\rangle)(W_{1}+\epsilon W_{2})$
		$\displaystyle\quad-\epsilon\partial_{z}W_{2}+\frac{1}{\epsilon}\left[\Sigma_{a}(W_{0})-\Sigma_{a}(\langle{W_{\epsilon}}\rangle)\right]W_{0},$
	$\displaystyle h_{2}(z,\mathbf{x},\mathbf{k})$	$\displaystyle=-\sum_{i=1}^{d}\frac{D_{i}(\mathbf{k})}{\Sigma_{s}}\frac{\partial W_{0}}{\partial x_{i}}-\epsilon W_{2},$
	$\displaystyle h_{3}(0,\mathbf{x},\mathbf{k})$	$\displaystyle=-W_{1}-\epsilon W_{2}\,.$

Let $\delta:=\frac{1}{2}\inf_{[0,T]\times\Omega}\Sigma_{a,0}(z,\mathbf{x})$ . Then by the assumtion in (30), $\delta>0$ . We then have that for any $z\in[0,T)$ ,

	$\displaystyle\\|\phi_{\epsilon}(z,\cdot,\cdot)\\|_{\infty}$	$\displaystyle\leq\epsilon\\|h_{3}\\|_{\infty}e^{-\delta z}+\epsilon\int_{0}^{z}e^{-\delta(z-s)}(\\|h_{1}\\|_{\infty}+\delta\\|h_{2}\\|_{\infty})ds$		(47)
		$\displaystyle\leq\epsilon T(\\|h_{1}\\|_{\infty}+\delta\\|h_{2}\\|_{\infty})+\epsilon\\|h_{3}\\|_{\infty}.$		(47)

It remains to show that $h_{1}$ , $h_{2}$ , and $h_{3}$ are bounded. We first observe from the equations in (42) that

\|W_{1}\|_{\infty}\leq C_{1}\|W_{0}\|_{C([0,T),C^{1}(\Omega))},\quad\|W_{2}\|_{\infty}\leq C_{2}\|W_{0}\|_{C([0,T),C^{2}(\Omega))}

We then take the derivative with respect to $z$ of the equations in (42) to deduce that

\|\partial_{z}W_{1}\|_{\infty}\leq C_{3}\|W_{0}\|_{C([0,T),C^{3}(\Omega))},\quad\|\partial_{z}W_{2}\|_{\infty}\leq C_{4}\|W_{0}\|_{C([0,T),C^{4}(\Omega))}

together with

\|\mathbf{k}\cdot\nabla W_{2}\|_{\infty}\leq C_{5}\|W_{0}\|_{C([0,T),C^{3}(\Omega))}\,.

Therefore, given that $W_{0}\in C([0,T),C^{4}(\Omega))$ , we have the following estimate for the diffusion approximation:

\|W_{\epsilon}-W_{0}\|_{\infty}\leq\|\phi_{\epsilon}\|_{\infty}+\epsilon\|W_{1}\|_{\infty}+\epsilon^{2}\|W_{2}\|_{\infty}=C(T,\|W_{0}\|_{C([0,T),C^{4}(\Omega))})\epsilon.

To ensure the regularity of the solution $W_{0}$ , at least for a short time, we simply need the initial condition $f$ to be smooth enough since $W_{0}$ would be smoother than the initial condition due to the diffusive nature.

We can prove the following result.

Theorem 3.4.

Assume that $f\in C^{4}_{0}(\Omega)$ is non-negative, and the absorption and scattering coefficients satisfy the condition

\Sigma_{a,l}\in C^{2}(\Omega),\quad\Sigma_{a,l}\geq 0,\quad\Sigma_{s}>0.

Then the diffusion equation (45) admits a unique strong solution $W_{0}\in C([0,T),C^{4}(\Omega))$ when $1\leq d\leq 3$ .

Proof.

Let $L=-\nabla\cdot\left(\frac{A}{\Sigma_{s}}\nabla\right)$ . Then $-L$ is the infinitesimal generator of an analytic semigroup $G(t)$ on $L^{2}(\Omega)$ and $\|G(t)\|\leq 1$ for all $t\geq 0$ . We denote

\mathcal{D}(L)=H^{2}(\Omega)\cap H^{1}_{0}(\Omega).

By [43, Theorem 8.4.4] and [43, Theorem 6.3.1], we have that when $f\in\mathcal{D}(L)$ , there exists a unique local strong solution $W_{0}$ to the equation (45) on $(0,T)\times\Omega$ , that is,

W_{0}\in C([0,T),L^{2}(\Omega))\cap C((0,T),H^{2}(\Omega)\cap H^{1}_{0}(\Omega))\cap C^{1}((0,T),L^{2}(\Omega))\,.

(48)

Moreover, $0\leq W_{0}\leq\|f\|_{\infty}$ by the comparison principle. Then the result of [43, Corollary 6.3.2] ensures that $\partial_{z}W_{0}$ is locally Hölder continuous for $z\in(0,T)$ . Hence $W_{0}(z,\mathbf{x})$ and $\partial_{z}W_{0}(z,\mathbf{x})$ are both continuous on $(0,T)\times\overline{\Omega}$ . This means that $\Sigma_{a}(W_{0})$ is also continuous. Therefore we must have $W_{0}(z,\cdot)\in C^{2}(\Omega)$ , which means $W_{0}$ is a classical solution.

Let $g(\mathbf{x}):=-Lf-\Sigma_{a}(f)$ , $f\in C^{2}(\Omega)\in\mathcal{D}(L)$ . By differentiating (45), we find that $\psi:=\partial_{z}W_{0}$ satisfies the following equation:

$\displaystyle\partial_{z}\psi+L\psi+\left(\Sigma_{a}^{\prime}(W_{0})W_{0}+\Sigma_{a}(W_{0})\right)\psi$	$\displaystyle=0\quad$	in	$\displaystyle(0,T)\times\Omega,$	(49)
$\displaystyle\psi(z,\mathbf{x})$	$\displaystyle=0\quad$	on	$\displaystyle(0,T)\times\partial\Omega,$
$\displaystyle\psi(0,\mathbf{x})$	$\displaystyle=g(\mathbf{x})\quad$	on	$\displaystyle\{0\}\times\Omega.$

Following a similar process, we can deduce that $\psi(z,\cdot)\in C^{2}(\Omega)$ . Since $W_{0}(z,\cdot)\in C^{2}(\Omega)$ and $\partial_{z}W_{0}(z,\cdot)\in C^{2}(\Omega)$ for $z\in(0,T)$ , we have $\partial_{z}W_{0}(z,\cdot)+\Sigma_{a}(W_{0}(z,\cdot))W_{0}(z,\cdot)\in C^{2}(\Omega)$ . By classical regularity theory for elliptic equations [26], $W_{0}(z,\cdot)\in C^{4}(\Omega)$ . ∎

Remark 3.5.

We have assumed so far that the initial condition $f$ is independent of the variable $\mathbf{k}$ . In fact, the case of $f$ depending on $\mathbf{k}$ , that is, $f=f(\mathbf{x},\mathbf{k})$ , can be treated in a similar manner by introducing another fast variable $\theta=\frac{z}{\epsilon^{2}}$ , as in [19, Section XXI.5.3]. We will not repeat the calculations here.

3.2 The case of degenerate coefficients

Let us now briefly consider the case when the problem is degenerate, that is, when the absorption coefficient can vanish in part of the domain of interest. More precisely, we relax the requirement that all $\Sigma_{a,l}>0$ to the following:

\Sigma_{a,l}\geq 0,\quad\Sigma_{a}(\|f\|_{\infty})>0,\quad\forall(z,\mathbf{x})\in[0,T]\times\Omega.

(50)

In this case, $\Sigma_{a}^{\prime}(\|f\|_{\infty})>0$ . When $\epsilon$ is sufficiently small, the scaled transport equation (40) admits a unique solution in $L^{\infty}(X)$ . Let $w_{\epsilon}$ be the solution to the following linear transport equation:

$\displaystyle\partial_{z}w_{\epsilon}+\frac{1}{\epsilon}\mathbf{k}\cdot\nabla w_{\epsilon}+\Sigma_{a}(\\|f\\|_{\infty})w_{\epsilon}+\frac{\Sigma_{s}}{\epsilon^{2}}(\mathcal{I}-\mathcal{K})w_{\epsilon}$	$\displaystyle=0\quad$	in	$\displaystyle X,$	(51)
$\displaystyle w_{\epsilon}(z,\mathbf{x},\mathbf{k})$	$\displaystyle=0\quad$	on	$\displaystyle\Gamma_{-},$
$\displaystyle w_{\epsilon}(0,\mathbf{x},\mathbf{k})$	$\displaystyle=f(\mathbf{x})\quad$	on	$\displaystyle X_{0}.$

Since the absorption coefficient $\Sigma_{a}(\|f\|_{\infty})\geq\Sigma_{a}(\langle{W_{\epsilon}}\rangle)$ , we conclude that $w_{\epsilon}\leq W_{\epsilon}$ for any $\epsilon>0$ . On the other hand, as $\epsilon\to 0$ , we have $w_{\epsilon}\to w_{0}$ , where $w_{0}$ is the solution to

$\displaystyle\partial_{z}w_{0}-\nabla\cdot\left(\frac{A}{\Sigma_{s}}\nabla w_{0}\right)+\Sigma_{a}(\\|f\\|_{\infty})w_{0}$	$\displaystyle=0\quad$	in	$\displaystyle(0,T)\times\Omega,$	(52)
$\displaystyle w_{0}(z,\mathbf{x})$	$\displaystyle=0\quad$	on	$\displaystyle(0,T)\times\partial\Omega,$
$\displaystyle w_{0}(0,\mathbf{x})$	$\displaystyle=f(\mathbf{x})\quad$	on	$\displaystyle\Omega.$

Hence $\|W_{\epsilon}\|\geq\|w_{\epsilon}\|\geq\inf_{[0,T]\times\Omega}w_{0}-C\epsilon$ for some $C>0$ . Because $\inf_{(0,T)\times\Omega}w_{0}>0$ strictly for $\epsilon$ sufficiently small, we conclude that $W_{\epsilon}$ is bounded from below by a positive number, which implies $\Sigma_{a}(\langle{W_{\epsilon}}\rangle)$ is also strictly positive. Then repeat the process in (47) by setting $\delta=\frac{1}{2}\inf_{(0,T)\times\Omega}\Sigma_{a}(\langle{W_{\epsilon}}\rangle)$ instead, we obtain the same conclusion that $\|W_{\epsilon}-W_{0}\|_{\infty}\leq C^{\prime}\epsilon$ for a constant $C^{\prime}$ independent of $\epsilon$ .

4 An application in inverse problems

We now consider the following inverse medium problem as a direct application of the transport model:

$\displaystyle\partial_{z}u+\mathbf{k}\cdot\nabla u+\Sigma_{a}(\langle{u}\rangle)u+\Sigma_{s}(\mathcal{I}-\mathcal{K})u$	$\displaystyle=0\quad$	in	$\displaystyle X,$	(53)
$\displaystyle u(z,\mathbf{x},\mathbf{k})$	$\displaystyle=0\quad$	on	$\displaystyle\Gamma_{-},$
$\displaystyle u(0,\mathbf{x},\mathbf{k})$	$\displaystyle=f(\mathbf{x})\quad$	on	$\displaystyle X_{0},$

where $f(\mathbf{x})\in L^{\infty}(X_{0})$ is a strictly positive source function. We assume that (53) has a unique positive solution $u\in\mathcal{W}_{\infty}$ .

We assume that the absorption coefficient $\Sigma_{a}$ is not known, but we have additional data that is the density of the solution, that is,

g(z,\mathbf{x})=\langle{u}\rangle:=\int_{\mathbb{S}^{d-1}}u(z,\mathbf{x},\mathbf{k})d\mathbf{k}\,.

(54)

The inverse problem amounts to finding the unknown absorption coefficients $\Sigma_{a,l}$ from the observed data $g$ from a given $f$ .

We can prove the following result.

Theorem 4.1.

Let $g$ and $\widetilde{g}$ be data of the form (54) generated from (53) with coefficients $\Sigma_{a}$ and $\widetilde{\Sigma}_{a}$ respectively. Then $g=\widetilde{g}$ implies $\Sigma_{a}(\langle{u}\rangle)=\widetilde{\Sigma}_{a}(\langle{\widetilde{u}}\rangle)$ .

Proof.

Let $\delta u=u-\widetilde{u}$ . We verify that for any $\delta u$ , we have the identity

\int_{\Omega\times\mathbb{S}^{d-1}}\left[\mathbf{k}\cdot\nabla\delta u\right]\frac{\delta u}{\tilde{u}}d\mathbf{x}d\mathbf{k}=\int_{\Omega\times\mathbb{S}^{d-1}}\mathbf{k}\cdot\nabla\frac{|\delta u|^{2}}{2\tilde{u}}d\mathbf{x}d\mathbf{k}-\int_{\Omega\times\mathbb{S}^{d-1}}\left[\mathbf{k}\cdot\nabla\frac{1}{\tilde{u}}\right]\frac{|\delta u|^{2}}{2}d\mathbf{x}d\mathbf{k}\,,

(55)

and the identity

\mathbf{k}\cdot\nabla\frac{1}{\tilde{u}}=\frac{1}{\tilde{u}^{2}}\partial_{z}\tilde{u}+\frac{\widetilde{\Sigma_{a}}(\langle{\tilde{u}}\rangle)}{\tilde{u}}+\frac{\Sigma_{s}(\mathcal{I}-\mathcal{K})\tilde{u}}{|\tilde{u}|^{2}}\,.

(56)

Using the fact that $g=\widetilde{g}$ , we can also conclude that

\langle{\delta u}\rangle=0\,.

(57)

It is also straightforward to check that $\delta u$ solves the transport equation:

\partial_{z}\delta u+\mathbf{k}\cdot\nabla\delta u+{\Sigma_{a}}(\langle{{u}}\rangle)\delta u+\Sigma_{s}(\mathcal{I}-\mathcal{K})\delta u=(\widetilde{\Sigma}_{a}(\langle{\tilde{u}}\rangle)-\Sigma_{a}(\langle{{u}}\rangle))\tilde{u}\,,

(58)

with zero initial and incoming boundary conditions. We multiply this equation by $\frac{\delta u}{\tilde{u}}$ and integrate over $\Omega\times\mathbb{S}^{d-1}$ to obtain

		$\displaystyle\int_{\Omega\times\mathbb{S}^{d-1}}(\partial_{z}\delta u)\frac{\delta u}{\tilde{u}}d\mathbf{x}d\mathbf{k}+\int_{\Omega\times\mathbb{S}^{d-1}}\mathbf{k}\cdot\nabla\frac{\|\delta u\|^{2}}{2\tilde{u}}d\mathbf{x}d\mathbf{k}-\int_{\Omega\times\mathbb{S}^{d-1}}\frac{\|\delta u\|^{2}}{2\|\tilde{u}\|^{2}}\partial_{z}\tilde{u}d\mathbf{x}d\mathbf{k}$		(59)
		$\displaystyle-\int_{\Omega\times\mathbb{S}^{d-1}}\frac{\widetilde{\Sigma_{a}}(\langle{\tilde{u}}\rangle)+\Sigma_{s}}{2\tilde{u}}\|\delta u\|^{2}d\mathbf{x}d\mathbf{k}+\int_{\Omega\times\mathbb{S}^{d-1}}\frac{(\Sigma_{s}\mathcal{K}\tilde{u})\|\delta u\|^{2}}{\|\tilde{u}\|^{2}}d\mathbf{x}d\mathbf{k}$
		$\displaystyle+\int_{\Omega\times\mathbb{S}^{d-1}}\frac{\Sigma_{a}(\langle{u}\rangle)}{\tilde{u}}\|\delta u\|^{2}d\mathbf{x}d\mathbf{k}+\Sigma_{s}\int_{\Omega\times\mathbb{S}^{d-1}}\left[(\mathcal{I}-\mathcal{K})\delta u\right]\frac{\delta u}{\tilde{u}}d\mathbf{x}d\mathbf{k}$
		$\displaystyle=\int_{\Omega\times\mathbb{S}^{d-1}}(\widetilde{\Sigma_{a}}(\langle{\tilde{u}}\rangle)-\Sigma_{a}(\langle{{u}}\rangle))\delta ud\mathbf{x}d\mathbf{k}.$

where we have used the identities (55) and (56).

We first observe that since $\Sigma_{a}(\langle{{u}}\rangle)$ and $\widetilde{\Sigma_{a}}(\langle{\tilde{u}}\rangle)$ do not depend on $\mathbf{k}$ , the right hand side of (59) vanishes due to (57).

To handle the left hand side of (59), we observe that

\int_{\Omega\times\mathbb{S}^{d-1}}\mathbf{k}\cdot\nabla\frac{|\delta u|^{2}}{2\tilde{u}}d\mathbf{x}d\mathbf{k}=\int_{\Gamma_{+}}\mathbf{k}\cdot\mathbf{n}\frac{|\delta u|^{2}}{2\tilde{u}}\geq 0\,,

(60)

(\partial_{z}\delta u)\frac{\delta u}{\tilde{u}}=\frac{1}{2}\partial_{z}\left[\frac{|\delta u|^{2}}{\tilde{u}}\right]+\frac{1}{2}\frac{|\delta u|^{2}}{|\tilde{u}|^{2}}\partial_{z}\tilde{u}\,,

(61)

and

\Sigma_{s}\int_{\Omega\times\mathbb{S}^{d-1}}\left[(\mathcal{I}-\mathcal{K})\delta u\right]\frac{\delta u}{\tilde{u}}d\mathbf{x}d\mathbf{k}=\Sigma_{s}\int_{\Omega\times\mathbb{S}^{d-1}}\frac{|\delta u|^{2}}{\tilde{u}}d\mathbf{x}d\mathbf{k}-\Sigma_{s}\int_{\Omega\times\mathbb{S}^{d-1}}(\mathcal{K}\delta u)\frac{\delta u}{\tilde{u}}d\mathbf{x}d\mathbf{k}\,.

(62)

We can also prove the following inequality (see Appendix A),

\int_{\Omega\times\mathbb{S}^{d-1}}(\mathcal{K}\delta u)\frac{\delta u}{\tilde{u}}d\mathbf{x}d\mathbf{k}\leq\int_{\Omega\times\mathbb{S}^{d-1}}(\mathcal{K}\tilde{u})\frac{|\delta u|^{2}}{|\tilde{u}|^{2}}d\mathbf{x}d\mathbf{k}+\nu_{d-1}\frac{\kappa^{2}}{4}\int_{\Omega\times\mathbb{S}^{d-1}}\frac{|\delta u|^{2}}{\tilde{u}}d\mathbf{x}d\mathbf{k},

(63)

where $\kappa=\frac{(\overline{\theta}-\underline{\theta})}{2\sqrt{\underline{\theta}}}$ , the constants $\overline{\theta}$ and $\underline{\theta}$ are defined in (29).

Let $M(z,\mathbf{x}):=\Sigma_{a}(\langle{u}\rangle)-\frac{\widetilde{\Sigma_{a}}(\langle{\tilde{u}}\rangle)}{2}+\Sigma_{s}\left(\frac{2-\nu_{d-1}\kappa^{2}}{4}\right)$ . We can then deduce from (59), using (60), (61), (62), and (63), that

\displaystyle\frac{1}{2}\partial_{z}\int_{\Omega\times\mathbb{S}^{d-1}}\frac{|\delta u|^{2}}{\tilde{u}}d\mathbf{x}d\mathbf{k}+\int_{\Omega\times\mathbb{S}^{d-1}}M(z,\mathbf{x})\frac{|\delta u|^{2}}{\tilde{u}}d\mathbf{x}d\mathbf{k}\leq 0.

(64)

Since the coefficients $\Sigma_{a,l}$ and $\Sigma_{s}$ are finite and both $\langle{u}\rangle,\langle{\tilde{u}}\rangle$ are bounded by $\|f\|_{L^{\infty}(X_{0})}$ , there exists a constant $\underline{M}=\inf_{(0,T)\times\Omega}M(z,\mathbf{x})$ that

\frac{1}{2}\partial_{z}\int_{\Omega\times\mathbb{S}^{d-1}}\frac{|\delta u|^{2}}{\tilde{u}}d\mathbf{x}d\mathbf{k}+\underline{M}\int_{\Omega\times\mathbb{S}^{d-1}}\frac{|\delta u|^{2}}{\tilde{u}}d\mathbf{x}d\mathbf{k}\leq 0.

Then by the Grönwall inequality and the initial condition $\delta u=0$ at $z=0$ , we must have that $\delta u\equiv 0$ . Therefore $u=\tilde{u}$ , which implies that

\Sigma_{a}(\langle{u}\rangle)=-\frac{\partial_{z}u+\mathbf{k}\cdot\nabla u+\Sigma_{s}(\mathcal{I}-\mathcal{K})u}{u}=-\frac{\partial_{z}\tilde{u}+\mathbf{k}\cdot\nabla\tilde{u}+\Sigma_{s}(\mathcal{I}-\mathcal{K})\tilde{u}}{\tilde{u}}=\widetilde{\Sigma_{a}}(\langle{\tilde{u}}\rangle)\,.

(65)

The proof is complete. ∎

The following corollary is a direct result of the comparison principle and Theorem 4.1.

Corollary 4.2.

Under the assumption of Theorem 4.1, the coefficients $\Sigma_{a,l}$ can be uniquely determined with finitely many data sets $\langle{u_{j}}\rangle$ , $j=1,2,\dots,L+1$ if the initial conditions satisfy $0<f_{1}<f_{2}<\dots<f_{L+1}$ on $\Omega$ .

5 Concluding remarks

This work describes the derivation of semilinear radiative transport models for wave propagation in highly-scattering media with nonlinear absorption. While the technical aspects of the derivation are relatively standard, we believe that our work provides a theoretical justification for the semilinear radiative transport models, as well as their diffusion approximations, used in applications such as multi-photon imaging [10, 36, 44, 45, 51, 55, 56, 57].

As we have remarked before, one concrete example of the quadratic absorption we considered here is two-photon absorption in nonlinear optics [38, 48]. The radiative transport equation we derive for this case, equation (20), is different from the two-photon radiative transport equation of [40] where the phase space intensity corresponds to a two-photon entangled state of light, not two-photon absorption.

The calculation we have presented here does not generalize to media with nonlinearities such as those that arise in the Kerr effect and second harmonic generation [3, 29, 46, 52]. The derivation of transport equations for wave propagation in such media is a topic of great interest, but is much more challenging due to the richness of the behavior of the corresponding wave equations [17, 18, 50]. We point to the derivation in [22] in the context of the nonlinear Schrödinger equation within the mean-field approximation, and leave further investigations in this direction to future work. We note that the acousto-optic effect, in which light undergoes a frequency shift due to interaction with an acoustic wave, has also been studied in random media. Although this effect is not nonlinear in the sense considered in this paper, it is possible to develop a suitable kinetic model and associated radiative transport equations [31].

Acknowledgment

The work of JK is supported in part by the Simons Foundation Math + X Investigator Award #376319 (to Michael I. Weinstein). KR is partially supported by the National Science Foundation through grants DMS-1913309 and DMS-1937254. JCS acknowledges support from the NSF Grant DMS-1912821 and the AFOSR Grant FA9550-19-1-0320.

Appendix A Proof of inequality (63)

Proof.

Using the AM-GM inequality, we deduce that

\int_{\Omega\times\mathbb{S}^{d-1}}(\mathcal{K}\delta u)\frac{\delta u}{\tilde{u}}d\mathbf{x}d\mathbf{k}\leq\int_{\Omega\times\mathbb{S}^{d-1}}(\mathcal{K}\tilde{u})\frac{|\delta u|^{2}}{|\tilde{u}|^{2}}d\mathbf{x}d\mathbf{k}+\frac{1}{4}\int_{\Omega\times\mathbb{S}^{d-1}}\left|\frac{\mathcal{K}\delta u}{\sqrt{\mathcal{K}\tilde{u}}}\right|^{2}d\mathbf{x}d\mathbf{k}.

Since $\underline{\theta}\leq p(\mathbf{k},\mathbf{k}^{\prime})\leq\overline{\theta}$ and $\langle{\delta u}\rangle=0$ , we have the preliminary estimate

\left|\frac{\mathcal{K}\delta u}{\sqrt{\mathcal{K}\tilde{u}}}\right|=\left|\frac{\mathcal{K}\delta u-\frac{\overline{\theta}+\underline{\theta}}{2}\langle{\delta u}\rangle}{\sqrt{\mathcal{K}\tilde{u}}}\right|\leq\frac{(\overline{\theta}-\underline{\theta})\langle{|\delta u|}\rangle}{2\sqrt{\underline{\theta}\langle{\tilde{u}}\rangle}}.

Denote by $\kappa$ the constant $\kappa:=\frac{(\overline{\theta}-\underline{\theta})}{2\sqrt{\underline{\theta}}}$ . We then have

\int_{\Omega\times\mathbb{S}^{d-1}}\left|\frac{\mathcal{K}\delta u}{\sqrt{\mathcal{K}\tilde{u}}}\right|^{2}d\mathbf{x}d\mathbf{k}\leq\kappa^{2}\int_{\Omega\times\mathbb{S}^{d-1}}\frac{\langle{|\delta u|}\rangle^{2}}{{\langle{\tilde{u}}\rangle}}d\mathbf{x}d\mathbf{k}=\nu_{d-1}\kappa^{2}\int_{\Omega}\frac{\langle{|\delta u|}\rangle^{2}}{{\langle{\tilde{u}}\rangle}}d\mathbf{x}.

Using the Cauchy-Schwartz inequality, we arrive at

\int_{\mathbb{S}^{d-1}}\frac{|\delta u(\mathbf{x},\mathbf{k})|^{2}}{\tilde{u}(\mathbf{x},\mathbf{k})}d\mathbf{k}\int_{\mathbb{S}^{d-1}}\tilde{u}(\mathbf{x},\mathbf{k})d\mathbf{k}\geq\left(\int_{\mathbb{S}^{d-1}}|\delta u(\mathbf{x},\mathbf{k})|d\mathbf{k}\right)^{2}=\langle{|\delta u|}\rangle^{2}\,.

Therefore, we have

\frac{\langle{|\delta u|}\rangle^{2}}{{\langle{\tilde{u}}\rangle}}\leq\int_{\mathbb{S}^{d-1}}\frac{|\delta u(\mathbf{x},\mathbf{k})|^{2}}{\tilde{u}(\mathbf{x},\mathbf{k})}d\mathbf{k}.

This implies that

\int_{\Omega\times\mathbb{S}^{d-1}}(\mathcal{K}\delta u)\frac{\delta u}{\tilde{u}}d\mathbf{x}d\mathbf{k}\leq\int_{\Omega\times\mathbb{S}^{d-1}}(\mathcal{K}\tilde{u})\frac{|\delta u|^{2}}{|\tilde{u}|^{2}}d\mathbf{x}d\mathbf{k}+\nu_{d-1}\frac{\kappa^{2}}{4}\int_{\Omega\times\mathbb{S}^{d-1}}\frac{|\delta u|^{2}}{\tilde{u}}d\mathbf{x}d\mathbf{k}.

This completes the proof. ∎

References

[1] J.-L. Akian and E. Savin, Kinetic modeling of multiple scattering of acoustic waves in randomly heterogeneous flows, Multiscale Modeling & Simulation, 19 (2021), pp. 1394–1424.
[2] H. Ammari, E. Bossy, J. Garnier, W. Jing, and L. Seppecher, Radiative transfer and diffusion limits for wave field correlations in locally shifted random media, J. Math. Phys., 54 (2013). 021501.
[3] Y. M. Assylbekov and T. Zhou, Direct and inverse problems for the nonlinear time-harmonic Maxwell equations in Kerr-type media, J. Spectr. Theor., (2020).
[4] G. Bal, Kinetics of scalar wave fields in random media, Wave Motion, 43 (2005), pp. 132–157.
[5] G. Bal, T. Komorowski, and L. Ryzhik, Kinetic limits for waves in a random medium, Kinet. Relat. Models, 3 (2010), pp. 529–644.
[6] G. Bal, G. Papanicolaou, and L. Ryzhik, Radiative transport limit for the random Schrödinger equation, Nonlinearity, 15 (2002), pp. 513–529.
[7] G. Bal and O. Pinaud, Accuracy of transport models for waves in random media, Wave Motion, 43 (2006).
[8] , Kinetic models for imaging in random media, SIAM Multiscale Model. Simul., 6 (2007), pp. 792–819.
[9] G. Bal and K. Ren, Transport-based imaging in random media, SIAM J. Appl. Math., 68 (2008), pp. 1738–1762.
[10] P. Bardsley, K. Ren, and R. Zhang, Quantitative photoacoustic imaging of two-photon absorption, J. Biomed. Opt., 23 (2018). 016002.
[11] I. Baydoun, É. Savin, R. Cottereau, D. Clouteau, and J. Guilleminot, Kinetic modeling of multiple scattering of elastic waves in heterogeneous anisotropic media, Wave Motion, 51 (2014), pp. 1325–1348.
[12] L. Borcea and J. Garnier, Derivation of a one-way radiative transfer equation in random media, Phys. Rev. E, 93 (2016). 022115.
[13] L. Borcea, J. Garnier, and K. Solna, Wave propagation and imaging in moving random media, Multiscale Modeling & Simulation, 17 (2019), pp. 31–67.
[14] R. W. Boyd, Nonlinear optics, Academic press, 2020.
[15] R. Carminati and J. C. Schotland, Principles of Scattering and Transport of Light, Cambridge University Press, 2021.
[16] S. Chen and Q. Li, Semiclassical limit of an inverse problem for the Schrödinger equation, Research in the Mathematical Sciences, 8 (2021), pp. 1–18.
[17] C. Conti, L. Angelani, and G. Ruocco, Light diffusion and localization in three-dimensional nonlinear disordered media, Phys. Rev. A, 75 (2007). 033812.
[18] C. Conti and L. Leuzzi, Complexity of waves in nonlinear disordered media, Phys. Rev. B, 83 (2011). 134204.
[19] R. Dautray and J.-L. Lions, Mathematical Analysis and Numerical Methods for Science and Technology, Vol VI, Springer-Verlag, Berlin, 1993.
[20] W. Denk, J. Strickler, and W. Webb, Two photon laser scanning fluorescence microscopy, Science, 248 (1990), pp. 73–76.
[21] L. Erdös and H. T. Yau, Linear Boltzmann equation as the weak coupling limit of a random Schrödinger Equation, Comm. Pure Appl. Math., 53(6) (2000), pp. 667–735.
[22] A. Fannjiang, G. Papanicolaou, and S. Jin, High frequency behavior of the focusing nonlinear Schrödinger equation with random inhomogeneities, SIAM J. Appl. Math., 63 (2003), pp. 1328–1358.
[23] A. Fannjiang and L. V. Ryzhik, Radiative transfer of sound waves in a random flow: Turbulent scattering, straining, and mode-coupling, SIAM J. Appl. Math., 61 (2001), pp. 1545–1577.
[24] A. C. Fannjiang, Radiative transfer limit of two-frequency Wigner distribution for random parabolic waves: an exact solution, Comptes Rendus Physique, 8 (2007), pp. 267–271.
[25] J. Garnier and K. Solna, Effective transport equations and enhanced backscattering in random waveguides, SIAM J. Appl. Math., 68 (2008), pp. 1574–1599.
[26] D. Gilbarg and N. S. Trudinger, Elliptic Partial Differential Equations of Second Order, Springer-Verlag, Berlin, 2000.
[27] F. Golse, P.-L. Lions, B. Perthame, and R. Sentis, Regularity of the moments of the solution of a transport equation, J. Func. Anal., 76 (1988), pp. 110–125.
[28] C. Gomez, Radiative transport limit for the random schrödinger equation with long-range correlations, J. Math. Pures Appl., 98 (2012), pp. 295–327.
[29] A. Goy and D. Psaltis, Digital reverse propagation in focusing Kerr media, Phys. Rev. A, 83 (2011). 031802.
[30] M. Guo and X.-P. Wang, Transport equations for a general class of evolution equations with random perturbations, J. Math. Phys., 40 (1999), pp. 4828–4848.
[31] J. Hoskins and J. Schotland, Acousto-optic effect in random media, Phys. Rev. E, 95 (2017). 033002.
[32] A. Ishimaru, Wave Propagation and Scattering in Random Media, IEEE Press, 1997.
[33] Y. Jing and N. Xiang, One-dimensional transport equation models for sound energy propagation in long spaces: Simulations and experiments, J. Acoust. Soc. Am., 127 (2010), pp. 2323–2331.
[34] R. Kellogg, Uniqueness in the schauder fixed point theorem, Proceedings of the American Mathematical Society, 60 (1976), pp. 207–210.
[35] J. Kraisler and J. C. Schotland, Collective spontaneous emission and kinetic equations for one-photon light in random media, J. Math. Phys., 63 (2022). 031901.
[36] R.-Y. Lai, K. Ren, and T. Zhou, Inverse transport and diffusion problems in photoacoustic imaging with nonlinear absorption, SIAM J. Appl. Math., 82 (2022), pp. 602–624. arXiv:2107.08118.
[37] J. Lukkarinen and H. Spohn, Kinetic limit for wave propagation in a random medium, Arch. Ration. Mech. Anal., 183 (2007), pp. 93–162.
[38] H. Mahr, Two-photon absorption spectroscopy, in Quantum Electronics: A Treatise, H. Rabin and C. L. Tang, eds., Academic Press, 2012.
[39] L. Margerin, Introduction to radiative transfer of seismic waves, in Seismic Earth: Array Analysis of Broadband Seismograms, A. Levander and G. Nolet, eds., American Geophysical Union, Washington, DC, 2005.
[40] V. A. Markel and J. C. Schotland, Radiative transport for two-photon light, Phys. Rev. A, 90 (2014). 033815.
[41] M. I. Mishchenko, L. D. Travis, and A. A. Lacis, Multiple scattering of light by particles: radiative transfer and coherent back scattering, Cambridge University Press, 2006.
[42] P. P. Mondal, G. Vicidomini, and A. Diaspro, Image reconstruction for multiphoton fluorescence microscopy, Appl. Phys. Lett., 92 (2008). 103902.
[43] A. Pazy, Semigroups of Linear Operators and Applications to Partial Differential Equations, Springer Science & Business Media, 2012.
[44] K. Ren and R. Zhang, Nonlinear quantitative photoacoustic tomography with two-photon absorption, SIAM J. Appl. Math., 78 (2018), pp. 479–503.
[45] K. Ren and Y. Zhong, Unique determination of absorption coeffients in a semilinear transport equation, SIAM J. Math. Anal., 53 (2021), pp. 5158–5184. arXiv:2007.09516.
[46] M.-L. Ren and Z.-Y. Li, Exact iterative solution of second harmonic generation in quasi-phase-matched structures, Optics Express, 18 (2010), pp. 7288–7299.
[47] L. Roux, P. Mareschal, N. Vukadinovic, J.-B. Thibaud, and J.-J. Greffet, Scattering by a slab containing randomly located cylinders: comparison between radiative transfer and electromagnetic simulation, J. Opt. Soc. Am. A, 18 (2001), pp. 374–384.
[48] M. Rumi and J. W. Perry, Two-photon absorption: an overview of measurements and principles, Adv. Opt. Photon., 2 (2010), pp. 451–518.
[49] L. Ryzhik, G. Papanicolaou, and J. B. Keller, Transport equations for elastic and other waves in random media, Wave Motion, 24 (1996), pp. 327–370.
[50] S. E. Skipetrov, Langevin description of speckle dynamics in nonlinear disordered media, Phys. Rev. E, 67 (2003). 016601.
[51] P. Stefanov and Y. Zhong, Inverse boundary problem for the two photon absorption transport equation, SIAM J. Math. Anal., 54 (2022), pp. 2753–2767.
[52] T. Szarvas and Z. Kis, Numerical simulation of nonlinear second harmonic wave generation by the finite difference frequency domain method, J. Opt. Soc. Am. B, 35 (2018), pp. 731–740.
[53] F. Voit, J. Schäfer, and A. Kienle, Light scattering by multiple spheres: comparison between Maxwell theory and radiative-transfer-theory calculations, Opt. Lett., 34 (2009), pp. 2593–2595.
[54] D. Xu, X. Wang, Z. Xu, W. Zhou, and J. Yin, A super-resolution reconstruction algorithm for two-photon fluorescence polarization microscopy, Optics Communications, 499 (2021). 127116.
[55] Y. Yamaoka, M. Nambu, and T. Takamatsu, Fine depth resolution of two-photon absorption-induced photoacoustic microscopy using low-frequency bandpass filtering, Optics Express, 19 (2011), pp. 13365–13377.
[56] C. S. Yelleswarapu and S. R. Kothapalli, Nonlinear photoacoustics for measuring the nonlinear optical absorption coefficient, Optics Express, 18 (2010), pp. 9020–9025.
[57] J. Ying, F. Liu, and R. R. Alfano, Spatial distribution of two-photon-excited fluorescence in scattering media, Appl. Opt., 38 (1999), pp. 224–229.

$\displaystyle\left\|\int_{\mathbb{R}^{d}}\frac{\|\mathbf{p}\|\|\widehat{S}_{0}(z,\mathbf{p})\|d\mathbf{p}}{(2\pi)^{d}}\right\|^{2}$	$\displaystyle\leq\frac{1}{(2\pi)^{2d}}\int_{\mathbb{R}^{d}}\frac{\|\mathbf{p}\|^{2}}{(1+\|\mathbf{p}\|^{2})^{s}}d\mathbf{p}\int_{\mathbb{R}^{d}}(1+\|\mathbf{p}\|^{2})^{\frac{s}{2}}\|\widehat{S}_{0}(z,\mathbf{p})\|^{2}d\mathbf{p}$	(16)
	$\displaystyle=\frac{1}{(2\pi)^{2d}}\left(\int_{\mathbb{R}^{d}}\frac{\|\mathbf{p}\|^{2}}{(1+\|\mathbf{p}\|^{2})^{s}}d\mathbf{p}\right)\\|S_{0}(z,\cdot)\\|_{H^{s}(\mathbb{R}^{d})}^{2}$
	$\displaystyle\leq C_{s}\\|K(z,\cdot)\\|_{H^{s}(\mathbb{R}^{d})}^{2}\left\\|\int_{\mathbb{R}^{d}}W_{0}(z,\cdot,\mathbf{k})d\mathbf{k}\right\\|_{H^{s}(\mathbb{R}^{d})}^{2}$

	$\displaystyle\\|w\\|_{\mathcal{W}_{2}((0,T)\times\Omega\times\mathbb{S}^{d-1})}$	$\displaystyle\leq C\\|(\Sigma_{a}(m_{2})-\Sigma_{a}(m_{1}))\phi_{2}\\|_{L^{2}((0,T)\times\Omega\times\mathbb{S}^{d-1})}$		(33)
		$\displaystyle\leq C\Sigma_{a}^{\prime}(\\|f\\|_{\infty}+\delta)\\|f\\|_{\infty}\\|m_{2}-m_{1}\\|_{L^{2}((0,T)\times\Omega)}$		(33)

	$\displaystyle\\|Fm_{1}-Fm_{2}\\|_{H^{1/2}((0,T)\times\Omega)}$	$\displaystyle=\\|\langle{w}\rangle\\|_{H^{1/2}((0,T)\times\Omega)}\leq C^{\prime}\\|w\\|_{\mathcal{W}_{2}((0,T)\times\Omega\times\mathbb{S}^{d-1})}$		(34)
		$\displaystyle\leq CC^{\prime}\Sigma_{a}^{\prime}(\\|f\\|_{\infty}+\delta)\\|f\\|_{\infty}\\|m_{2}-m_{1}\\|_{L^{2}((0,T)\times\Omega)}.$		(34)

$\displaystyle\|w(z,\mathbf{x},\mathbf{k})\|$	$\displaystyle\leq\int_{0}^{z\wedge\tau_{-}(\mathbf{x},\mathbf{k})}\Sigma_{t}(z-s,\mathbf{x}-s\mathbf{k})e^{-\int_{0}^{s}\Sigma_{t}(z-l,\mathbf{x}-l\mathbf{k})dl}\left\|\frac{R(z-s,\mathbf{x}-s\mathbf{k},\mathbf{k})}{\Sigma_{t}(z-s,\mathbf{x}-s\mathbf{k})}\right\|ds$	(37)
	$\displaystyle\leq\left(1-e^{-\int_{0}^{z\wedge\tau_{-}(\mathbf{x},\mathbf{k})}\Sigma_{t}(z-l,\mathbf{x}-l\mathbf{k})dl}\right)\sup_{0\leq s\leq z\wedge\tau_{-}(\mathbf{x},\mathbf{k})}\left\|\frac{R(z-s,\mathbf{x}-s\mathbf{k},\mathbf{k})}{\Sigma_{t}(z-s,\mathbf{x}-s\mathbf{k})}\right\|$
	$\displaystyle\leq\gamma\sup_{(0,T)\times\Omega\times\mathbb{S}^{d-1}}\left\|\frac{R(z,\mathbf{x},\mathbf{k})}{\Sigma_{t}(z,\mathbf{x})}\right\|$

	$\displaystyle\\|\phi_{\epsilon}(z,\cdot,\cdot)\\|_{\infty}$	$\displaystyle\leq\epsilon\\|h_{3}\\|_{\infty}e^{-\delta z}+\epsilon\int_{0}^{z}e^{-\delta(z-s)}(\\|h_{1}\\|_{\infty}+\delta\\|h_{2}\\|_{\infty})ds$		(47)
		$\displaystyle\leq\epsilon T(\\|h_{1}\\|_{\infty}+\delta\\|h_{2}\\|_{\infty})+\epsilon\\|h_{3}\\|_{\infty}.$		(47)