Inverse diffusion problems with redundant internal information

This paper concerns the reconstruction of a scalar diffusion coefficient $\sigma(x)$ from redundant functionals of the form $H_i(x)=\sigma^{2\alpha}(x)|\nabla u_i|^2(x)$ where $\alpha\in\Rm$ and $u_i$ is a solution of the elliptic problem $\nabla\cdot \sigma \nabla u_i=0$ for $1\leq i\leq I$. The case $\alpha=\frac12$ is used to model measurements obtained from modulating a domain of interest by ultrasound and finds applications in ultrasound modulated electrical impedance tomography (UMEIT) as well as ultrasound modulated optical tomography (UMOT). The case $\alpha=1$ finds applications in Magnetic Resonance Electrical Impedance Tomography (MREIT). We present two explicit reconstruction procedures of $\sigma$ for appropriate choices of $I$ and of traces of $u_i$ at the boundary of a domain of interest. The first procedure involves the solution of an over-determined system of ordinary differential equations and generalizes to the multi-dimensional case and to (almost) arbitrary values of $\alpha$ the results obtained in two and three dimensions in \cite{CFGK} and \cite{BBMT}, respectively, in the case $\alpha=\frac12$. The second procedure consists of solving a system of linear elliptic equations, which we can prove admits a unique solution in specific situations.

