A Fourier approach to the inverse source problem in an absorbing and anisotropic scattering medium

We revisit the inverse source problem in a two dimensional absorbing and scattering medium and present a direct reconstruction method, which does not require iterative solvability of the forward problem, using measurements of the radiating flux at the boundary. The attenuation and scattering coefficients are known and the unknown source is isotropic. The approach is based on the Cauchy problem for a Beltrami-like equation for the sequence valued maps, and extends the original ideas of Bukhgeim from the non-scattering to the scattering media. We demonstrate the feasibility of the method in a numerical experiment in which the scattering is modeled by the two dimensional Henyey–Greenstein kernel with parameters meaningful in optical tomography.


Introduction
This work concerns a Fourier approach to the inverse source problem for radiative transport in a strictly convex domain Ω in the Euclidean plane. For simplicity we assume Ω is the unit disc. This assumption simplifies some of the analytical arguments but does not restrict generality, as the domain can be embedded in a larger disc and the coefficients extended by preserving all regularities. The attenuation and scattering coefficients are known real valued functions. Generated by an unknown source f , in the steady state case, the density of particles u(z, θ) at z traveling in the direction θ solve the stationary transport equation where the measure dθ is the usual measure on the unit sphere S 1 normalized such that Let Γ ± := {(z, θ) ∈ Γ × S 1 : ±ν(z) · θ > 0} be the incoming (−), respectively outgoing (+), unit tangent sub-bundles of the boundary; where ν(z) is the outer unit normal at z ∈ Γ. The (forward) boundary value problem for (1) assumes a given incoming flux u on Γ − , In here we assume that there is no incoming radiation from outside the domain, u| Γ− = 0. The boundary value problem is know to be well-posed under various admissibility and subcritical assumptions, e.g. in [1,[9][10][11]22], with the most general result for a generic pair of coefficients obtained in Stefanov and Uhlmann [39]. In here we assume that the forward problem is wellposed, and that the outgoing radiation u| Γ+ is measured, and thus the trace u Γ×S 1 is known.
In here we show how to recover f from knowledge of u on the torus Γ × S 1 . We also provide some error and stability estimates.
When a = k = 0, this is the classical x-ray tomography problem of Radon [32], where f is to be recovered from its integrals along lines, see also [15,21,27]. For a = 0 but k = 0, this is the problem of inversion of the attenuated Radon transform in two dimensions, solved successfully in Arbuzov et al [2], and Novikov in [29,30]; see [4,5,7,17,22,24,28] for later approaches.
The inverse source problem in an absorbing and scattering media, a, k = 0, has also been considered (e.g. [20,37]) in the Euclidean setting, and in [36] in the Riemannian setting. The most general result (k may vary with two independent directions) on the stable determination of the source was obtained in Stefanov and Uhlmann [39]. The reconstruction of the source based on [39] is yet to be realized. When the anisotropic part of scattering is sufficiently small, a convergent iterative method for source reconstruction was proposed in [6]. Based on a perturbation argument to the non-scattering case in [29], the method in [6] does not extend to strongly anisotropic scattering. In addition, it requires solving one forward problem (a computationally extensive effort) at each iteration. A recent result in [38] provides a reconstruction method for the slab geometry, where the angle is limited. The method therein, however, is numerically tested only on a weakly absorbing and scattering medium.
The main motivation of this work is to provide a source reconstruction method that applies to the anisotropic scattering media, with non-small anisotropy. In here we propose a direct method, which does not require iterative solvability of the forward problem. Our approach extends the original ideas in [2] from the non-scattering to the scattering media, and makes essential use of the isotropy of the source. Key to our reconstruction method is the realization that any finite Fourier content in the angular variable of the scattering kernel splits the problem into a non-scattering one and a boundary value problem for a finite elliptic system. The role of the finite Fourier content has been independently recognized in [25].
Throughout we assume that a ∈ C 2,s (Ω), s > 1/2 and k and its angular derivative are periodic in the angular variable, k ∈ C 2 per ([−1, 1]; C 2 (Ω)), and that the forward problem is well posed. It is known from [39] that for scattering kernels k in a vicinity of k in C 2 per ([−1, 1]; C 2 (Ω)), and for any f ∈ L 2 (Ω), there is a unique solution u ∈ L 2 (Ω × S 1 ) to the forward boundary value problem. This will be particularly important in section 6, where the truncation of the Fourier modes yield a complex valued scattering kernel, and thus uniqueness based on the subcritical assumption, e.g. [11], may not be valid. Moreover, our approach requires a smooth solution u ∈ H 1 (Ω × S 1 ). As a direct consequence of [39, proposition 3.4] the regularity of u is dictated by its ballistic term. In particular, if f ∈ H 1 (Ω), then u ∈ H 1 (Ω × S 1 ). With the exception of the numerical examples in section 7, we assume that the unknown source f ∈ H 1 (Ω), and thus the unknown solution u ∈ H 1 (Ω × S 1 ). (2) In the numerical experiment we use a discontinuous source, whose successful quantitative reconstruction indicates robustness of the method. Let u(z, θ) = ∞ −∞ u n (z)e inθ be the formal Fourier series representation of u in the angular variable θ = (cos θ, sin θ). Since u is real valued, u −n = u n and the angular dependence is completely determined by the sequence of its nonpositive Fourier modes Let k n (z) = 1 2π π −π k(z, cos θ)e −inθ dθ, n ∈ Z, be the Fourier coefficients of the scattering kernel. Since k(z, cos θ) is both real valued and even in θ, k n (z) are real valued and k n (z) = k −n (z).
In particular, the sequence valued map (3) solves the Beltrami-like equation where Lu(z) = L(u 0 (z), u −1 (z), u −2 (z), ...) := (u −1 (z), u −2 (z), ...) denotes the left translation, and is a Fourier multiplier operator determined by the scattering kernel. Our data u Γ×S 1 yields the trace of a solution of (6) on the boundary, Bukhgeim's original theory in [8] concerns solutions of (6) for a = 0 and K = 0. Solutions of (called L 2 -analytic) satisfy a Cauchy-like integral formula, which recovers u in Ω from its trace u| Γ . In the explicit form in [12], for each n 0, In section 2 we review the absorbing, non-scattering case. While we follow the treatment in [33], it is in this section that the new analytical framework and notation is introduced. Section 3 describe the reconstruction method for scattering kernels of polynomial depend ence in the angular variable. Except for the numerical section in the end, the remaining of the paper analyzes the error made by the polynomial approximation of the scattering kernel. In section 4 we exhibit the gain in smoothness due to the scattering, in particular, the 1/2-gain in (65) below has been known (with different proofs) see [22], and in a more general case than considered here in [39].
The key ingredient in our analysis is an a priori gradient estimate for solutions of the inhomogeneous Bukhgeim-Beltrami equations, see theorem 4.2 in section 4. Our starting point is an energy identity, an idea originated in [26], and an equivalent of Pestov's identity [31,36,40] for the Bukhgeim-Beltrami equation. The proof of the gradient estimate in theorem 4.2 uses essentially 1 derivative gain in smoothness due to scattering. As a consequence, in theorem 6.1 we establish an error estimate, which yields a stability result for scattering with polynomial angular dependence (see corollary 6.1). Furthermore, in a weakly anisotropic scattering medium and for the exact data, the method is convergent (see theorem 6.2).
The feasibility of the proposed method is implemented in two numerical experiments in section 7. Among the several models for the scattering kernel used in optical tomography [3], we work with the two dimensional version of the Henyey-Greenstein kernel for its simplicity. In this kernel we chose the anisotropy parameter to be 1/2 ( half way between the ballistic and isotropic regime), and a mean free path of 1/5 units of length e.g. (meaningful value for fluorescent light scattering in the fog).

