A Fourier approach to the inverse source problem in an

www.oeaw.ac.at

www.ricam.oeaw.ac.at

A Fourier approach to the inverse source problem in an

absorbing and anisotropic scattering medium

H. Fujiwara, K. Sadiq, A. Tamasan

RICAM-Report 2019-24

A FOURIER APPROACH TO THE INVERSE SOURCE PROBLEM IN AN ABSORBING AND ANISOTROPIC SCATTERING MEDIUM

HIROSHI FUJIWARA, KAMRAN SADIQ, AND ALEXANDRU TAMASAN

ABSTRACT. We revisit the inverse source problem in a two dimensional absorbing and scattering medium and present a non-iterative reconstruction method 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 A. Bukhgeim from the non-scattering to scattering me- dia. 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.

1. INTRODUCTION

This work concerns a Fourier approach to the inverse source problem for radiative transport in a strictly convex domainΩin the Euclidean plane. The attenuation and scattering coefficients are known real valued functions. Generated by an unknown source f, in the steady state case, the density of particlesupz,θqatztraveling in the directionθsolve the stationary transport equation

θ¨∇upz,θq `apzqupz,θq “ ż

S¹

kpz,θ¨θ¹qupz,θ¹qdθ¹`fpzq, pz,θq PΩˆS¹, (1)

whereS¹denotes the unit sphere.