(Communicated by the associate editor name) Abstract. This paper concerns the reconstruction of a scalar diffusion coefficient σ(x) from redundant functionals of the form H i (x) = σ 2α (x)|∇u i | 2 (x) where α ∈ R and u i is a solution of the elliptic problem ∇ · σ∇u i = 0 for 1 ≤ i ≤ I. The case α = 1 2 is used to model measurements obtained from modulating a domain of interest by ultrasound and finds applications in ultrasound modulated electrical impedance tomography (UMEIT), ultrasound modulated optical tomography (UMOT) as well as impedance acoustic computerized tomography (ImpACT). The case α = 1 finds applications in Magnetic Resonance Electrical Impedance Tomography (MREIT).
We present two explicit reconstruction procedures of σ for appropriate choices of I and of traces of u i at the boundary of a domain of interest. The first procedure involves the solution of an over-determined system of ordinary differential equations and generalizes to the multi-dimensional case and to (almost) arbitrary values of α the results obtained in two and three dimensions in [10] and [5], respectively, in the case α = 1 2 . The second procedure consists of solving a system of linear elliptic equations, which we can prove admits a unique solution in specific situations.
1. Introduction. Medical imaging modalities aim to combine high resolution with high contrast between healthy and unhealthy tissues. Optical Tomography and Electrical Impedance Tomography display such high contrasts but often suffer from poor resolution. Ultrasound Tomography and Magnetic Resonance Imaging are high resolution modalities that sometimes suffer from low contrast. The ultrasound modulation of electrical or optical properties of tissues and the combination of simultaneous electrical and magnetic resonance measurements both offer the possibility to combine high resolution with high contrast. For the acquisition of ultrasoundmodulated measurements, we refer the reader to, e.g., [1,5,10,12,14] for works in the mathematical literature. For the acquisition of internal information on electrical conductivities by magnetic resonance imaging, we refer the reader to, e.g., [13,17,18,19].
Mathematically, we aim to reconstruct a scalar diffusion coefficient σ in an elliptic equation from knowledge of internal information of the form H ij (x) = σ 2α (x)∇u i (x) · ∇u j (x) for α ∈ R and 1 ≤ i, j ≤ m, where u i and u j are solutions of the elliptic problem with different boundary conditions; see (1) and (2) below. Coupling impedance (or diffusion) with acoustic waves or magnetic resonance correspond to the cases α = 1 2 and α = 1, respectively. Such information can be obtained from functionals of the form σ 2α (x)|∇u i | 2 (x) by standard polarization (expressions of the form 4ab = (a + b) 2 − (a − b) 2 ).
This problem was first solved in the two dimensional setting in [10] in the case m = 2 and α = 1 2 . The three dimensional setting was addressed in [5] with m = 4 and α = 1 2 . In these papers, the elliptic equation is recast as a system of equations for quantities of the form S i = σ α ∇u i using the elliptic equation and the fact that ∇u i is curl free. This strategy allows one to eliminate σ from the system of equations and solve for the vectors S i . The stable reconstruction of σ is then straightforward. The case α = 1 2 in the setting of non-redundant measurements, i.e., with m = 1 and measurements of the form H = σ|∇u| 2 is considered in [4]. It is shown in that paper that the stable reconstruction of σ may not be possible from such non-redundant measurements. This justifies the analysis of redundant measurements.
The objectives of this paper are twofold. We first generalize the reconstruction of σ to the case of arbitrary space dimension n and almost arbitrary α ∈ R. Assuming that the vectors S i form a frame, we obtain a system of equations for the vectors S i that involves H ij = S i · S j but no longer σ. The resulting system of equations may be seen as an overdetermined nonlinear system of equations. By appropriately choosing the boundary conditions used to construct the internal functionals H ij (x), we obtain a global uniqueness and stability result for the reconstruction of the scalar quantity σ(x). Although several portions of the algorithm generalize to the reconstruction of anisotropic diffusion tensors, we restrict ourselves to the scalar case in this paper. We also describe and investigate the compatibility conditions associated with such a redundant system.
The second objective of the paper is to present a system of elliptic equations for the solutions u i with constitutive parameters that depend on the measurements H ij but not on the unknown diffusion coefficient σ. We show that the system is uniquely solvable when a Fredholm alternative holds. We obtain existence and uniqueness results for the proposed system for all but a discrete number of values of the dimension n and the coefficient α ∈ R.
Both algorithms require boundary conditions for the elliptic solutions that ensure that n of the vectors S i form a frame in R n at each point of the domain of interest. Whereas such a condition is easy to meet in two space dimensions, in dimensions three and higher, the only available technique that guarantees such an independence is based on using complex geometrical optics (CGO) solutions. We generalize here the CGO construction of [5] to the multi-dimensional setting and for almost all values of α.
The inverse diffusion problems with internal functionals considered here are examples of hybrid inverse problem where two imaging modalities are combined to provide both high resolution and high contrast. For recent works on the mathematics of hybrid inverse problems and their many applications in medical imaging, we refer the reader to the articles in the book [20] and to the recent review paper [3].
The rest of the paper is structured as follows. Section 2 presents the main results of the paper on the stable reconstruction of σ from available internal functionals.
The elimination of σ from the system of equations for the vectors S i and the corresponding differential calculus is explained in section 3. The redundant system of equations for the vectors S i is addressed in section 4 while the system of linear equations for the solutions u i is given in section 5. Finally, section 6 presents further reconstruction algorithms in the two dimensional case and analyzes the compatibility conditions satisfied by the redundant data and their potential use.
2. Statement of the main results. Let X be an open convex bounded domain of R n with n ≥ 2. In the following, we address the reconstruction of the scalar conductivity (or diffusion) coefficient σ in the equation where m ≥ n, from knowledge of the interior functionals where α ∈ R is fixed and such that (n − 2)α + 1 = 0. The derivation of the internal functionals (2) in the case α = 1 2 is detailed in [5,14] as examples of synthesized focusing, in [1] in a setup of temporal, physical, focusing, and in [12] by considering thermoelastic effects. The case α = 1 with m = 1 related to MREIT and CDII (Current Density Impedance Imaging) is addressed in [13,17,18,19].
Following a similar approach to [5,10], we first perform the change of unknown functions S i = σ α ∇u i for every i and define We also equip X ⊂ R n with its Euclidean metric g ij = δ ij in the canonical basis (e 1 , . . . , e n ). For a given vector field V = V i e i defined on X, we define the corresponding one-form V ♭ := V i dx i (i.e., by means of the flat operator ). With this notation, we obtain that the vector fields S j satisfy the system of equations where ∧ and d denote the usual exterior product and exterior derivative, respectively. The first equation stems directly from (1) whereas the second one states that the one-form σ −α S ♭ j = du j is exact, therefore closed, and hence d(σ −α S ♭ j ) = 0. When n = 2, 3, equation (5) is recast as: where in dimension n = 2, we define J := 0 −1 1 0 and ∇ ⊥ := J∇, and in dimension n = 3, × denotes the standard cross-product. The available information becomes A crucial hypothesis for our reconstruction procedure is that the m gradients have maximal rank in R n at every point x ∈ X. This hypothesis can be formalized by the somewhat stronger statement: there exists a finite open covering This assumption is equivalent to imposing the following condition on the data