A brief review of the absorbing non-scattering medium
In the case for a = 0 and K = 0, the Beltrami equation (6) can be reduced to (9) via an integrating factor. While this idea originates in [2], in here (as in [33]) we use the special integrating factor proposed in [12], which enjoys the crucial property of having vanishing negative Fourier modes. This special integrating factor is e −h , where  (11). Then h ∈ C p,s (Ω × S 1 ) and the following hold (ii) h has vanishing negative Fourier modes yielding the expansions (iv) For any z ∈ Ω ∂β 0 (z) = 0, (v) For any z ∈ Ω ∂α 0 (z) = 0, (vi) The Fourier modes α k , β k , k 0 satisfy The Fourier coefficients of e ±h define the integrating operators e ±G u component-wise for each m 0 by where α k and β k are the Fourier modes of e −h and e h in (13), and α, β as in (14), respectively, (15). Note that e ±G can also be written in terms of left translation operator as where L k = L • · · · • L k is the kth composition of left translation operator. It is important to note that the operators e ±G commute with the left translation, [e ±G , L] = 0. Different from [35], in this work we carry out the analysis in the Sobolev spaces l 2,p (N; H q (Ω)) with the respective norm · p,q given by In proposition 2.1 below we revisit the mapping properties of e ±G relative to these new spaces. Throughout, in the notation of the norms, the first index p ∈ {0, 1 2 , 1} refers to the smoothness in the angular variable (expressed as decay in the Fourier coefficient), while the second index q ∈ {0, 1} shows the smoothness in the spatial variable. The most often occurring is the space l 2 (N; L 2 (Ω)), when p = q = 0. To simplify notation, in this case we drop the double zero subindexes, The traces on the boundary Γ of functions in l 2,p (N; H 1 (Ω)) are in l 2,p (N; H 1 2 (Γ)), endowed with the norm Since Γ is the unit circle, the H 1/2 (Γ)-norm can be defined in the Fourier domain as follows. For each integer j 0, if we consider the Fourier expansion of the trace u −j Γ , In view of (25), if g ∈ l 2, 1 2 (N; H 1 2 (Γ)), then In the estimates we need the following variant of the Poincaré inequality obtained by component-wise summation: If u ∈ l 2 (N; H 1 (Ω)), then where µ is a constant depending only on Ω; for the unit disc µ = 2. Consider the Banach space Proposition 2.1. Let a ∈ C 2,s (Ω), s > 1/2. Then α, ∂α, β, ∂β ∈ C(Ω; l 1,1 ), and for any p ∈ {0, 1 2 , 1}, q ∈ {0, 1}, the operators e ±G : l 2,p (N; H q (Ω)) → l 2,p (N; H q (Ω)) (29) are bounded, and satisfy the following estimates The same estimates work for e G u with α replaced by β.
where in the last equality we have used (18) and (19). An analogue calculation using the properties in lemma 2.1 (iv) shows the converse. □ We also need estimates for the trace e −G u Γ in terms of u Γ . Let us consider the Fourier expansion of e −h Γ in the spatial variable over the boundary: and let λ(θ) denote the corresponding sequence valued map The Banach space C(S 1 ; l 1,1/2 ) is defined by those λ such that , and the following estimates hold The proof of the proposition 2.2 can be found in the appendix.