Let Γ_˘ :“ tpz,θq P Γ ˆ S¹ : ˘νpzq ¨ θ ą 0u be the incoming (-), respectively outgoing (+), unit tangent sub-bundles of the boundary; whereνpzqis the outer unit normal atz P Γ. 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 [12, 10, 11, 1, 26], with the most general result for a generic pair of coefficients obtained by Stefanov and Uhlmann [46]. In here we assume that the forward problem is well-posed, and that the outgoing radiationu|_Γ` is measured, and thus the trace u|_ΓˆS¹ is known. Without loss of generalityΩis the unit disc.

In here we show how to recoverffrom knowledge ofuon the torusΓˆS¹and provide an error, and stability estimates.

Whena “k “ 0, this is the classicalX-ray tomography problem of Radon [37], wheref is to be recovered from its integrals along lines, see also [32, 18, 25]. Fora ‰ 0butk “ 0, this is the problem of inversion of the Attenuated Radon transform in two dimensions, solved successfully by Arbuzov, Bukhgeim and Kazantsev [2], and Novikov [34]; see [33, 8, 5] for later approaches.

Date: July 15, 2019.

2010Mathematics Subject Classification. Primary 35J56, 30E20; Secondary 45E05.

Key words and phrases. AttenuatedX-ray transform, Attenuated Radon transform, scattering, A-analytic maps, Hilbert transform, Bukhgeim-Beltrami equation, optical tomography, optical molecular imaging.

1

The inverse source problem in an absorbing and scattering media,a, k ‰ 0, has also been con- sidered (e.g., [24, 44]) in the Euclidean setting, and in [42] in the Riemannian setting. The most general result (k may vary with two independent directions) on the stable determination of the source was obtained by Stefanov and Uhlmann [46]. The reconstruction of the source based on [46] is yet to be realized. When the anisotropic part of scattering is sufficiently small, a convergent iterative method for source reconstruction was proposed in [7]. Based on a perturbation argument to the non-scattering case in [34], the method does not extend to strongly anisotropic scattering.

In addition, it requires solving one forward problem (a computationally extensive effort) at each iteration.

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 such a non-iterative method. Our approach extends the original ideas in [2] from the non-scattering to the scattering media.

Throughout we assume that a P C^2,spΩq, and k and its angular derivative are periodic in the angular variable,k P C_per^1,spr´1,1s;LippΩqq, s ą1{2, and that the forward problem is well posed.

It is known from [46] that for pairs of coefficientspa, kqin an open and dense sets ofC²ˆC², and for anyf PL²pΩq, there is a unique solutionuPL²pΩˆS¹qto the forward boundary value problem.

However, our approach requires a smooth solutionu P H¹pΩˆS¹q. As a direct consequence of [46, Proposition 3.4] the regularity ofuis dictated by its ballistic term. In particular, iff P H¹pΩq, thenuP H¹pΩˆS¹q. With the exception of the numerical examples in Section 7, we assume that the unknown sourcef PH¹pΩq, and thus the unknown solution

(2) uPH¹pΩˆS¹q.

In the numerical experiment we use a discontinuous source, whose successful quantitative recon- struction indicates robustness of the method.

Let upz,θq “ ř8

´8u_npzqe^inθ be the formal Fourier series representation of u in the angular variable θ “ pcosθ,sinθq. Since u is real valued, u_´n “ u_n and the angular dependence is completely determined by the sequence of its nonpositive Fourier modes

ΩQz ÞÑupzq:“ xu₀pzq, u_´1pzq, u_´2pzq, ...y.

(3)

Letknpzq “ _2π¹ şπ

´πkpz,cosθqe^´inθdθ, n P Z, be the Fourier coefficients of the scattering kernel.

Sincekpz,cosθqis both real valued and even inθ,k_npzqare real valued andk_npzq “ k_´npzq.

Throughout this paper the Cauchy-Riemann operatorsB “ pBx `iB<sub>y</sub>q{2 andB “ pBx ´iB<sub>y</sub>q{2 refer to derivatives in the spatial domain. By using the advection operatorθ¨∇“e<sup>´iθ</sup>B `e^iθB, and identifying the Fourier coefficients of the same order, the equation (1) reduces to the system:

Bu₁pzq ` Bu<sub>´1</sub>pzq `apzqu₀pzq “ k₀pzqu₀pzq `fpzq, (4)

and

Bunpzq ` Bun´2pzq `apzqun´1pzq “kn´1pzqun´1pzq, n‰1.

(5)

In particular, the sequence valued map (3) solves the Beltrami-like equation Bupzq `L<sup>2</sup>Bupzq `apzqLupzq “ LKupzq, zP Ω, (6)

where Lupzq “ Lpu₀pzq, u_´1pzq, u_´2pzq, ...q :“ pu_´1pzq, u_´2pzq, ...q denotes the left translation, and

Kupzq:“ pk₀pzqu₀pzq, k_´1pzqu_´1pzq, k_´2pzqu_´2pzq, ...q (7)

is a Fourier multiplier operator determined by the scattering kernel.

Our datau|_ΓˆS¹ yields the trace of a solution of (6) on the boundary, g “u|_Γ“ xg₀, g_´1, g_´2, ...y.

(8)

Bukhgeim’s original theory in [9] concerns solutions of (6) fora“0andK “0. Solutions of Bu`L²Bu“0,

(9)

(called L²-analytic) satisfy a Cauchy-like integral formula, which recovers u in Ω from its trace u|_Γ. In the explicit form in [15], for eachně0,

u_´npζq “ 1 2πi

ż

Γ

u´npzq

z´ζ dz` 1 2πi

ż

Γ

"

dz

z´ζ ´ dz z´ζ

* 8

ÿ

j“1

u_´n´2jpzq

ˆz´ζ z´ζ

˙^j

, ζ PΩ.

(10)

In Section 2 we review the absorbing, non-scattering case. While we follow the treatment in [38], 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 dependence 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, the1{2-gain in (61) below has been known (with different proofs) see [26], and in a more general case than considered here in [46].

The key ingredient in our analysis is an a priori gradient estimate for solutions of the inhomoge- neous Bukhgeim-Beltrami equations, see Theorem 4.2 in Section 4. Our starting point is an energy identity, an idea originated in the work of Mukhometov [30], and an equivalent of Pestov’s identity [42, 49] for the Bukhgeim-Beltrami equation. The proof of the gradient estimate in Theorem 4.2 uses essentially1derivative 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 the method is convergent (see Theorem 6.2), thus recovering the result in [7].

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 be1{2( half way between the ballistic and isotropic regime), and a mean free path of1{5units of length e.g. (meaningful value for fluorescent light scattering in the fog).

2. ABRIEF REVIEW OF THE ABSORBING NON-SCATTERING MEDIUM

In the case fora‰0andK “0, the Beltrami equation (6) can be reduced to (9) via an integrating factor. While this idea originates in [2], in here (as in [38]) we use the special integrating factor proposed by Finch in [15], which enjoys the crucial property of having vanishing negative Fourier modes. This special integrating factor ise^´h, where

hpz,θq:“Dapz,θq ´ 1

2pI´iHqRapz¨θ^K,θ^Kq, (11)

and θ^K is orthogonal to θ, Dapz,θq “ ż8

0

apz `tθqdt is the divergent beam transform of the attenuationa, Raps,θ^Kq “

ż8

´8

a`

sθ^K`tθ˘

dt is the Radon transform of the attenuationa, and the classical Hilbert transformHhpsq “ 1

π ż8

´8

hptq

s´tdtis taken in the first variable and evaluated

at s “ z¨θ^K. The function happeared first in the work of Natterer [32]; see also [8] for elegant arguments that show howhextends analytically (in the angular variable on the unit circleS¹) inside the unit disc. We recall some properties ofhfrom [40, Lemma 4.1], while establishing notations.

Lemma 2.1. [40, Lemma 4.1]AssumeΩĂ R² isC^2,s, s ą 1{2, convex domain. Forp “1,2, let aPC^p,spΩq,są1{2, andhdefined in(11). ThenhPC^p,spΩˆS¹qand the following hold

(i)hsatisfies

θ¨∇hpz,θq “ ´apzq, pz,θq PΩˆS¹. (12)

(ii)hhas vanishing negative Fourier modes yielding the expansions e^´hpz,θq :“

8

ÿ

k“0

α_kpzqe^ikϕ, e^hpz,θq :“

8

ÿ

k“0

β_kpzqe^ikϕ, pz, θq PΩˆS¹, (13)

with (iii)

zÞÑαpzq:“ xα0pzq, α1pzq, α2pzq, α3pzq, ...,y PC^p,spΩ;l1q XCpΩ;l1q, (14)

zÞÑβpzq:“ xβ₀pzq, β₁pzq, β₂pzq, β₃pzq, ...,y PC^p,spΩ;l₁q XCpΩ;l₁q.

(15)

(iv) For anyz PΩ

Bβ₀pzq “0, Bβ₁pzq “ ´apzqβ₀pzq, (16)

Bβ_{k`2</sub>pzq ` Bβ_kpzq `apzqβ<sub>k`1}pzq “0, kě0.

(17)

(v) For anyz PΩ

Bα₀pzq “ 0, Bα₁pzq “ apzqα₀pzq, (18)

Bα<sub>k`2</sub>pzq ` Bα_kpzq ´apzqα_k`1pzq “ 0, k ě0.

(19)

(vi) The Fourier modesα_k, β_k, kě0satisfy α₀β₀ “1,

k

ÿ

m“0

α_mβ_k´m “0, k ě1.

(20)

The Fourier coefficients ofe^˘h define the integrating operatorse^˘Gucomponent-wise for each mď0by

pe^´Guq_m “ pα˚uq_m “

8

ÿ

k“0

α_ku_m´k, and pe^Guq_m “ pβ˚uq_m“

8

ÿ

k“0

β_ku_m´k, (21)

whereα_kandβ_kare the Fourier modes ofe^´hande^hin (13), andα,βas in (14), respectively, (15).

Note thate^˘Gcan also be written in terms of left translation operator as e^´Gu “

8

ÿ

k“0

α_kL^ku, and e^Gu “

8

ÿ

k“0

β_kL^ku, (22)

whereL^k “ L˝ ¨ ¨ ¨ ˝L l jh n

k

is thek-th composition of left translation operator. It is important to note that the operatorse^˘Gcommute with the left translation,re^˘G, Ls “0.

Different from [40], in this work we carry out the analysis in the Sobolev spacesl^2,ppN;H^qpΩqq with the respective norm}¨}_p,qgiven by

l^2,ppN;H^qpΩqq:“

#

u: }u}²_p,q :“

8

ÿ

j“0

p1`j²q^p}u_´j}²_HqpΩq ă 8 +

. (23)

In Proposition 2.1 below we revisit the mapping properties ofe^˘Grelative to these new spaces.

Throughout, in the notation of the norms, the first index p P t0,¹₂,1u refers to the smooth- ness in the angular variable (expressed as decay in the Fourier coefficient), while the second in- dex q P t0,1u shows the smoothness in the spatial variable. The most often occurring is the spacel²pN;L²pΩqq, whenp “ q “ 0. To simplify notation, in this case we drop the double zero subindexes,

}u}² :“ }u}²_0,0 “

8

ÿ

j“0

}u´j}²_L2pΩq.

The traces on the boundary Γ of functions in l^2,ppN;H¹pΩqq are in l^2,ppN;H¹²pΓqq, endowed with the norm

}g}²_p,1 2 :“

8

ÿ

j“0

p1`j²q^p}g_´j}²_H1{2pΓq. (24)

SinceΓis the unit circle, theH^1{2pΓq-norm can be defined in the Fourier domain as follows. For each integerj ě0, if we consider the Fourier expansion of the traceu_´j|_Γ,

u´j|_Γpe^iβq “

8

ÿ

k“´8

u´j,ke^ikβ, for e^iβ PΓ, then

}u_´j}²_H1{2pΓq “

8

ÿ

k“´8

p1`k²q

1

2|u_´j,k|². (25)

In view of (25), ifgP l^2,12pN;H¹²pΓqq, then }g}²¹

2,¹₂ “

8

ÿ

j“0 8

ÿ

n“´8

p1`j²q

1

2p1`n²q

1

2|g_´j,n|². (26)

In the estimates we need the following variant of the Poincar´e inequality obtained by component- wise summation: IfuPl²pN;H¹pΩqq, then

}u}² ďµ

´

}Bu}²` }u|Γ}²_0,¹

