GAUGE EQUIVALENCE IN STATIONARY RADIATIVE TRANSPORT THROUGH MEDIA WITH VARYING INDEX OF REFRACTION

Three dimensional anisotropic attenuating and scattering media sharing the same albedo operator have been shown to be related via a gauge transformation. Such transformations define an equivalence relation. We show that the gauge equivalence is also valid in media with non-constant index of refraction, modeled by a Riemannian metric. The two dimensional model is also investigated.

1. Introduction. We consider the inverse boundary value problem of radiation transport in which the attenuation and the scattering kernel are to be determined from the albedo operator. In contrast to the existing works, in this paper the medium has a variable index of refraction and the attenuation is anisotropic.
The radiation is transported in a bounded domain M in R n , n ≥ 2, with smooth boundary and endowed with a Riemannian metric g. If u(x, v) represents the density of particles at position x with velocity vector v in the unit tangent sphere at x, S x M , then the stationary linear transport equation which describes the propagation of particles is The operator D is the derivative along the geodesic flow (see (11) below), which in the case of g being Euclidean is simply Du(x, v) = v · ∇ x u(x, v). The measure 2 STEPHEN MCDOWALL, PLAMEN STEFANOV AND ALEXANDRU TAMASAN dω x (v ) is the volume form on S x M induced from the volume form on T x M determined by g at x; here T x M is the full tangent space to M at x. (See [16] for further details.) The resulting form on SM is the Liouville form and is preserved under the geodesic flow of g. The parameter a(x, v) in (1) is the attenuation coefficient and quantifies the rate at which particles are lost from the point (x, v) in phase space due to absorption and scattering into new directions. The kernel k(x, v , v) is referred to as the scattering coefficient and is proportional to the probability that a particle at position x with velocity v ∈ S x M will scatter to have new velocity v ∈ S x M . The albedo operator, defined more precisely below, maps the flux of particles entering M to the flux of particles exiting M .
In the pioneering paper [5] on inverse transport theory it was pointed out that there can be non-uniqueness to the inverse problem if the attenuation is anisotropic, that is, it depends on direction as well as position. Under the assumption that a = a(x) is isotropic, there is a large collection of uniqueness results to the stationary inverse transport problem under varying assumptions on the parameters. We refer to the recent review paper [2] for an exhaustive account on the subject and mention here some of the results most closely related to the present article. Assuming a Euclidean geometry, and with minimal restrictions guaranteeing only that the forward problem is well-posed, uniqueness of a(x) (dimensions two and above) and k(x, v , v) (dimensions three and above) was proven in [6]; for sufficiently small k this was extended to dimension two in [18], see also [20,21]. The work [8] considers the problem when less information is available, namely one no longer has the angular resolution of the out-going flux, and the spatial part of k is shown to be determined from the corresponding albedo operator. See also [4]. Stability results are proven in [18,8,3,15,22].
The case of a Euclidean metric corresponds to transport in material with a constant index of refraction. If the index of refraction is isotropic, but varying, then (1) can be derived as a limiting case of Maxwell's equations with non-constant (but isotropic) permeability, resulting in a metric which is conformal to the Euclidean metric (see [1] for example). For a general metric, we consider (1) as a model for transport in a medium with varying, anisotropic index of refraction. When the attenuation is assumed isotropic, uniqueness results in Euclidean geometry are extended to the Riemannian setting in [9,10,11,12].
The recent paper [17] deals with anisotropic attenuation a(x, v). In that work it is shown how it is possible to have media of differing attenuation and scattering which yield the same albedo operator, and the non-uniqueness is completely characterized via a gauge transformation as described below. In this paper we extend the results in [17] to media with variable index of refraction. Moreover, we also extend those results to transport in two dimensional domains, a case which was not addressed in [17].
The attenuation a(x, v) is a combination of absorption σ and loss of particles due to scattering. Precisely, a(x, v) = σ(x, v) + SxM k(x, v, v ) dω x (v ). We point out that, even if the absorption coefficient is assumed to be isotropic, σ = σ(x), if k depends on two independent directions, the resulting attenuation a = a(x, v) is anisotropic. In many previous works on the subject, a is assumed to be a function of x alone whereas k is allowed to be fully anisotropic; this requires one to interpret a in a slightly different manner.
Define the "incoming" and "outgoing" bundles where ν x is the unit outer normal vector to the boundary ∂M at x and ·, · is the inner product, each with respect to g at x. The medium is probed with the given radiation and the exiting radiation is detected on Γ + , thus defining the albedo operator A which takes the incoming flux to the outgoing flux at the boundary: Au − = u| Γ+ . In Section 2 we recall the well-posedness results on the forward problem (1) and (2), the mapping properties of the albedo operator and the singular decomposition of its Schwartz kernel. These results appeared in [10,11] in various forms. For the sake of completeness we include in Section 2 some of their proofs in a new and simpler presentation. The parameters (a, k) are assumed admissible, i.e.
Let T be the operator defined by the left hand side of (1). Admissibility (3) guarantees that the second and third terms of T are bounded operators on both L 1 (SM ) and L ∞ (SM ); the first term is unbounded. Since the boundary value problem T u = 0, u| Γ− = u − is equivalent to a non-homogeneous problem with a homogeneous boundary condition, the forward problem is well-posed if T −1 exists as a bounded operator between appropriate spaces. We view T as an unbounded operator on L 1 (SM ) with the domain Given (x, v) ∈ SM , we define the "distance to boundary" functions and set τ = τ + + τ − . The forward problem is well-posed under either one of the two subcritical conditions: or see, e.g., [3,6,7,13,14]. It is proven further in [17,19] that the forward problem is well-posed for a generic set of parameters, namely for an open dense set of (a, k) ∈ L ∞ (SM ) × C M , L ∞ (S x M, L 1 (S x M )) . The proof in [17] is presented in the context of the Euclidean metric, but it remains unchanged in the presence of a general Riemannian metric. For two dimensional domains we consider T defined on Provided that where diam M is the diameter of M with respect to g, such T has a bounded inverse in L ∞ (SM ), see [11]. Here we introduce the somewhat unconventional notation The key observation in [17], which motivates the definition of the gauge equivalence below, is valid in the Riemannian setting as well: Let φ ∈ L ∞ (SM ) be positive with 1/φ ∈ L ∞ (SM ), Dφ ∈ L ∞ (SM ) and such that φ = 1 on ∂SM . Set Then u satisfies (1) if and only ifũ = φu solves Since φ = 1 on Γ, u =ũ there, so the albedo operator A for the parameters (a, k) is indistinguishable from the albedo operatorÃ for the pair (ã,k), i.e. A =Ã.
Gauge equivalent pairs generate the same albedo operator. In this paper we show that the converse also holds. The precise assumptions we make on (M, g) are the following. x (M ) → M is a diffeomorphism. If the dimension of M is two, let κ 0 be the maximum sectional curvature of (M, g); if κ 0 > 0, we assume that diam M < π/ √ κ 0 .
Remark 1. If (M, g) is two dimensional and has constant positive curvature κ 0 , then the restriction on the diameter is equivalent to simplicity of (M, g). If the curvature is not constant, (M, g) may still be simple without satisfying diam M < π/ √ κ 0 , but we must assume this to use the results of [11]. The constraint on the diameter in [11] is imposed to facilitate comparison of geodesic triangles on M to geodesic triangles on constant curvature manifolds where we must guarantee that vertices of these triangles are not conjugate to one another. In fact, in [11] more is assumed, namely diam M < π/(2 √ κ 0 ). The proofs in [11] are written in terms of the stricter diameter assumption, but they can be modified easily to accommodate Assumption 1.2. See Appendix 2 for details.
The two dimensional case differs from the case of dimensions three and higher; we begin by presenting the latter. Throughout, A (respectivelyÃ) are the albedo operators corresponding to parameter pairs (a, k) (respectively (ã,k)). As in [17], one can obtain new uniqueness results for parameter determination if the coefficients have additional symmetry as follows. Suppose that the attenuation depends only on the line determined by ±v: and the scattering coefficient k > 0 obeys: A common situation where (9) holds is in the scattering of light in tissue where it is assumed to depend only on the angle of scattering,  1. If k,k > 0 satisfy (9) then k =k and a =ã + Dw(x) for some function w(x) vanishing on ∂M . In particular, the total attenuations are equal: (9) and further a,ã satisfy (8), then a =ã. Theorem 1.3 and its corollary above are proved in Section 3. In Section 5 we also give a reconstruction formula for the coefficients in the symmetric case.
The two dimensional case is closely related to the work in [11] and it relies on the following new idea. Suppose that A =Ã. Then it follows, just as in higher dimensions, that a andã are related via some gauge φ as in (7). Now define a new scattering kernel k 1 to be the gauge transformation (as in (7)) of k via the gauge φ. We now have (ã, k 1 ) ∼ (a, k) and, as observed above, this implies the corresponding albedo operators are equal. Now we have changed the question into a simpler one: if the albedo operators for the pairs (ã, k 1 ) and (ã,k) are equal, show thatk = k 1 (for then k andk are also related via φ). The advantage is given by the fact that the attenuation is identical. We prove this result in Theorem 4.1, Section 4.
For the two dimensional case, define the class Notice that (a, k) ∈ U Σ,ε automatically implies admissibility and well-posedness of the forward problem. We prove the following results. Corollary 2. Let (M, g) be a two dimensional Riemannian manifold satisfying Assumption 1.2. Given Σ > 0 there exists δ > 0 such that the following holds: Suppose that (a, k), (ã,k) ∈ U Σ,δ are gauge equivalent.
A =Ã then the geodesic ray transform of a−ã, along all geodesics joining boundary points of M , equals zero. Given the simplicity assumption on (M, g), this implies a =ã ( [16]) and the gauge function relating a toã is φ ≡ 1. As in (7), this then yields k =k.
2. The albedo operator and its kernel. Preparatory results. Identifiability of parameters is based on an expansion of the distribution kernel of the albedo operator into terms corresponding to: ballistic particles which do not scatter in M ; particles which undergo a single scattering event in M ; and particles which undergo multiple scattering events. In the current section we present this expansion (Proposition 1, n ≥ 3, and Proposition 2, n = 2), thus indicating the relative strengths of the singularities in each term. We then show how limiting arguments may be applied to extract information from the expansion (Propositions 3 and 4). First we must establish some notation.
We temporarily denote by dσ the volume form on Γ ± , that is the natural restriction of the volume form on SM to Γ ± , see [16] for more details, where this form is denoted by dΣ 2n−2 . In particular, extending dσ as a homogeneous form in |v | of order n − 1, we have that d|v | dσ(x , v ) coincides with the volume form on SM . Then we set The newly defined form dµ has the following invariant property. Fix (x 0 , v 0 ) ∈ Γ ± . Let ∂Ω be another surface so that the geodesic issued from (x 0 , v 0 ) hits it transversally. Then the geodesic flow defines a natural local "projection" near (x 0 , v 0 ), of Γ ± ontoΓ ± , where the latter is related to ∂Ω. Let dμ be the measure onΓ ± defined above. Then the pull back of dμ is dµ.
Given (x, v) ∈ SM , we denote by γ (x,v) (·) the geodesic uniquely determined by We will use the shorthand notation The operator D in (1) is differentiation along the geodesic flow and is defined by If (x i , ξ i ) n i=1 are local coordinates for SM with the (ξ i ) with respect to the natural basis ∂ ∂x i then in these coordinates where Γ i jk are the Christoffel symbols of the Levi-Civita connection of g. If x, y ∈ M denote (momentarily) by v(x, y) ∈ S x M the tangent vector at x of the unit speed geodesic joining x to y. Let d(x, y) be the Riemannian distance between x and y. Define Note thatγ (y,v(y,x)) (d(y, x) − s) = −γ (x,v(x,y)) (s), so when a depends on direction, E(x, y) = E(y, x). It will also be convenient to define, inductively, Proposition 1 below, which describes the terms in the expansion of the kernel of A, is proven in [10], the Euclidean equivalent appearing in [6]. We denote by δ {x,v} (x , v ) the delta-distribution on Γ − with respect to the measure dµ defined by Similarly, δ {x} (y) is the delta distribution on M supported at x. Proposition 1.
[10] Let (M, g) be a smooth Riemannian manifold of dimension n ≥ 2 satisfying Assumption 1.2. Assume that (a, k) are admissible and that the forward problem is well-posed. Then the albedo operator A : Note that k(z(t),ż(t),ẏ(s)) is only defined on the support of the integrand, namely when y(s) = z(r).
As we shall see below, the singular nature of α 0 enables determination of the function E on Γ + × Γ − in all dimensions by way of a limiting argument using an approximate identity. However, it is only in dimensions n ≥ 3 that the same idea can be applied to determine the function k from α 1 . For determination of k in dimension n = 2, we need Proposition 2 which first appeared in [18] for the Euclidean case and is in [11] in the Riemannian case. We first introduce some notation. Given Here, J is a function uniformly bounded 0 < m 1 ≤ J ≤ m 2 < ∞ on Γ + × Γ − (see [11,Proposition 4]).
We begin by showing that α 0 can be extracted from α. This is a re-presentation of [10, Proposition 4.1], presented in a slightly cleaner manner. We include the proof 8
3. Gauge equivalence in dimensions three or higher. In this section we prove Theorem 1.3 and its Corollary 1.
Sincek > 0, φ(y, w ) = φ(y, w) for all w ∈ S y M and φ = φ(y) is independent of w; (15) then implies k =k. By (14) and the above, 4. Gauge equivalence in two dimensions. In this section we prove Theorem 1.4 and its Corollary 2. We begin by proving the intermediate result of unique identifiability of k when a is already determined. This is very closely related to [11,Theorem 1]. Recall the class U Σ,ε , (10). Proof. We sketch the ideas of the proof simply to point out that once we know that the attenuation coefficients are the same (as is assumed in the hypothesis) the proof of [11,Theorem 1], where a is assumed to depend only on x, carries through unchanged when a is allowed to be anisotropic.
Since the attenuation coefficients are the same, α 0 =α 0 . Thus from Proposition 2 we have, for ( (See Proposition 2 and the discussion preceding it for definitions of the terms involved.) Thus where C 1 is a uniform constant depending only on (M, g), y is the point of intersec- and ψ is the angle between w and w. The theorem is proved once we obtain the estimate for then we have The proof of (16) follows from estimates for K 2 φ 0 and K 3 φ 0 , where, with [11]). Upon analysis of the proofs in [11] it is clear that the only manner in which the attenuation coefficient a appears is within the function E, and that the only property of E used is the fact that |E| ≤ 1, uniformly in x ∈ M and in v ∈ S x M . Thus the proof of estimate (16) presented in [11] remains valid and unaltered when a is allowed to depend on direction as well as position.
Proof of Theorem 1.4. The first of (7) follows as in the first part of the proof of Theorem 1.3; let ϕ(x, v) be the gauge function. Now define and let A 1 be the albedo operator corresponding to the pair (ã, k 1 ). We must show that k 1 =k. By definition, (a, k) and (ã, k 1 ) are gauge equivalent. As in Section 1, this implies that A = A 1 ; thusÃ = A 1 , but now the corresponding pairs have the same attenuation: (ã,k) and (ã, k 1 ) respectively. From Theorem 4.1, there exists ε > 0 such that (ã,k), (ã, k 1 ) ∈ U Σ,ε implies that k 1 =k. We now choose δ = δ(ε) to guarantee this. Since we get the uniform bound Choosing δ = e −4Σ diam M we have δ ≤ ε, so (ã,k) ∈ U Σ,δ implies (ã,k) ∈ U Σ,ε , and (a, k) ∈ U Σ,δ implies The proof of Corollary 2 is now identical to that of Corollary 1.