FRANÇ OIS MONARD AND GUILLAUME BAL
where H τ (i) stands for the n×n matrix of elements H While one can always find illuminations such that (6) holds in two dimensions with m = n = 2 and O = {X} (the most preferrable case) by virtue of [2,Theorem 4], higher dimensions can be dealt with using complex geometrical optics solutions provided that σ has enough regularity, as the following lemma shows.
The first step toward an inversion is to express the source term F in terms of a local frame: where H ij denotes the element (i, j) of the matrix H −1 .
Formula (9) was first proved in [5] in the two-and three-dimensional cases with α = 1 2 and is here proved for general n and α ∈ R such that (n − 2)α + 1 = 0. This formula gives us a way to reconstruct F locally from n linearly independent solutions. Assuming condition (7), one can then reconstruct F globally over X.
From lemma 2.2, one can follow two directions to reconstruct the conductivity, which we now describe in more detail in the next two paragraphs. The ODE-based reconstruction procedure. The first approach consists in plugging equation (9) back into the system (4)-(5) and obtain a closed system for the vectors S j . We then show that the resulting system leads to a gradient system, which can then be solved for the vectors S j by ODE integration. Once the vectors S j are reconstructed, one recovers σ from the knowledge of its value at a given point and the fact that ∇ log σ is now known by equation (9). This approach is a generalization of the results of [5] to higher-dimensional settings and general α ∈ R such that (n − 2)α + 1 = 0, and leads to well-posed reconstructions as stated in the following: m ≥ n solutions of (1) generate measurements (H, H ′ ) whose components belong to W 1,∞ (X), and who jointly satisfy condition (7) with the same triple (O, τ, c 0 ). Let also x 0 ∈ Ω i0 ⊂ X and σ(x 0 ), σ ′ (x 0 ) and {S τ (i0)i (x 0 ), S ′ τ (i0)i (x 0 )} 1≤i≤n be given. Let σ and σ ′ be the conductivities corresponding to the measurements H and H ′ , respectively. Then we have the stability estimate: where ε 0 is the error committed at the point x 0 : The solution for the vectors S i and then for log σ requires the solution of full gradient equations of the form ∇u = f (u), where u stands for either unknown. These overdetermined PDEs require compatibility conditions on f if we wish to ensure that their solution does not depend on the path of integration. Theorem 2.3 shows that the reconstruction is unique and stable with respect to the data once a fixed family of integration curves is chosen. The compatibility conditions addressed in section 6 are shown to depend quadratically on the unknown frame. It is therefore difficult to enforce them while solving for the frame S. Nonetheless, depending on the value of α, they may lead to algebraic (i.e. pointwise) reconstructions of all or part of the unknown frame, and may also provide further conditions on the data H ij . Such analyses are carried out in section 6.
Remark 1. Solving a system of equations for the unknown S i may not be efficient numerically. Let S be the matrix whose columns are the n linearly independent vectors S j at a given x. Then S T S = H is known. By the Gram-Schmidt (GS) orthonormalization procedure or by setting R = SH − 1 2 , we can write an equation for an oriented orthonormal frame R; see section 4.3 below. This approach requires that we reconstruct n(n − 1)/2 = dim SO n (R) scalar functions instead of the n × m components of the vector fields {S j }. The only additional constraint is that the transition matrix from S to R satisfies a certain stability property with respect to the data H, see Section 4.3 for details.
simplifies in the sense that the elimination of F is not necessary. Indeed, we show in the next section that knowledge of S i · S j and the constraints dS ♭ j = 0 for 1 ≤ j ≤ m uniquely determine the vectors S j provided they are known at one point x 0 . Once ∇u i is known, the reconstruction of σ may proceed from using (9). Note that, alternatively, the equation (1) may be seen as a transport equation for σ 1−α when α = 1 once the vector field σ α ∇u is known. The stability properties of such a reconstruction are established in [6,7].
The elliptic-based reconstruction procedure. The second approach is novel and consists in injecting equation (9) back into the initial conductivity equations and obtain a strongly coupled elliptic system of the form where the vector fields W ij are known from the data and where the illuminations g i were prescribed in the first place. Here and below, we use the Einstein convention of summation over repeated indices. The vector fields W ij satisfy stability conditions of the form whenever two data sets H and H ′ jointly satisfy condition (7) with the same triple (O, τ, c 0 ). After proving solvability of this system, one is able to reconstruct the functions u i and then to reconstruct σ as described below. Uniqueness and stability of the solution to (11) with respect to the drift fields W ij relies on the fact that is not an eigenvalue of the operator P W : H → H defined by where ∆ −1 D denotes the inverse of the Dirichlet Laplacian on X, and where we have defined the space H : When the coefficients W ij are bounded, we show that the operator P W is compact and its operator norm satisfies the estimate . As a consequence, the system (11) satisfies a Fredholm alternative which will provide uniqueness and stability as stated in the next proposition, for all α ∈ R when n = 2, and for all α but possibly a discrete set (possibly converging to −(n − 2) −1 ) in the Proposition 1 (Stability of the strongly coupled elliptic system). Let vector fields {W ij , W ′ ij } 1≤i,j≤m belong to L ∞ (X) and such that −c −1 F is an eigenvalue of neither P W nor P W ′ . Let u, u ′ be the unique solutions to (11) with same illumination g and respective drift terms W , W ′ . Then we have that u − u ′ ∈ H and satisfies the stability estimate Remark 3. In the case n = 2 or (α = 0 with m = n), we can recast (11) as a coercive system in divergence form, the injectivity of which follows immediately. These cases correspond to c F = 1.
Once the solutions u i are reconstructed, one may reconstruct σ using a formula of the form σ = H 11 /|∇u 1 | 2 . However, such a formula may not offer the best stability estimates. Another reconstruction strategy is deduced from (9), which can be recast locally as where ∇u 1 , . . . , ∇u n denote the n linearly independent gradients. As in the ODEbased reconstruction procedure, we can devise an ODE-based algorithm to reconstruct σ locally from formula (17). We then arrive at the following stability result.
Theorem 2.4. Let the conditions of Proposition 1 be satisfied. Then the corresponding σ, σ ′ satisfy the estimate Note that a necessary and sufficient condition for the unique solvability of (17) on a simply connected domain is that the exterior derivative of the right-hand side (seen as a one-form) vanish by an application of the Poincaré lemma. The compatibility condition that arises here takes the form of a quadratic equation in the components of ∇u i that is difficult to ensure as it depends on the unknowns.
3. Geometric setting and proofs of lemmas 2.2 and 2.1. Defining geometric notation for now, let us first denote the Euclidean orthonormal frame e i = ∂ x i and e i = dx i . For 0 ≤ k ≤ n, Λ k denotes the space of k− forms. We recall the definition of the Hodge star operator ⋆ : We recall the following useful identities, see e.g., [23]: We now prove Lemma 2.2 which is the cornerstone of our explicit reconstructions.
Proof of Lemma 2.2. Because S 1 (x), . . . , S n (x) is a basis of R n at any point x ∈ X, a vector V can be represented in this basis by the following representation (x is implicit here) For j = 1, . . . , n, let us introduce the following 1-forms: where the hat indicates an omission and σ j = (−1) j−1 is the signature of the permutation (1, 2 . . . , n) → (j, 1, . . . , j − 1, j + 1, . . . , n). At each x ∈ Ω, the vector X j (x) obtained from X ♭ j (x) by "raising an index" can also be seen as the unique vector obtained by the Riesz representation lemma that corresponds to the linear form D j : R n → R such that for any V ∈ R n ,