2

¯ , (27)

whereµis a constant depending only onΩ; for the unit discµ“2.

Consider the Banach space l₈^1,1pΩq:“

#

α:“ xα0, α1, α2, ...,y:}α}_l^1,1

8 pΩq :“sup

zPΩ 8

ÿ

j“1

j|αjpzq|ă 8 +

. (28)

Proposition 2.1. LetaPC^2,spΩq,są1{2. Thenα,Bα,β,BβP l^1,1₈ pΩq, and for anypP t0,¹₂,1u, qP t0,1u, the operators

e^˘G :l^2,ppN;H^qpΩqq Ñl^2,ppN;H^qpΩqq (29)

are bounded, and satisfy the following estimates

›

›e^´Gu›

› ď }α}_l^1,1

8 pΩq}u}, (30)

›

›e^´Gu›

›1,0 ď2}α}_l^1,1

8 pΩq}u}_1,0, (31)

››e^´Gu›

›0,1 ď

´ }α}_l^1,1

8 pΩq` }Bα}_l^1,1

8 pΩq

¯

}u}_0,1, (32)

››e^´Gu›

›1,1 ď2

´ }α}_l^1,1

8 pΩq` }Bα}_l^1,1

8 pΩq

¯

}u}_1,1. (33)

The same estimate works fore^Guwithαreplaced byβ.

The proof of the Proposition 2.1 can be found in the Appendix.