Source reconstruction for scattering of polynomial type
This section contains the basic idea of reconstruction in the special case of scattering kernel of polynomial type, for some fixed integer M 1. Recall that since k(z, cos θ) is both real valued and even in θ, k n (z) are real valued and k n (z) = k −n (z), 0 n M . We stress here that no smallness assumption on k 0 , k −1 , k −2 , · · · , k −M is assumed. Let u (M) be the solution of (1) with k as in (38) and u (M) denote the sequence valued map (39) Let also K (M) denote the corresponding Fourier multiplier operator (40) The transport equation (1) reduces to the system In sequence valued notation, the system (42) and (43) rewrites: where K (M) as in (40).
The trace of the boundary v (M) Γ is determined by the trace of u (M) From the uniqueness of an L 2 -analytic map with a given trace, we recovered for n 0, where β k 's as in (13). In particular we recovered u where the right hand side of (49a) is known. Since by construction u .

H Fujiwara et al Inverse Problems 36 (2020) 015005
where the constant C depends only on Ω and max 1, max . Successively all the other modes u We applying lemma 3.1 to (50) and estimate The source f (M) is computed by and we estimate This method is implemented in the numerical experiments in section 7. Next we analyze the error introduced by truncation.

Gradient estimates of solutions to nonhomogeneous Bukhgeim-Beltrami equation
When applying the reconstruction method to the data arising from a general scattering kernel k(z, cos θ) = ∞ n=0 k −n (z) cos(nθ) an error is made due to the truncation in the Fourier series of k. This error is controlled by the gradient of the solution to the Cauchy problem for the inhomogeneous Bukhgeim-Beltrami equation for some specific f and operator coefficient B. Estimates of the gradient for solutions of (55) may be of separate interest, reason for which we treat them here independently of the transport problem.