FRANÇ OIS MONARD AND GUILLAUME BAL
We now show that the vector fields X j satisfy a simple divergence equation. We compute The decomposition of X j in the basis S 1 , . . . , S n may be obtained by computing its dotproducts with S 1 , . . . , S n . Indeed, for k = j, there is an l such that i l = k and we have Using formula (20), we deduce that X j admits the expression Plugging this expression into equation (22), and using ∇ · (ϕV ) = ∇ϕ · V + ϕ∇ · V , we obtain We can also recast the previous expression as follows and the proof is complete.
We now give a proof of Lemma 2.1, which guarantees the existence of illuminations that ensure condition (7) and thus justifies the two global reconstruction approaches. The CGO constructions, introduced in [7] in this context, generalize those defined in [5].
Consider the problem ∇ · σ(x)∇u = 0 on R n with σ(x) extended in a continuous manner outside of X and such that σ equals 1 outside of a large ball. The construction requires sufficient smoothness of σ in order to be valid. Let q(x) = − ∆ √ σ σ on R n . We assume that q ∈ H n 2 +1+ε (R n ), which holds if σ − 1 ∈ H n 2 +3+ε (R n ) for some ε > 0, i.e., the original σ |X ∈ H n 2 +3+ε (X). Note that by Sobolev imbedding, σ is of class C 3 (X) while q is of class C 1 (X). With the above hypotheses, we can apply [7, Corollary 3.2] which states the following.
Conclusion. In each of the previous cases, let {g l } 1≤l≤m be the traces of the solutions defined above with m = 2⌊ n+1 2 ⌋. These illuminations generate solutions that satisfy the desired properties of maximal rank and positive determinants. By continuity arguments, any boundary conditionsg l in an open set sufficiently close to g l will ensure that the maximum of the determinants stay bounded from below by c 0 > 0. This concludes the proof of the lemma. 4. The ODE-based method. In this section, we extend the results presented in [5] to general dimension and for a more general class of measurements (described by the coefficient α). We first need to introduce standard geometric notation, without which the derivations become quickly intractable. 4.1. Definitions, notation and identities. We work on a convex set Ω ⊂ R n with the Euclidean metric g(X, Y ) ≡ X · Y = δ ij X i Y j on R n . Following [15], we denote by ∇ the Euclidean connection, i.e. the unique connection that is torsionfree, and compatible with the Euclidean metric in the sense that for smooth vector fields X, Y, Z. On zero-and one-forms, this connection takes the expression: for given vector fields X = X i e i and Y = Y i e i . An important identity for the sequel is the following characterization of the exterior derivative of a one-form ω or equivalently in the Euclidean metric, writing ω = Z ♭ for some vector field Z, where the Lie bracket (commutator) of X and Y coincides with (and thus may be "defined" here as) [X, Y ] = ∇ X Y − ∇ Y X by virtue of the torsion-free property.
A frame refers to an oriented family E = (E 1 , . . . , E n ) of n vector fields over Ω such that for every x ∈ Ω, (E 1 (x), . . . E n (x)) is a basis of T x Ω ≡ R n . For a given frame E, we define the Christoffel symbols (of the second kind) with respect to this frame, by the relations The following very useful identity allows us to compute the Christoffel symbols from inner products and Lie brackets of a given frame (see e.g. [15, Eq. 5.1 p. 69]): where X, Y, Z are smooth vector fields. For a vector X = X j e j , we want to form the matrix of partial derivatives (∂ j X i ) i,j . Geometrically, gradients generalize to tensors via the total covariant derivative, which maps a vector field X to a tensor of type (1, 1) defined by In a given frame E, we may express ∇E i in the basis {E j ⊗ E ♭ k } n j,k=1 of such tensors by writing ∇E i = a ijk E j ⊗ E ♭ k and identifying the coefficients a ijk by writing Equating the two, we obtain the representation The theory of the following sections proves that all partial derivatives of a frame (given in (31)) are uniquely determined by inner products g ij and by Lie brackets, as (29) indicates, or equivalently by exterior derivatives, as (27) expresses. These derivations will be carried out first for the S frame and second for the R frame with values in the space of rotations SO(n, R).