We remark that casespp “ 0 “ qq, pp “1, q “0q, andpp “0, q “1qin Proposition 2.1 hold fora PC^1,spΩq,s ą1{2, and these are the only properties needed for the lemma below.

Lemma 2.2. LetaPC^1,spΩq,są1{2, ande^˘Gas defined in(21).

(i) If u P l²pN;H¹pΩqqsolvesBu`L<sup>2</sup>Bu`aLu “ 0, thenv “ e^´Gu P l²pN;H¹pΩqqsolves Bv`L²Bv“0.

(ii) Conversely, ifvPl²pN;H¹pΩqqsolvesBv`L<sup>2</sup>Bv“0, thenu“e<sup>G</sup>vPl<sup>2</sup>pN;H<sup>1</sup>pΩqqsolves Bu`L²Bu`aLu “0.

Proof. (i) Let v “ e^´Gu “

8

ÿ

k“0

α_kL^ku. Since u P l²pN;H¹pΩqq, then from Proposition 2.1, vPl²pN;H¹pΩqq. Thenvsolves

Bv`L²Bv“

8

ÿ

k“0

Bα_kL^ku`

8

ÿ

k“0

α_kL^kBu`

8

ÿ

k“0

Bα_kL<sup>k`2</sup>u`

8

ÿ

k“0

α_kL^k`2Bu

“ Bα0u` Bα1Lu`

8

ÿ

k“0

`Bαk`2` Bαk

˘L<sup>k`2</sup>u`

8

ÿ

k“0

αkL^k`

Bu`L²Bu˘

“ Bα0u` Bα1Lu`

8

ÿ

k“0

`Bαk`2` Bαk

˘L<sup>k`2</sup>u`

8

ÿ

k“0

αkL^kp´aLuq

“ Bα₀u``

Bα₁´aα₀˘ Lu`

8

ÿ

k“0

`Bα<sub>k`2</sub> ` Bα<sub>k</sub>´aα<sub>k`1</sub>˘

L^k`2u“0,

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.

3. 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,

kpz,cosθq “

M

ÿ

n“0

k´npzqcospnθq, (34)

for some fixed integerM ě1. Recall that sincekpz,cosθqis both real valued and even inθ,k_npzq are real valued andk_npzq “ k_´npzq,0ď n ďM. We stress here that no smallness assumption on k₀, k_´1, k_´2,¨ ¨ ¨ , k_´M is assumed. Letu^pMqbe the solution of (1) withkas in (34) andu^pMqdenote the sequence valued map

ΩQz ÞÑu^pMqpzq:“ xu^pM₀ ^qpzq, u^pMq_´1 pzq, u^pM_´2^qpzq,¨ ¨ ¨, u^pM_´M^qpzq, u^pMq_´M´1pzq,¨ ¨ ¨ y.

(35)

Let alsoK^pMq denote the corresponding Fourier multiplier operator

K^pMqu^pMqpzq “ pk₀pzqu^pMq₀ pzq, k_´1pzqu^pM_´1^qpzq,¨ ¨ ¨, k_´Mpzqu^pMq_´Mpzq,0,0,¨ ¨ ¨ q.

(36)

The transport equation (1) reduces to the system

Bu^pMq_´1 pzq ` Bu<sup>pM</sup><sub>´1</sub><sup>q</sup>pzq `apzqu^pM₀ ^qpzq “ k₀pzqu^pMq₀ pzq `f^pMqpzq, (37)

Bu^pMq_´n pzq ` Bu<sup>pM</sup><sub>´n´2</sub><sup>q</sup> pzq `apzqu^pMq_´n´1pzq “ k´n´1pzqu^pMq_´n´1pzq, 0ďn ďM´1, (38)

Bu^pMq_´n pzq ` Bu<sup>pM</sup><sub>´n´2</sub><sup>q</sup> pzq `apzqu^pMq_´n´1pzq “ 0, n ěM.

(39)

In sequence valued notation, the system (38) and (39) rewrites:

Bu^pMq`L<sup>2</sup>Bu<sup>pMq</sup>`apzqLu^pMq “LK^pMqu^pMq, (40)

whereK^pMqas in (36).

Sincef P H¹pΩq, the solutionu P H¹pΩˆS¹q, and consequentlyu^pMq P l^2,1pN;H¹pΩqq. We note that in our method we only useu^pMq Pl^2,12pN;H¹pΩqq, indicating that it may apply to rougher sources.

Let the transformation v^pMq “ e^´GL^Mu^pMq, then by Proposition 2.1, v^pMq P l^2,12pN;H¹pΩqq, andv^pMqisL²-analytic:

(41) Bv^pMq`L²Bv^pMq “0.

The trace of the boundaryv^pMq|_Γ is determined by the trace of u^pMq|_Γ“ g “ pg₀, g_´1, g_´2, ...q P l^2,12pN;H^1{2pΓqq, by

v^pMq|_Γ“e^´GL^Mu^pMq|_Γ“e^´GL^Mg.

(42)

By Proposition 2.1,v^pMq|_ΓPl^2,12pN;H^1{2pΓqq.

The Bukhgeim-Cauchy integral formula (10) extends v^pMq from Γ to Ω as L²-analytic map.

From the uniqueness of anL²-analytic map with a given trace, we recovered forn ě0, v^pM_´n^qpzq “ 1

2πi ż

Γ

v_´n^pMqpζq

ζ´z dζ` 1 2πi

ż

Γ

"

dζ

ζ´z ´ dζ ζ´z

* 8

ÿ

j“1

v_´n´2j^pMq pζq

ˆζ´z ζ´z

˙j

, zP Ω.

(43)

Thusv^pMq“ xv^pM₀ ^q, v_´1^pMq, v_´2^pMq, ...yis recovered inl^2,12pN;H¹pΩqq.

We recoverL^Mu^pMq “ xu^pMq_´M, u^pMq_´M´1, u^pM_´M´2^q , ...yinΩby using the convolution formula (21) u^pMq_´n´Mpzq “

8

ÿ

k“0

β_kpzqv^pM_´n´k^q pzq, z PΩ, ně0, (44)

whereβ_k’s as in (13). In particular we recoveredu^pMq_´M´1, u^pMq_´M PH¹pΩq.

By applying 4B to (38), the mode u^pM_´M^q_`1 is then the solution to the Dirichlet problem for the Poisson equation

4u^pM_´M^q_`1 “ ´4B²u^pM_´M´1^q ´4B

”

pa´k_´Mqu^pM_´M^q ı (45a) ,

u^pM_´M^q_{`1</sub>|<sub>Γ</sub> “g<sub>´M`1}, (45b)

where the right hand side of (45a) is known.

Since by constructionu^pMq_´M, u^pMq_´M´1 PH¹pΩq, we have

›

›

›B²u^pMq_´M´1` B

´

pa´k´Mqu^pMq_´M

¯›

›

›

2

H^´1pΩq ď

›

›

›Bu^pMq_´M´1

›

›

›

2

L²pΩq`

›

›

›

´

pa´k´Mqu^pMq_´M

¯›

›

›

2

L²pΩq

ď

›

›

›Bu^pMq_´M´1

›

›

›

2

L²pΩq` }a´k_´M}²_L8pΩq

›

›

›u^pMq_´M

›

›

›

2

L²pΩq

ď

›

›

›u^pMq_´M´1

›

›

›

2

H¹pΩq

` }a´k_´M}²_L8pΩq

›

›

›u^pM_´M^q

›

›

›

2

H¹pΩq

. Sinceg´M`1 PH<sup>1{2</sup>pΓq, the solutionu<sup>pM</sup><sub>´M</sub><sup>q</sup><sub>`1</sub> PH¹pΩqand

›

›

›u^pMq_´M`1

›

›

›

2

H¹pΩq ďC ˆ›

›

›u^pM_´M^q_´1

›

›

›

2

H¹pΩq`

›

›

›u^pM_´M^q

›

›

›

2

H¹pΩq` }g´M`1}²_H1{2pΓq

˙ , (46)

where the constant C depends only on Ωand max

"

1, max

0ďjďM}a´k_´j}²_L8pΩq

*

. Successively all the other modes u^pMq_´M`j for j “ 2,¨ ¨ ¨ , M are computed by solving the corresponding Dirichlet problem for the Poisson equation. To account for the successive accumulation of error we note the following result which can be proven by induction.

Lemma 3.1. Letta_nuandtb_nube sequences of nonnegative numbers, such that a_{n`2</sub> ďcpa<sub>n`1}`a<sub>n</sub>`b_n`2q, ně0,

wherecą0is a constant, then

a<sub>n`2</sub> ď p1`cq^n`2

˜

a₁`a<sub>0</sub>`

n

ÿ

k“0

b_k`2

¸

, ně0.

We applying Lemma 3.1 to (46) and estimate

›

›

›u^pMq_´1

›

›

›

2

H¹pΩq ď p1`Cq^M´1

˜

›

›

›u^pM_´M^q_´1

›

›

›

2

H¹pΩq`

›

›

›u^pM_´M^q

›

›

›

2

H¹pΩq`

M´1

ÿ

j“1

}g´M`j}²_H1{2pΓq

¸ , (47)

and

›

›

›u^pMq₀

›

›

›

2

H¹pΩq ď p1`Cq^M

˜

›

›

›u^pMq_´M´1

›

›

›

2

H¹pΩq`

›

›

›u^pM_´M^q

›

›

›

2

H¹pΩq`

M

ÿ

j“1

}g´M`j}²_H1{2pΓq

¸ . (48)

The sourcef^pMqis computed by f^pMqpzq “ 2Re

´

Bu^pM_´1^qpzq

¯

` papzq ´k₀pzqqu^pMq₀ pzq, (49)