5.
Explicit reconstruction under symmetry assumptions. We demonstrate how one is able to take advantage of the symmetry assumptions (8) and (9) to obtain explicit reconstruction formulae for a and k which are uniquely determined in this case, as shown in Corollary 1. The formulae hold only in dimensions three and higher, the reason for this being that we employ the limiting formula of Proposition 4, the equivalent of which is not known in dimension two. As in Section 2, the albedo operator determines integrals of a along geodesics joining boundary points, but, even under the symmetry condition (8), this does not yet determine the function a(x, v). Next, we use the fact that when n ≥ 3, A also determines This too, on its own, is not sufficient to determine k since the function E is the attenuation along a "broken geodesic" which has not been determined. As we demonstrate below, we can use the symmetry assumptions to first recover k, and then recover a. We must assume higher regularity on the coefficients, namely Recalling Proposition 3, for ε > 0, (x, v) ∈ SM , define Note that J ε is constant along We now extend this to (x, v) ∈ SM .
From the regularity assumption (18), this equality extends pointwise to hold in all of S 2 M . We now set w = −w. We make use of the obvious facts that τ + (y, −w) = τ − (y, w) and γ (y,−w) (s) = (γ (y,w) (−s), −γ (y,w) (−s)). Making the change of variables s ↔ s − τ − (y, w) in the first of the integrals in (24) and s ↔ −s in the second, the left hand side of (24) becomes 0 −τ−(y,w) a(γ (y,w) (s),γ (y,w) (s))+a(γ (y,w) (s), −γ (y,w) (s)) ds a(γ (y,w) (s),γ (y,w) (s)) ds since we are assuming the symmetry (8) for a. This is clearly differentiable in y along the geodesic γ (y,w) and so the right hand side of (24) is also; denote by D w this derivative along the geodesic flow at y in the direction w. Performing this derivative one obtains the reconstruction formula for a(y, w): a(y, w) = 1 4 D w log I 0 (y, w , w) I 0 (y, w, w )J 0 (y, w)J 0 (y, w ) .
Again, ν, ∆t are distances between points on the manifold. If one uses instead sin x ≥ sin( √ κ 0 A)/( √ κ 0 A) x then (27) simply becomes The modification to the proof of Proposition 8 is similar.