4.2.
The S frame. We now study the properties of the S frame. S is a frame provided that the determinant condition inf x∈Ω det S ≥ c 0 > 0 holds. Our objective in this section is to find an expression for ∇S i that allows us to solve for S i by the method of characteristics. We have seen in the preceding section that this involved calculating the Lie brackets (commutators) of the vectors composing the frame. For 1 ≤ i < j ≤ n, we have Now using (27) we write Plugging this into (32) and using that H kl H kj = δ lj , we obtain the Lie brackets [S i , S j ] for 1 ≤ i < j ≤ n : Returning to the computation of ∇S i using (31), we combine (33) with (29) to arrive at

FRANÇ OIS MONARD AND GUILLAUME BAL
Plugging this expression into (31) (expressed in the S frame), and using H ij H jk = δ ik , we obtain where we have used ∇H jk = −H jp (∇H pq )H qk and have defined Using formulas H jk S j ⊗S ♭ k = I n := e i ⊗e i and H kl (V ·S k )S l = V for any smooth vector field V , we obtain for 1 ≤ i ≤ n Using (23), we observe that ∇S i is equal to a polynomial of degree at most three in the frame S with coefficients involving the known inner products H ij . For each 1 ≤ i, k ≤ n, ∂ k S i is nothing but ∇ e k S i = ∇S i (·, e k ), which can be obtained from (35). Denoting S := [S 1 | . . . |S n ], we are then able to construct the system of equations where Q k β depends only on the data and β is an n 2 -index. This redundant system can then be integrated along any curve (where it becomes a system of ordinary differential equations with Lipschitz right-hand sides ensuring uniqueness of the solution) in order to solve for the matrix-valued function S.