and we estimate

›

›f^pMq›

›

2

L²pΩq ď2

›

›

›u^pMq_´1

›

›

›

2

H¹pΩq`C

›

›

›u^pMq₀

›

›

›

2

H¹pΩq

ď2p1`Cq<sup>M`1</sup>

˜

›

›

›u^pM_´M^q_´1

›

›

›

2

H¹pΩq`

›

›

›u^pMq_´M

›

›

›

2

H¹pΩq`

M

ÿ

j“1

}g_´M`j}²_H1{2pΓq

¸ . (50)

This method is implemented in the numerical experiments in Section 7. Next we analyze the error introduced by truncation.

4. GRADIENT ESTIMATES OF SOLUTIONS TO NONHOMOGENEOUSBUKHGEIM-BELTRAMI EQUATION

When applying the reconstruction method to the data arising from a general scattering kernel kpz,cosθq “

8

ÿ

n“0

k_´npzqcospnθqan 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

Bv`L<sup>2</sup>Bv“Bv`f, (51)

for some specificf and operator coefficientB. Estimates of the gradient for solutions of (51) may be of separate interest, reason for which we treat them here independently of the transport problem.

We start with an energy identity (see [49] forf “0), `a la Mukhometov [30] or Pestov [42].

Theorem 4.1. (Energy identity) Letf P l^2,1pN;L²pΩqqand letB be a bounded operator such that B :l²pN;L²pΩqq Ñl^2,1pN;L²pΩqqandB :l²pN;H¹pΩqq Ñ l^2,1pN;H¹pΩqq.

IfvP l^2,12pN;H¹pΩqqis a solution to(51), then ż

Ω

}Bv}²_l

2dx“ ´2

ż

Ω 8

ÿ

j“1

RexL^2jBv, L^2j´2Bvydx` ż

Ω 8

ÿ

j“0

›

›L^2jBv›

›

2 l2dx (52)

´2 ż

Ω 8

ÿ

j“1

RexL^2jBv, L^2j´2fydx`2 ż

Ω 8

ÿ

j“0

RexL^2jBv, L^2jfydx

` ż

Ω 8

ÿ

j“0

›

›L^2jf›

›

2

l2dx` i 2

ż

Γ 8

ÿ

j“0

xL^2jv,BsL^2jvyds.

Proof. Using the Green’s identity2 ż

Ω

}Bv}²_l

2dx “2 ż

Ω

›

›Bv›

›

2

l2dx`i ż

Γ

xv,Bsvyds,whereBsis the tangential derivative at the boundary, it follows that

ż

Ω

}Bv}²_l