whenever one of the sides in (i) and (ii) is finite.
Proof.
is a solution to the inhomogeneous Bukgheim-Beltrami equation (55), then Proof. We estimate each term on the right hand side of the energy identity (56). For brevity I. We estimate the first term in (56): where in the first inequality we use lemma 4.1 part (i), in the second inequality we use Cauchy-Schwarz inequality, in the third inequality we use (1 + x) 2 2(1 + x 2 ), in the fifth inequality we use Bv 1,0 C B v , and in the last inequality we use the Poincaré inequality (27). II. We estimate the second term in (56): where in the first inequality we use lemma 4.1 part (i), and then Cauchy-Schwarz. III. We estimate the third term in (56): where in the first inequality we use lemma 4.1 part (i), in the next to the last inequality we use Bv 1,0 C B v , while in the last we use the Poincaré inequality (27). IV. We estimate the fourth term in (56): where in the first inequality we have used lemma 4.1 part (i).
V. We estimate the next to the last term in (56): where in the first inequality we have used Cauchy-Schwarz, in the second inequality we have used the fact xy 1 2 (x 2 + y 2 ), and in the last inequality we have used estimates from (III) and (IV). VI. To estimate the last term in (56), we use the fact that Γ is the unit circle and consider the Fourier expansion of the traces of the Using the parametrization, the last term in (56) becomes where in the second equality we have used lemma 4.1 part (ii), in the first inequality we have used the fact 1 + j 2 (1 + j 2 ) 1/2 , and in the last equality we have used the definition of the norm (26). For the case when B = 0 we obtain the immediate corollary.
then ∂v 2 12 Proof. This is the case = 0 and a = 1 in (59). We also use b

Smoothing due to scattering
In this section we explicit the smoothing properties of the Fourier multiplier operator K in (7) as determined by the appropriate decay of the Fourier coefficients of the scattering kernel k n (z) = 1 2π π −π k(z, cos θ)e −inθ dθ, n ∈ Z. The gain of 1/2 smoothness in the angular variable (see (65) below) has been shown before by different methods in [22], and in a more general case than considered here in [39].

Lemma 5.1 (Smoothing due to scattering).
Let M 0 be an integer and K be the Fourier multiplier in (7). Assume that k is such that its Fourier coefficients starting from index M onward satisfy for p > 1/2.
To prove (ii), let v ∈ l 2 (N; H 1 (Ω)). Then We estimate the first term in (68), where in the first inequality we have used decay of k j 's. Similarly, following the same proof as of the first term (69), the term and the last term