4.3.
The orthonormal R frame. The above system (36) involves a priori n 2 unknowns since the matrix S does not necessarily have any useful symmetries. However, we know the inner products H = S T S, i.e., a matrix of dimension 1 2 n(n + 1). We therefore hope to be able to find a closed-form system involving 1 2 n(n − 1) dimensions. This is the dimension of the orthonormal R frame.
We now provide the details of remark 1. From the frame S, we build an oriented orthonormal frame R = [R 1 | . . . |R n ] (or equivalently, an SO n (R)-valued function) from a matrix-valued function T (x) = {t ij (x)} 1≤i,j≤n that satisfies the relations T T T = H −1 and det T > 0 at every x ∈ Ω, as well as a stability property of the form where C T > 0 depends only on the way we construct T from H. T can either be constructed by the GS procedure or by setting T = H − 1 2 , the positive square root of H −1 . The stability statement (37), first proved in the GS case for n = 2, 3 in [5], can be obtained for both GS and T = H − 1 2 , see [16] for proofs of these statements.
The function R := ST T satisfies everywhere R T R = I n and det R = 1, hence R is an SO n (R)-valued function. The column vectors of S and R transform according to: We also define for 1 ≤ i, k ≤ n We are brief on the derivation of the gradient system for R as it is very similar to that of the S frame. The system of equations (4)-(5) together with the transformation rules (38) allow us to derive the following system of equations for the R frame: From this system, we express F in the R frame as Equation (42) can also be derived directly from (9) and the transformation rules (38). Then, using equation (41) and formula (27), the Lie brackets [R i , , R j ] of the vectors take the form, for 1 ≤ i < j ≤ n: From (43) we deduce the Christoffel symbols relative to the R frame: Finally, in the orthonormal case, the expression of the gradient reduces to ∇R i = Γ j ki R j ⊗ R ♭ k , from which we deduce that (45) As for the S frame, the R.H.S. of (45) depends polynomially on R and on the data. This system can thus be solved for the vectors R i via ODE integration along any curve in a connected domain and provided that we know the R frame at one point. In practice, this system is less expensive to integrate than (36) since the R frame can be locally parameterized with n(n − 1)/2 scalar functions (such as the Euler angles) whereas the S frame requires n 2 scalar functions.

4.4.
Global reconstruction algorithm. The proof of the stability theorem 2.3 can be found in [5] in dimension n = 3 with α = 1 2 (although the proof would be identical in arbitrary dimension). In that paper, the theorem is proved using the system for the rotation matrix R and thus requires the extra stability condition (37). This condition is necessary only if we reconstruct σ via the R frame. The same stability result can be obtained without this requirement if we reconstruct σ via the S frame directly. In the latter setting, the proof is quite similar to the one in [5] with the further simplification that we do not need to change bases when switching subdomain Ω i . The system of ODEs that one must solve based on the gradient system (35) is well-posed since the function S satisfies a priori the uniform bound FRANÇ OIS MONARD AND GUILLAUME BAL and the right-hand side of (36) is Lipschitz in S over the set {S : X → R nm , S ∞ ≤ m H ∞ } as a polynomial of the components of S, and using the fact that the polynomial Q k β are bounded; see [5] for additional details, which we do not reproduce here.
Rewriting the conductivity equation (1) as and plugging (46) into it yields the coupled elliptic system of equations where we have have defined From the last form of W ik , we derive (12) and (13) since the denominators only involve D which is bounded away from zero and the rest is polynomial in the H ij 's and their derivatives. Multiplying (47) by DH pi and writing it in divergence form, one obtain the following equivalent formulation to (47) in variational form: 5.1.2. The case m > n. In the case where we have m > n solutions, we can still define ∇ log σ over the entire domain X using a partition of unity that is subordinate to the open cover O, call it {ϕ i } N i=1 . Then we can define ∇ log σ globally over X by writing where the restrictions are constructed from the n solutions of positive determinant on each Ω i . Each of these restrictions can still be written in the form Thus we can patch these formulas together into a globally defined Plugging this expression into the conductivity equation yields the coupled elliptic system and one arrives at a system of the form (11) by setting W ij := H ik F jk for every 1 ≤ i, j ≤ m. In this case, the stability inequalities (12) and (13) can be derived using the fact that and noticing that on each Ω i , W ij is either zero or locally defined by (48) (and weighed by ϕ i ) whose expression has been proved to be stable. A similar argument holds for proving (13) thanks to the fact that the partition of unity {ϕ i } is the same for two data sets H, H ′ that jointly satisfy (7) with the same triple (O, τ, c 0 ).