2dx“ ż

Ω

›

›L²Bv´Bv´f›

›

2

l2dx` i 2

ż

Γ

xv,Bsvyds

“ ż

Ω

››L²Bv›

›

2

l2dx´2 ż

Ω

RexL²Bv, Bvydx` ż

Ω

}Bv}²_l

2dx

´2 ż

Ω

RexL²Bv,fydx`2 ż

Ω

RexBv,fydx` ż

Ω

}f}²_l2dx` i 2

ż

Γ

xv,Bsvyds.

For eachn PN, ż

Ω

››L²ⁿBv›

›

2 l2dx“

ż

Ω

››L^2n`2Bv›

›

2

l2dx´2 ż

Ω

RexL<sup>2n`2</sup>Bv, L<sup>2n</sup>Bvydx` ż

Ω

››L²ⁿBv›

›

2 l2dx

´2 ż

Ω

RexL<sup>2n`2</sup>Bv, L<sup>2n</sup>fydx`2 ż

Ω

RexL²ⁿBv, L²ⁿfydx

` ż

Ω

›

›L²ⁿf›

›

2

l2dx` i 2

ż

Γ

xL²ⁿv,BsL²ⁿvyds.

By summing inn, and using lim

nÑ8

ż

Ω

››L²ⁿBv›

›

2

l2dx“0, we conclude the theorem.

We note the general identity [38, Lemma 2.1], for a sequence of nonnegative numbers:

Lemma 4.1. Lettc_nube a sequence of nonnegative numbers. Then piq

8

ÿ

m“0 8

ÿ

n“0

cm`n “

8

ÿ

j“0

p1`jqcj, piiq

8

ÿ

m“0 8

ÿ

n“0

cm`2n“

8

ÿ

j“0

ˆ 1`

Zj 2

^˙

cj, whenever one of the sides in (i) and (ii) is finite.

Proof. (i) By changing the indexj “m`n, formě0, (j´ně0,andn ďj), we get

8

ÿ

m“0 8

ÿ

n“0

c_m`n “

8

ÿ

j“0 j

ÿ

n“0

c_j “

8

ÿ

j“0

c_j

j

ÿ

n“0

1“

8

ÿ

j“0

p1`jqc_j.

(ii) Similarly, by changing the indexj “ m`2n, form ě 0, ˆ

j´2ně0, andn ď Zj

2

^˙

, we get

8

ÿ

m“0 8

ÿ

n“0

c_m`2n“

8

ÿ

j“0

t^j2u ÿ

n“0

c_j “

8

ÿ

j“0

ˆ 1`

Zj 2

^˙

c_j.

Theorem 4.2(Gradient estimate). Letf P l^2,1pN;L²pΩqq, andBbe a smoothing operator such that B :l²pN;L²pΩqq Ñ l^2,1pN;L²pΩqqandB : l²pN;H¹pΩqq Ñl^2,1pN;H¹pΩqqare bounded, and let C_Bą0be such that

}Bv}_1,0 ďCB}v}, @vP l²pN;L²pΩqq, (53a)

}Bv}_1,1 ďC_B}v}_0,1, @vP l²pN;H¹pΩqq.

(53b)

Assume thatBis such that

:“a

2µC_B ă? 2´1, (54)