Thus the expression in (68) becomes
To prove (iii), assume that k is such that (63) holds for p 1. Let v ∈ l 2 (N; L 2 (Ω)), then where in the first inequality we have used decay of k j 's.
To prove the last part, let v ∈ l 2 (N; H 1 (Ω)). Then Following the similar proof of (ii), with (1 + j 2 ) instead of (1 + j 2 ) 1/2 , we have Recall that the integrating operators e ±G commute with the left translation, [e ±G , L] = 0. From the equation (40) it is easy to see that the translated sequence solves For a ∈ C 2,s (Ω), s > 1/2, let α, β be as in proposition 2.1, and let λ be as in proposition 2.2. For brevity, let us define Under sufficient regularity on k, its Fourier coefficients satisfy (63) with p = 1.
for any M 0. (75) Proof. The proof follows directly by Bernstein's lemma [16, chapter I, theorem 6.3]. □ As a consequence we obtain the following a priori error estimate for the source reconstructed by the method in section 3. Theorem 6.1 (Error estimate). Let M 0 be an integer, a ∈ C 2,s (Ω), k ∈ C 1,s per ([−1, 1]; Lip(Ω)), for s > 1/2, be known, and f ∈ H 1 (Ω) and u ∈ H 1 (Ω × S 1 ) be unknown functions satisfying (1). Let u be the unknown sequence of nonpositive Fourier modes of u as in (3), and f (M) in (53) be the reconstructed source using the boundary data g Γ ∈ l 2, 1 2 (N; H 1/2 (Γ)). Let δg = u Γ − g be the error in the data. The error δf = f − f (M) in the reconstructed source is estimated by (74), γ M in (75), and C depends only on Ω and max 1, max Proof. Since u ∈ H 1 (Ω × S 1 ), u ∈ l 2,1 (N; H 1 (Ω)), in particular u ∈ l 2, 1 2 (N; H 1 (Ω)). From linearity of the problem, the error δf also satisfy (54): where the second inequality uses (31) in proposition 2.1, and the third inequality uses (66) with p = 1. By Poincare, proposition 2.2, and using (78), From (79) and (77), the expression in (77) becomes   Proof. This is the case γ M = 0 in (76). So c 1 = 0 and result (80) follows. □ We stress that if k is not of polynomial type, then (1 + C) M+1 L M+1 u 2 may not be made arbitrarily small. However, if the anisotropic part of scattering is small enough (same case as in [6]), then we prove the following convergence result for the exact boundary data.. 1]; Lip(Ω)), for s > 1/2, and f ∈ H 1 (Ω). Assume that where µ as in Poincaré inequality (27) and σ as in (74). For each arbitrarily fixed M, let f (M) in (53) be the reconstructed source using the exact boundary data. Then Proof. Since f ∈ H 1 (Ω), the solution u ∈ H 1 (Ω × S 1 ), and consequently u ∈ l 2,1 (N; H 1 (Ω)), in particular u ∈ l 2, 1 2 (N; H 1 (Ω)). The error we make in the source reconstruction is controlled by the sequence valued map q (M) = e −G (u − u M ) ∈ l 2, 1 2 (N; H 1 (Ω)), which solves The operators B M and A M above are operators given by with e ±G as in (21), K as in (7) and K (M) is the truncated Fourier multiplier as in (40).
where in the fourth inequality we use proposition 2.1, in the last inequality we use the Poincaré inequality (27) where the third inequality uses proposition 2.1, and the last inequality uses lemma 5.1 (iii) for p = 1. We can apply theorem 4.
Since u is fixed, the expression L M+1 u → 0, as M → ∞. Therefore, both b M and c M tends to zero and lim

Numerical experiments
In this section we demonstrate the numerical feasibility of the method on two numerical examples. We focus on the numerical results, and leave their realization details for a separate discussion [14]. Although the theoretical results assume a source of square integrable gradient, we shall see below that the numerical reconstruction works even in the case of a discontinuous source.
The numerical experiments consider with the two dimensional Henyey-Greenstein (Poisson) model of scattering kernel and the attenuation coefficient where µ s is the scattering coefficient, and µ a is the absorption coefficient. The parameter g in (86) models a degree of anisotropy, with g = 1 for the ballistic regime and g = 0 for the isotropic scattering regime. In our numerical experiments we work with g = 1/2 and µ s (z) = 5, the latter value shows that, on average, the particle scatters within 1/5 units of length along the path. Let R = (−0.25, 0.5) × (−0.15, 0.15) be the rectangular region, and be the circular region inside the unit disc Ω as shown in figure 1. We consider two examples for the attenuation coefficient a(z) for two different functions µ a (z), while µ s = 5 is the same. In the first example, we work with a C 2 smooth absorption, whereas in the second example we consider attenuation to be discontinuous. In the first example, the attenuation a is C 2 -smooth via some scaled and translated polynomial −(|z| − 1) 3 (6|z| 2 + 3|z| + 1), 0 |z| 1, in the -neighborhoods of ∂B 1 and ∂B 2 : In both examples the data is generated with the same source This is the function we reconstruct by implementing the method in section 3. We generate the boundary measurement by the numerical computation of the forward problem by piecewise constant approximation following the method in [13]. In the forward problem the triangular mesh has 1966 690 triangles, while the velocity direction is split into 360 equi-spaced intervals.
To obtain the data, we disregard the value of the solution inside and only keep the boundary values. The boundary data, u(z, θ) on ∂Ω × S 1 , is represented in figure 2: For each z ∈ ∂Ω (indicated by a cross), we represent the graph { u(z, θ), θ : θ ∈ S 1 } as the (red) closed curve in the polar coordinates z + 2u(z, θ)θ : θ ∈ S 1 . Note that, since u| Γ− = 0 (there is no incoming radiation), the red curve lies outside the domain Ω. The nearly tangential behavior at the boundary is the numerical evidence that the regime is far from ballistic, while the irregular shapes show that it is far from isotropic.
The reconstruction starts with solving (43) with M = 8 via the Bukhgeim-Cauchy integral formula (10), where the infinite series is replaced by a finite sums of up to 64 terms.
In the reconstruction procedure, all the steps, including the integral transforms D