5.2.
Uniqueness and stability results. where as well as the stability of its solution with respect to the vector fields W ij . H −1 (X) denotes the dual space of H 1 0 (X). As described in section 2, we apply the inverse of the Dirichlet Laplacian to (51) and obtain the system of integral equations which can be recast in vector notation as Because ∆ −1 D is continuous in the functional setting H −1 (X) → H 1 0 (X) (see [11]), it is clear that f belongs to the space H. We now have the following: Lemma 5.1. Assuming that the vector fields W ij ∈ L ∞ (X), the operator P W : H → H defined in (53) is compact, and its norm satisfies Proof. As can be seen in [11] for instance, the operator ∆ −1 D : L 2 (X) → H 2 (X) is bounded. Therefore, by the Rellich compactness theorem, the operator ∆ −1 D : L 2 (X) → H 1 0 (X) is compact and of norm denoted by ∆ −1 D . Now P is also compact since each of its components is the composition of the continuous operator H ∋ v → W ij · ∇v j ∈ L 2 (X) with the compact operator ∆ −1 D : L 2 (X) → H 1 0 (X). Moreover, for v ∈ H and every 1 ≤ i, j ≤ m, we have the obvious bounds and thus Summing over i proves (54). The proof is complete.
As a consequence of lemma 5.1 and by virtue of standard compact operator theory (e.g. [11,Theorem 6 p 643]), we have the following facts: • 0 is eigenvalue of P W , which corresponds to the case (n − 2)α = −1, a value for α that we exclude from our analysis, • the remaining spectrum of P W is point spectrum and consists of at most a discrete sequence of values that is either finite or converges to zero. Finally, the operator I + c F P W ∈ L(H) satisfies a Fredholm alternative. Therefore it suffices that −c −1 F / ∈ sp (P W ) in order to obtain uniqueness and stability of the solution of (53) and therefore of the solution of (11) as well. The proof of Proposition 1 makes these statements more precise.
Proof of Proposition 1. Let W, W ′ have their coefficients in L ∞ (X) and such that −c −1 F / ∈ sp (P W )∪sp (P W ′ ), and let v, v ′ ∈ H solve the system (51) with respective drift terms W , W ′ and same illumination g. By virtue of the Fredholm alternative, the operators I + c F P W and I + c F P W ′ are invertible with continuous inverses in L(H).
Applying the inverse Dirichlet Laplacian to both systems, we obtain the systems Taking the difference of both systems, the resulting system reads The first difference in the right-hand side of (55) may be bounded by . Applying lemma 5.1 to the operator P W −W ′ , the second difference in the right-hand side of (55) may be bounded by Combining the last two estimates with (55) and the fact that We now conclude with the proof of theorem 2.4.
Proof of theorem 2.4. We focus on the case α = 0 and (n − 2)α = −1. The proof for α = 0 is identical up to small changes in notation. Let H, H ′ have their components in W 1,∞ (X) and jointly satisfy (7) with the same triple (O, τ, c 0 ). Then the families of vector fields W and W ′ have their coefficients in L ∞ (X) and we further assume that −c −1 F / ∈ sp (P W ) ∪ sp (P W ′ ). Let v, v ′ ∈ H solve the system (51) with respective drift terms W and W ′ and same illumination g, and let σ, σ ′ be the corresponding conductivities. Without loss of generality, we work on one of the open sets Ω i ∈ O and renumber the n solutions whose gradients are linearly independent from 1 to n. The result will then hold provided that we have X ⊂ ∪ N i=1 Ω i and thus For a given Ω i ∈ O, and defining V ij := −2α cF D ∇(DH ij ), we write, using equality (17) Similarly to the vector fields W ij (48), the vector fields V ij satisfy estimates of the form Since H is bounded and σ, σ ′ are assumed to be bounded from below by a constant σ 0 > 0, each of the ∇u i is uniformly bounded by H ii /σ 0 ≤ H ∞ /σ 0 . Taking L 2 norms over Ω i and using the triangle inequality in (56), we obtain which by virtue of proposition 1 and estimates (57) yields an estimate of the form Further, from the pointwise relations σ −2α = |∇u 1 | 2 /H 11 and similarly for σ ′ , we write Taking L 2 norms and using the triangle inequality, we obtain that which again yields an estimate of the form Combining (58) and (59), we arrive at for every Ω i ∈ O. This concludes the proof.