whereµis the factor in Poincar´e inequality(27).

IfvPl^2,12pN;H¹pΩqqis a solution to the inhomogeneous Bukgheim-Beltrami equation(51), then 0ď }Bv} ď b`?

b²`4ac

2a ,

(55) where

a “1´2´² ą0, b “2}v|Γ}_0,¹

2 `2?

2}f}_1,0, c“²}v|Γ}²_0,¹

2 `π}v|Γ}²¹

2,¹₂ `2}f}²_1,0. Proof. We estimate each term on the right hand side of the energy identity (52). For brevity if we denoteBv“ xb₀, b_´1, b_´2,¨ ¨ ¨ y, then}Bv}²_1,0 “

8

ÿ

j“0

p1`j²q ż

Ω

|b_´j|². I. We estimate the first term in (52):

ˇ ˇ ˇ ˇ ˇ ż

Ω 8

ÿ

j“1

xL^2jBv, L^2j´2Bvydx ˇ ˇ ˇ ˇ ˇ

“ ˇ ˇ ˇ ˇ ˇ ż

Ω 8

ÿ

j“1 8

ÿ

k“0

Bv_´2j´kb_´2j`2´k ˇ ˇ ˇ ˇ ˇ

ď ż

Ω 8

ÿ

j“1

|Bv_´j||jb_´j`1| ď

˜ż

Ω 8

ÿ

j“1

|Bv_´j|²

¸1{2˜ ż

Ω 8

ÿ

j“0

p1`jq²|b_´j|²

¸1{2

ď

?2}Bv}

˜ż

Ω 8

ÿ

j“0

p1`j²q|b´j|²

¸1{2

ď

?2}Bv} }Bv}_1,0

ď

?2C_B}Bv} }v} ď?

2C_B}Bv}? µ

´

}Bv} ` }v|_Γ}_0,¹

2

¯

“a

2µC_B}Bv}²`a

2µC_B }Bv} }v|_Γ}_0,¹

2 ,

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 usep1`xq<sup>2</sup> ď2p1`x²q, in the fifth inequality we use}Bv}_1,0 ď C_B}v}, and in the next to last inequality we use the Poincar´e inequality (27).

II. We estimate the second term in (52):

ż

Ω 8

ÿ

j“1

|xL^2jBv, L^2j´2fy| “ ż

Ω 8

ÿ

j“1 8

ÿ

k“0

|Bv_´2j´kf_´2j`2´k| ď ż

Ω 8

ÿ

j“1

|Bv_´j||jf_´j`1|

ď

˜ż

Ω 8

ÿ

j“1

|Bv´j|²

¸1{2˜ ż

Ω 8

ÿ

j“0

p1`jq²|f´j|²

¸1{2

ď

?2}Bv} }f}_1,0,

where in the first inequality we use Lemma 4.1 part (i), and then Cauchy-Schwarz.

III. We estimate the third term in (52):

ż

Ω 8

ÿ

j“0

›

›|L^2jBv›

›

2 l2 “

ż

Ω 8

ÿ

j“0 8

ÿ

k“0

|b´2j´k|² ď ż

Ω 8

ÿ

j“0

p1`jq|b´j|² ď ż

Ω 8

ÿ

j“0

p1`j<sup>2</sup>q|b´j|<sup>2</sup> ď }Bv}<sup>2</sup><sub>1,0</sub> ďC<sub>B</sub><sup>2</sup> }v}<sup>2</sup> ďµC<sub>B</sub><sup>2</sup> }Bv}<sup>2</sup> `µC_B² }v|Γ}²_0,¹

2 ,

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´e inequality (27).

IV. We estimate the fourth term in (52):

ż

Ω 8

ÿ

j“0

›

›L^2jf›

›

2 l2 “

ż

Ω 8

ÿ

j“0 8

ÿ

k“0

|f_´2j´k|² ď ż

Ω 8

ÿ

j“0

p1`jq|f_´j|² ď ż

Ω 8

ÿ

j“0

p1`j²q|f_´j|² “ }f}²_1,0, where in the first inequality we have used Lemma 4.1 part (i).

V. We estimate the next to the last term in (52):

2 ż

Ω 8

ÿ

j“0

|xL^2jBv, L^2jfy| ď2 ż

Ω 8

ÿ

j“0

›

›L^2jBv›

›l2

›

›L^2jf›

›l2 ď ż

Ω 8

ÿ

j“0

›

›L^2jBv›

›

2 l2`

ż

Ω 8

ÿ

j“0

›

›L^2jf›

›

2 l2

ďµC_B² }Bv}²`µC_B² }v|_Γ}²_0,¹

2 ` }f}²₁,

where in the first inequality we have used Cauchy-Schwarz, in the second inequality we have used the factxyď ¹₂px²`y²q, and in the last inequality we have used estimates from (III) and (IV).

VI. To estimate the last term in (52), we use the fact that Γ is the unit circle and consider the Fourier expansion of the traces of the modesv_´j|_Γ. Fore^iβ PΓ,

v_´jpe^iβq “

8

ÿ

k“´8

v_´j,ke^ikβ, Bβv_´jpe^iβq “

8

ÿ

m“´8