[a], R[a] and H R[a]
, have been processed numerically [14]. In solving the Poisson equations (49a) and (49b), the standard P 1 finite element method (constant and piecewise linear elements) is employed.
The triangulation used in the reconstruction is different from that in the forward problem. In particular the reconstruction mesh consists of 6998 triangles (much less than the 1966 690 triangles used in the forward problem), and is generated without any information of the location of the subsets R, B 1 , and B 2 . In the first numerical experiment the attenuation is given by (87). The reconstructed source f (z) is shown in figure 3 on the left. On the right we show the cross section of f (z) along the dotted diameter y = − √ 3x passing through the origin and the center of B 2 . The reconstructed f (z) shows a quantitative agreement with the exact source in (88). The singularities in the ballistic part of the data carries the singularities of both the source and the attenuation, yielding the same artifacts as in the inversions of the attenuated x-ray transform. In the second example, we apply the proposed algorithm to the case of a discontinuous attenu- otherwise.  While different, the boundary data in the second experiment is graphically indistinguishable from the data in figure 2 in the smooth case. The numerically reconstructed source, shown in figure 4, agrees well with the exact source (88), implying that the proposed algorithm has a potential in reconstruction even if the attenuation admits the discontinuity. In both reconstructions two types of artifacts appear, one type due to the singular support of f and the other due to the singular support of a. Observe how the effect of the latter type of singularities diminishes with an increase in the smoothness of a.  At the level of singular support, the image is similar to that in the non-scattering case (k = 0). This is expected due to the smoothing effect of scattering. A quantitative understanding of the effect of singularities in the attenuation, in the presence of scattering is subject to future work.

Acknowledgment
The authors are grateful to the anonymous referees for their useful comments. The authors thank D Kawagoe for useful discussions of regularity of solution of the forward problem. The work of H Fujiwara was supported by JSPS KAKENHI Grant Nos. 16H02155, 18K18719, and 18K03436. The work of K Sadiq was supported by the Austrian Science Fund (FWF), Project P31053-N32. The work of A Tamasan was supported in part by the NSF Grant DMS-1907097.

Appendix. New mapping properties of the integrating operator e ±G
In this section we prove propositions 2.1 and 2.2.
The discrete convolution below is with respect to the index, while ∂ acts on the spatial variable.
Using Plancherel and the notations above, we estimate the last term in (A.9) as follows. (A.10) Next we estimate each term on the right hand side of (A.10) separately.
I. Estimate the first term in (A.10): where in the first equality we use Plancherel, in the first inequality we use Young's inequality, in the second equality we have used the fact that = sup z∈Ω ∞ k=1 k |∂α k (z)| 2 , and in the last equality we use the norm in (28) and the fact that u −n = u n as u is real valued. II. Estimate the second term in (A.10): where in the first equality we have use Plancherel, in the first inequality we use Young's inequality, in the second equality we have used the fact that sup Ω ∂e −h 2 l1(Z) = sup z∈Ω ∞ k=0 |∂α k (z)| 2 , and in the last equality we use the fact that u −n = u n as u is real valued. III. Estimate the third term in (A.10): where in the first equality we have use Plancherel, in the first inequality we use Young's inequality, in the second equality we have used the fact that Proof of proposition 2.2. For brevity we use the notation n = (1 + |n| 2 ) 1/2 . From convexity of n → n it is easy to see that, for all j, n ∈ Z, n n − j + j , and thus, n 1 2 n − j 1 2 + j 1 2 .
For the proof of the estimate (36) we use ∂ θ (e −h u) = ∂ θ (e −h )u + e −h ∂ θ u, and estimate each of the terms similarly to (A.12). The first term will be bounded by a factor of the product of ∂ θ λ ∞,1, 1 2 and u Γ 0, 1 2 , while the second term will be bounded by a factor of the product of λ ∞,1, 1 2 and u Γ 1, 1 2 . The estimate (37) follows from interpolation of (35) and (36).