5.2.2.
Discussion on the spectrum of P W .
The two-dimensional case. In the case n = 2, [2, Theorem 4] guarantees that we can pick m = n. In so doing and using the form (49) of the elliptic system together with the fact that c F = 1 for all α ∈ R, we arrive at the system ∇ · (DH pi ∇u i ) = 0, u p | ∂X = g p , p = 1, 2.
The weak formulation of the corresponding problem with homogeneous Dirichlet conditions involves the bilinear form Since H −1 is uniformly elliptic over X and inf X D ≥ c 0 , this bilinear form is coercive over H as seen from the following calculation where λ M stands for the largest eigenvalue of H, for which we have, pointwise (x M designates a unit eigenvector associated with λ M ) and hence the estimate Therefore by virtue of the Lax-Milgram theorem, the system (47) admits a unique solution in H. In particular, this shows that −c −1 F is not an eigenvalue of P W in this case for any α ∈ R. The case n ≥ 3. Using the fact that the spectrum of P W is bounded in norm by P W and that P W is compact, the elliptic system admits a unique and stable solution, except for a discrete set of values −c −1 F ∈ [− P W , P W ] possibly converging to zero. In terms of α, this corresponds to almost all values of α ∈ R except a sequence {α k } taking values in the interval [ − PW −1 n−2 , PW −1 n−2 ] and possibly converging to −(n − 2) −1 . The special case α = 0. In this case we have c F = 1. This implies that whenever one can ensure the positivity condition (7) with only m = n solutions (e.g. in even dimension and using lemma 2.1), one can rewrite the system into the form (49) with term 1 − c F = 0, that is Using the same arguments as in the two-dimensional case, this system is coercive and therefore ensures that −1 is not an eigenvalue of P W . Conclusion. As a conclusion of this discussion, the following statements hold: , PW −1 n−2 ] and possibly converges to −(n − 2) −1 . In the case m = n, the value 0 is excluded from the latter interval. 6. Constraints, reconstructions, and compatibility conditions. The ODEbased reconstructions use the full redundancy of the data to construct an overdetermined system of equations for the vectors S i (or R i ) and the vector ∇ log σ. The PDE-based method defines a well-posed system of equations for the scalar quantities u i and an overdetermined system of equations for vector ∇ log σ. Each of these overdetermined systems needs to satisfy compatibility conditions in order to admit a solution. In this section, we aim to extract information from the over-determinacy of the system. We first revisit the two-dimensional case and use the redundancy to extract explicit reconstruction algorithms in the setting α = 1 2 . We then consider the case of arbitrary dimension and show that the compatibility conditions that data must satisfy in order for the aforementioned systems to have solutions take the form of vanishing appropriately defined curvatures together with the cancellation of a given two-form. These conditions generate quadratic functionals of the unknown vectors S i or R i whose pratical applicability is discussed below.
6.1. Reconstructions in two dimensions. In this section, we revisit the twodimensional case which was first solved in [5,10] and generalize the approach to the case α = 1 2 . In that approach, the reconstruction of F = ∇ log σ requires the reconstruction of a function θ : X → S 1 that characterizes the unknown information about the frames S or R. We consider the SO 2 (R)-valued R = (R 1 , R 2 ) frame and parameterize it as R 1 = (cos θ, sin θ) T and R 2 = JR 1 , with J := 0 −1 1 0 . With the notation Φ ij := R i ⊗ R j for i, j = 1, 2, we recast equations (42) and (40) for α ∈ R as follows: We next derive an equation for ∇θ, which by construction is nothing but [R 2 , R 1 ].
We have: where the Christoffel symbols Γ 2 11 and Γ 1 22 are given by By orthonormality the other Christoffel symbols are given by whose expression matches the one given in [5,10] when α = 1 2 . Since F is a function of θ, the above equation is then a non-linear PDE whenever α = 1 2 . This is to be contrasted with the seemingly much nicer case α = 1 2 , whose r.h.s. is independent of θ.
A right-hand side independent of θ can, however, be obtained by taking divergence of both sides of (65) since F = ∇ log σ and ∇ · (J∇) = 0. The equation we obtain is This elliptic PDE requires knowledge of θ at the domain's boundary. Assume that we know u i | ∂X = g i , J i = σ∂ ν u i for i = 1, 2, and σ at the boundary. In this setting, we find that θ| ∂X = arg(t 11 ∇u 1 + t 12 ∇u 2 | ∂X ) = arg((t 11 ∂ t g 1 + t 12 ∂ t g 2 ) t + σ −1 (t 11 J 1 + t 12 J 2 ) ν), with ν and t = Jν the unit outgoing normal vector and its direct orthogonal vector, respectively. Once θ is reconstructed, we know the r.h.s. of (60) and solve for log σ, either by integrating (60) along a curve, or by taking the divergence of both sides of (60) and solving a Poisson equation provided that σ| ∂X is known. Note that the inversion for θ and log σ by means of the elliptic equations (66) and "divergence of (60)" with Dirichlet conditions is unique and Lipschitz-stable in H 2 (X) w.r.t. the data H ij . The details are left to the reader.
We now discuss the compatibility conditions for the gradient equations (60) and (65), which admit a solution only if their respective r.h.s. are curl-free. Such conditions lead to a better understanding of the range of the measurement operator and are necessary to ensure that reconstructions based on ODE integrations do not depend on the choice of integration path. 6.2. Compatibility conditions in two dimensions of space. The compatibility conditions for (65) and (60) are that ∇ · (J∇θ) = 0 and ∇ · (J∇ log σ) = 0, respectively. For α = 1 2 , these equations not only provide constraints on the redundant data, but in fact give us direct information about the unknown coefficients. In the two-dimensional case, they allow us to solve algebraically for cos(2θ), sin(2θ), which in turn characterizes F in terms of the data (and therefore does not require the prior resolution of θ).
Let us first simplify the expression of F as follows: As a result, we are able to express F in a rather compact way