pimqv_´j,me^imβ. Using the parametrization, the last term in (52) becomes

ˇ ˇ ˇ ˇ ˇ ż

Γ 8

ÿ

j“0

xL^2jv,BsL^2jvyds ˇ ˇ ˇ ˇ ˇ

“ ˇ ˇ ˇ ˇ ˇ ż

Γ 8

ÿ

j“0 8

ÿ

k“0

v_´2j´kBsv_´2j´kds ˇ ˇ ˇ ˇ ˇ

“ ˇ ˇ ˇ ˇ ˇ ż

Γ 8

ÿ

j“0

ˆ 1`

Zj 2

^˙

v_´jBsv_´jds ˇ ˇ ˇ ˇ ˇ

“ ˇ ˇ ˇ ˇ ˇ

ż2π 0

8

ÿ

j“0

ˆ 1`

Zj 2

^˙ 8

ÿ

k“´8

v_´j,ke^ikβ

8

ÿ

m“´8

p´imqv_´j,me^´imβdβ ˇ ˇ ˇ ˇ ˇ

“ ˇ ˇ ˇ ˇ ˇ

8

ÿ

j“0

ˆ 1`

Zj 2

^˙ 8

ÿ

k“´8

v_´j,k

8

ÿ

m“´8

p´imqv_´j,m ż2π

0

e^ipk´mqβdβ ˇ ˇ ˇ ˇ ˇ

ď2π

8

ÿ

j“0

p1`j²q¹²

8

ÿ

k“´8

|k| |v_´j,kv_´j,k|

ď2π

8

ÿ

j“0 8

ÿ

k“´8

p1`j<sup>2</sup>q<sup>12</sup>p1`k²q¹² |v´j,k|² “2π}v|Γ}²¹

2,¹₂ , (56)

where in the second equality we have used Lemma 4.1 part (ii), in the first inequality we have used the fact`

1`X_j

2

\˘ ď p1`j²q^1{2, and in the last equality we have used the definition of the norm (26).

Using the above estimates (I)-(VI) for the expressions in (52), we have proved forτ “ ˆż

Ω

}Bv}²_l

2

˙1{2

, thataτ²´bτ ´cď0. Assumption onas in (54) yielda ą0and we have the estimate (55).

For the case whenB “0we obtain the immediate corollary.

Corollary 4.1. Letf Pl^2,1pN;L²pΩqq. IfvPl^2,12pN;H¹pΩqqsolves Bv`L²Bv“f,

(57) then

}Bv}² ď12}f}²_1,0`2π}v|Γ}²¹

2,¹₂ . (58)

Proof. This is the case“0anda“1in (55). We also use

ˆb`?

b²`4c 2

˙2

ďb²`2c.

5. 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_npzq “ _2π¹ şπ

´πkpz,cosθqe^´inθdθ, n P Z. The gain of 1{2 smoothness in the angular variable (see (61) below) has been shown before by different methods in [26], and in a more general case than considered here in [46].

Lemma 5.1. (Smoothing due to scattering) LetM ě1be a positive integer andK be the Fourier multiplier in(7). Assume thatk is such that its Fourier coefficients starting from indexM onward satisfy

γ :“ sup

jěM

p1`jq^pmax }k_´j}₈,}∇_xk_´j}₈(

ă 8, forpą1{2.

(59)

(i) IfvPl²pN;L²pΩqq, then

8

ÿ

n“0

››L^n`MKv›

›

2 ď γ²

pM `1q^2p´1

››L^Mv›

›

2. (60)

(ii)K :l²pN;H¹pΩqq Ñ l^2,12pN;H¹pΩqqis bounded. More precisely,

›

›L^MKv›

›₁

2,1 ď

?2γ pM `1q^p´1{2

›

›L^Mv›

›0,1, @vPl²pN;H¹pΩqq.

(61)

(iii) Moreover, if (59)holds forpě1, thenK :l²pN;L²pΩqq Ñl^2,1pN;L²pΩqqand K :l²pN;H¹pΩqq Ñl^2,1pN;H¹pΩqqare bounded, and

›

›L^MKv›

›1,0 ď γ pM `1q^p´1

›

›L^Mv›

›, @vPl²pN;L²pΩqq, (62)

›

›L^MKv›

›1,1 ď

?2γ pM `1q^p´1

›

›L^Mv›

›0,1, @vPl²pN;H¹pΩqq.

(63)

Proof. (i) LetvP l²pN;L²pΩqq. Then

8

ÿ

n“0

›

›L^n`MKv›

›

2 “

8

ÿ

n“0 8

ÿ

m“0

ż

Ω

|k_´n´m´Mv_´n´m´M|² “

8

ÿ

j“0

p1`jq ż

Ω

|k_´j´Mv_´j´M|²

ďγ²

8

ÿ

j“0

1`j p1`j`Mq^2p

ż

Ω

|v_´j´M|² ďγ²

8

ÿ

j“0

1

p1`j`Mq^2p´1 ż

Ω

|v_´j´M|²

ď γ²

pM `1q^2p´1

8

ÿ

j“0

ż

Ω

|v´j´M|² ď γ² pM `1q^2p´1

›

›L^Mv›

›

2,