Reconstruction of a fully anisotropic elasticity tensor from 5 knowledge of displacement fields 1 0 2 Guillaume Bal∗ Franc¸ois Monard† Gunther Uhlmann‡ l u J July 6, 2015 2 ] P A Abstract . h We present explicit reconstruction algorithms for fully anisotropic unknown elasticity t a tensors from knowledge of a finite number of internal displacement fields, with applications m to transientelastography. Under certainrank-maximalityassumptions satifiedby the strain [ fields,explicitalgebraicreconstructionformulasareprovided. Adiscussionensuesonhowto fulfill theseassumptions,describingthe rangeofvalidityofthe approach. We alsoshowhow 2 the general method can be applied to more specific cases such as the transversely isotropic v one. 1 6 3 1 Introduction 5 0 . 1 We consider the reconstruction of a fully anisotropic elasticity tensor C = {C } from ijkl 1≤i,j,k,l≤3 0 knowledge of a finite number of displacement fields {u(j)} , solutions of the system of linear j∈J 5 elasticity 1 : v ∇·(C : (∇u+(∇u)T)) = 0 (X), u| = g (prescribed). (1) i ∂X X r Applications for such a theory include the medical imaging modality called elastography. a Elastographyisconcernedwiththereconstructionoftheelasticpropertiesinbiologicaltissues. It has been observed experimentally that certain biological tissues (e.g., muscle fiber [19], or white matter inside the brain [18]) display anisotropic mechanical properties, and the present article aimsatgivingaccesstotheseanisotropicfeatures. Thepresentapproachtoelastographyconsists of two steps. A first step is the reconstruction of the internal elastic displacements, which are ∗Department of Applied Physics and Applied Mathematics, Columbia University, New York NY, 10027; [email protected] †Department of Mathematics, University of Washington; [email protected] ‡Department of Mathematics, University of Washington; Hong Kong University of Science and Technology; [email protected] 1 accurately modeled as solutions of the above system of linear elasticity. Two different methods are typically used to reconstruct such displacements. One such method, called Ultrasound Elastography, consists of probing the elastic displacements by sound waves. We refer the reader to [20, 32,33,34] for moredetails on the method. A recent analysis of such reconstructions from ultrasound measurements is performed in [11]. A second method, called Magnetic Resonance Elastrography,leveragesthedisplacementofprotonsbypropagatingelasticsignalstoreconstruct the elastic waves by magnetic resonance imaging [21, 22, 30]. In this paper, we assume the first step done and consider the quantitative reconstruction of elastic coefficients from knowledge of such elastic displacements. For references on this quanti- tative step of elastography, we refer the reader to [15, 6, 27]. Here, we generalize earlier work performed for scalar second-order equations [13, 12] to the case of systems, and we generalize recent work on the Lam´e system [6, 24] to the fully anisotropic setting. Other hybrid inversions for elliptic systems can be found in the case of Maxell’s system in [8, 14] and in a more general framework in [5]. The approach consists in deriving explicit algebraic formulas reconstructing the unknown parameters locally, based on hypotheses of linear independence (or rank maximality) of func- tionals of the measurements and their partial derivatives. A discussion then follows on what regularity or property (e.g., the Runge approximation property) is required a priori on the un- known parameters so that the hypotheses of reconstructibility may be fulfilled. The idea of constructing local solutions fulfilling certain maximality conditions, which are then controlled from the boundary of the domain via Runge approximation, was also used in the context of reconstruction of conductivity tensors fromknowledge of so-called power density functionals [29] or current density functionals [9]. The above problem of reconstruction of coefficients in partial differential equations from knowledge of internal functionals is sometimes referred to as a hybrid inverse problems; see [2,3,4,5,23,35]forreferenceforhybridinverseproblemsinotherimagingmodalities. Notethat we do not consider here the inverse boundary elasticity problem, for which Lipschitz stability estimates may be obtained only when the Lam´e coefficients are piece-wise constant; see for instance [16] for a recent reference on such a topic not covered here. We now give the main results of the paper in the next section and give an outline of the remainder of the paper there. 2 Main results Preliminary notation and definitions. In what follows, we denote by M (R) the vector 3 space of 3 × 3 real matrices with inner product A : B := tr (ABT) = 3 A B , with i,j=1 ij ij respect to which we recall the orthogonal decomposition M3(R) = S3(R)⊕PA3(R), where the first (second) summand denotes (skew-)symmetric matrices. Let us fix X ⊂ R3 a boundeddomain with smooth boundary for the remainder of the paper. 2 An elasticity tensor C is a fourth-order tensor satisfying the following symmetries C = C , C = C , C = C , 1≤ i,j,k,l ≤ 3, ijkl jikl ijkl ijlk ijkl klij characterized by 21 independent components (instead of 81), where the latter symmetry corre- spondstotheassumptionthat C is hyperelastic. We assumebelow thatC isuniformly pointwise stable over X [26, Ch. 6, Def. 1.9] in the sense that there is a κ> 0 such that 1 ǫ : C(x): ǫ ≥ κ ǫ :ǫ, ∀ x ∈X, ∀ ǫ ∈ S (R). (2) 3 2 For C = {C } an elasticity tensor and ǫ ∈ S (R), we denote C : ǫ := { C ǫ } . ijkl i,j,k,l 3 k,l ijkl kl i,j In some sections below, the tensor C will be represented in the non-tensorPial Voigt notation, an S (R)-valued function c = {c } , where the correspondence of elements c in terms 6 αβ 1≤α,β≤6 αβ of the coefficients C is obtained via the double index mapping 11 7→ 1, 22 7→ 2, 33 7→ 3, ijkl 23,32 7→ 4, 13,31 7→ 5 and 12,21 7→ 6 (for instance, c = C and so on). Hooke’s law, 11 1111 relating stress σ to strain tensors ǫ via the relation σ = C : ǫ, now reads in Voigt notation σ = c ǫ , where we define V V ǫ = (ǫ ,ǫ ,ǫ ,2ǫ ,2ǫ ,2ǫ )T, σ = (σ ,σ ,σ ,σ ,σ ,σ )T. (3) V 11 22 33 23 31 12 V 11 22 33 23 31 12 Most often, we will drop the subscript “V” below as the context will tell us naturally what representation to pick. For (ǫ(1),...,ǫ(6)) in S (R), we define below 3 det(ǫ(1),...,ǫ(6)) = det(ǫ(1),...,ǫ(6)). (4) V R6 V V A crucial fact for further derivations is the following. Lemma 2.1. Under assumption (2), there exists κ′ > 0 such that detc(x) ≥ κ′, ∀ x ∈ X. (5) Main results. Provided linear independence conditions that can be satisfied locally and checked directly on the available measurements, we first provide an explicit reconstruction al- gorithm for fully anisotropic elasticity tensors. Here and below, a typical displacement field is denoted u : X → R3 with corresponding strain tensor ǫ : X → S (R) with components 3 ǫ = 1(∂ u +∂ u ), or, in coordinate-free notation, ǫ = 1(∇u+(∇u)T). ij 2 i j j i 2 We now formulate our main hypotheses in order to set up our reconstruction procedure. These hypotheses, based on algebraic redundancies of various elasticity solutions, force the unknown tensor tolieontheorthogonal ofaspacegenerated byarich enoughsetofdata. Thehypotheses below formulate how some functionals of the data set can be made to generate a hyperplane of S (R), a normal of which can be explicitely constructed and proved to be proportional to C, 6 3 the proportionality factor being reconstructed at the end. As seen below, fulfilling this agenda requires 6+N solutions, where N ≥ d−1 and d denotes the number of scalar components of 3 C (the factor 3 accounts for the fact that each additional solution after the sixth one provides three orthogonality constraints on C). Hypothesis 2.2. Let Ω ⊂ X. A. There exist 6 solutions u(1),...,u(6) of (1) whose strain tensors form a basis of S (R) at 3 every x ∈ Ω. This condition can be summarized as (det is defined in (4)) V inf det(ǫ(1)(x),...,ǫ(6)(x)) ≥ c > 0, for some constant c . (6) 0 0 x∈Ω V B. Assuming A fulfilled, there exists N additional solutions u6+1,...,u6+N giving rise to a family M of 3N matrices whose expressions are explicit in terms of {ǫ(j),∂ ǫ(j), 1 ≤ α ≤ α 3, 1 ≤ j ≤ 6+N} (see (12) for detail), and such that they span a hyperplane of S (R) at 6 every x ∈ Ω. This condition can be summarized as inf N(M′) :N(M′) ≥ c > 0, for some constant c , (7) 1 1 x∈Ω M′⊂MX, #M′=20 where N is an operator generalizing the cross-product, defined in (13). Remark 2.3. It should be noted that these hypotheses are stable under smooth perturbation of the boundary conditions generating the displacement fields u(1),...,u(6+N). This is because Hypotheses 2.2.A-B are expressed in terms of functionals which depend polynomially on the components of displacement fields and their derivatives up to second order, which are in turn continuous functionals of their boundary conditions (see, e.g., Sec. 5.2.1 for appropriate topolo- gies). Remark 2.4. In Hypothesis 2.2.B, the number N depends on the type of isotropy of the tensor C. In the most general, 21-parameter case, spanning a hyperplane of S (R) with 3N additional 6 constraints suggests N ≥ 7, and the total number of displacement fields needed is then 6+7= 13. On an open subset Ω ⊂ X where Hypotheses 2.2.A-B hold, we then derive an explicit reconstruction algorithm, reconstructing C over Ω. This is done as follows: in Voigt notation, decomposethe S (R)-valued function C as the productof a scalar function τ times a normalized 6 anisotropic structure C such that C = τC (C is such that its Voigt counterpart ˜c has determi- 1 nant 1, so τ = (detc˜)e6). Algebraic maneipuelations allow us to obtain pointwise orthogonality constraints which, under hypotheses (6)-(7), are numerous enough to reconstruct C pointwise. Once C is reconstructed, an equation for ∇logτ is derived over Ω, leading to the reconstruction e of τ by solving either a transport or a Poisson equation. Finally, once C = τC is known, ad- e ditional stability is recovered on certain components of C by deriving equations reconstructing e 4 the third-order tensor div C (whose components may be written as ∂ C , 1 ≤ j,k,l ≤ 3). i ijkl Such components, although linear combination of first-orderderivatives of the components of C, are reconstructed with the same stability as C itself. Such stability improvements have already been observed in e.g. [10, 9]. The approach just described yields unique and stable reconstructions over Ω in the sense of the following theorem. Theorem 2.5. Suppose that over some open set Ω ⊂ X, hypotheses 2.2.A-B hold for two families of displacement fields {u(j)}6+N and {u′(j)}6+N corresponding to elasticity tensors C j=1 j=1 and C′. Then C and C′ can each be uniquely reconstructed over Ω from knowledge of their corresponding solutions, with the following stability estimate for every integer p ≥ 0 N+6 kC −C′k +kdiv C −div C′k ≤ K kǫ(j)−ǫ′(j)k . (8) Wp,∞(Ω) Wp,∞(Ω) Wp+1,∞(Ω) Xj=1 Remark 2.6. Note that if the elasticity tensor is split into the product of an unknown scalar function τ times a known anisotropic tensor, the stability of the problem of reconstructing τ from fields u(j) is better-posed, i.e. involves the loss of one derivative instead of two, according to the statement N+6 kτ −τ′kWp+1,∞(Ω) ≤ K kǫ(j)−ǫ′(j)kWp+1,∞(Ω). (9) Xj=1 In the context of this reconstruction approach, an elasticity tensor is then reconstructible on some given open set Ω (or globally on X) if there exist displacement field solutions of (1) fulfillingHypotheses 2.2.A-B throughoutΩ. In this context, classifying reconstructible elasticity tensorsthenconsistsinfindingoutinwhatsituationsonecanconstructsuchdisplacementfields. Two such situations are described in the following theorem: (i) a fully anisotropic tensor that is close to constant is globally reconstructible from well-chosen displacement fields whose traces are explicitly given, and (ii) an elasticity tensor with smooth enough components which satisfies the Runge approximation property (see Section 5.2.3, in particular, Eq. (39)) is, in principle, locally reconstructible from knowledge of its displacement fields. Theorem 2.7 (Reconstructibility of elasticity tensors). In either of the following cases, there exists a non-empty open set of smooth enough boundary conditions generating displacement fields characterizing C uniquely and stably in the sense of Theorem 2.5. (i) C is C3-close to constant. (ii) C is smooth (at least of class C3) and satisfies the Runge approximation property. 5 An application of particular interest is the reconstruction of transversely isotropic (TI) elas- ticity tensors. Such tensors have, as of yet, the highest type of anisotropy for which the Runge approximation has been proved (see [31]), so that they are covered in case (ii) of Theorem 2.7. We also explain in Section 5.1 how to generate explicit boundary conditions reconstructing a constant TItensor which, by case (i) of Theorem 2.7, can beusedto reconstruct a near-constant TI tensor. Outline. The rest of the paper is organized as follows. Section 3 covers the proof of Lemma 2.1. Section 4 covers the reconstruction procedure as well as its stability property (proof of Theorem 2.5). Section 5 covers the two ways described in Theorem 2.7 to fulfill Hypotheses 2.2.A-B for certain classes of tensors, thereby establishing their unique and stable (in the sense of Theorem 2.5) reconstructibility from displacement fields. 3 Proof of Lemma 2.1 Suppose the elasticity tensor C satisfies (2). In the study of elastic eigentensors in [28], it is established that C can have at most six distinct eigencouples (ǫ,λ) such that the relation C : ǫ = λǫ is satisfied. The symmetries of C and hypothesis (2) imply that the λ’s are real and that all of them satisfy λ ≥ κ. It is also mentioned in [28] that the eigenvalues λ are precisely the eigenvalues of the S (R)-valued representation of C in the form 6 c′ = {2χ(i)+2χ(j)cij}1≤i,j≤6, where χ(i) = (cid:26) 10 iiff ii == 41,,52,,63,, and where c denotes the Voigt representation of C presented in the introduction. It is then ij immediate that detc′ = 8detc. Moreover, as mentioned above, all eigenvalues of c′ match the eigenvalues of C and therefore satisfy the estimate λ ≥ κ. Therefore 1 κ6 detc= detc′ ≥ , 8 8 hence (5) holds with κ′ = κ6. 8 4 Reconstruction algorithm and its stability As seen in Theorem 2.5, the left hand side of (8) contains two terms whose product forms the elasticity tensor C, reconstructed with stability estimates in different norms: we decompose C into the product τC, where C contains the rescaled anisotropic structure of C, defined by a normalizing condition (specifically, if c˜ describes C in Voigt notation, then detc˜ = 1) and τ is e e the remaining scalar factor. That this is possible comes from estimate (5), which states that e detc is uniformly bounded away from zero throughout X. The next two sections focus on the successive reconstruction of C, then τ upon assuming that Hypotheses 2.2.A-B holds. e 6 4.1 Reconstruction of the anisotropy In this section, we will use the Voigt notation, so that strain tensors will be represented as R6- valued funtions ǫ ≡ ǫ . In this notation, theelasticity tensor becomes an S (R)-valued function, V 6 characterized by 21 scalar functions c = {c } satisfying c = c and where Hooke’s αβ 1≤α,β≤6 αβ βα law is expressed as a regular matrix-vector product σ = c ǫ , with (σ ,ǫ ) as in (3). This V V V V makes the problem tractable via algebraic manipulations on matrices instead of 4-tensors. Using the Voigt notation, the elasticity system (1) takes the form ∂ 0 0 0 ∂ ∂ 1 3 2 DV ·(c ǫ)= 0, DV :=  0 ∂2 0 ∂3 0 ∂1 , (10) 0 0 ∂ ∂ ∂ 0 3 2 1   and where, for a M (R)-valued function A and a scalar function f, we have the identity 6 D ·(fA)= (D f)·A+fD ·A. (11) V V V As the reconstruction approach is local (even pointwise for the anisotropic part), we assume to have 6 elasticity solutions {u(j)} whosecorrespondingstrain tensors {ǫ(j)(x)} form 1≤j≤6 V 1≤j≤6 a basis of R6 for every x of some subdomain Ω⊂ X. Any additional solution u(p) with p ≥ 7 is such that for x ∈ Ω, ǫ(p)(x) decomposes uniquely into the basis above as j 6 det (ǫ(1)(x),...,ǫ(p)(x),...,ǫ(6)(x)) ǫ(p)(x)= µpj(x)ǫ(j)(x), where µpj := V z }| { , 1 ≤ j ≤ 6. det (ǫ(1)(x),...,ǫ(6)(x)) Xj=1 V Plugging this equality into the elasticity equation, we obtain that 6 6 0= D ·(cǫ(p)) = D ·(µ cǫ(j)(x )) = (D µ )·cǫ(j)+µ D ·(cǫ(j)), V V pj 0 V pj pj V Xj=1 Xj=1 which, since D ·(cǫ(j)) = 0, implies that V 6 (D µ )·cǫ(j) = 0. V pj Xj=1 This last equation can be seen as three scalar orthogonality constraints on the tensor c in the inner productstructureof S (R) which we denote by A:B := tr (ABT)= 6 A B , where 6 i,j=1 ij ij the matrices that are othogonal to C are directly known from available meaPsurements. The last 7 equation is equivalent to c: M(p),1 = c: M(p),2 = c: M(p),3 = 0, where 6 M(p),1 = (∂ µ 0 0 0 ∂ µ ∂ µ )⊗ǫ(j), 1 pj 3 pj 2 pj Xj=1 6 (12) M(p),2 = (0 ∂ µ 0 ∂ µ 0 ∂ µ )⊗ǫ(j), 2 pj 3 pj 1 pj Xj=1 6 M(p),3 = (0 0 ∂ µ ∂ µ ∂ µ 0)⊗ǫ(j). 3 pj 2 pj 1 pj Xj=1 Note that since c is orthogonal to A (R), one could replace the matrices M(p),i with their 6 symmetrized versions. Notice that the components of these matrices are first partial derivatives of rational functions of strain tensors. If enough linear constraints of the form (12) are available from a rich enough set of measure- ments, that isto say, ifenough suchmatrices areavailable andformahyperplaneofS (R)at x , 6 0 then the tensor c(x ), constrained to be perpendicular to this hyperplane, will be determined 0 up to a multiplicative constant. The reconstruction can be done via a generalization of the cross-product, as used for instance in [29] for similar purposesin the context of the conductivity equation. Define {m }21 a basis of S (R) and given M = {M }20 ⊂ S (R), define the oper- j j=1 6 j j=1 6 ator N : S (R)20 → S (R) by expanding the formal determinant below with respect to its last 6 6 row M :m ··· M :m 1 1 1 21 N(M) := 1 (cid:12)(cid:12) ... ... ... (cid:12)(cid:12). (13) (cid:12) (cid:12) det(m1,··· ,m21)(cid:12) M20 : m1 ··· M20 :m21 (cid:12) (cid:12) (cid:12) (cid:12) m ··· m (cid:12) 1 21 (cid:12) (cid:12) (cid:12) (cid:12) N is a 20-linear, alternating map that does not(cid:12)depend on the choice of basi(cid:12)s for S (R). N(M) 6 is a vector thatis normalto thehyperplanespannedby M whenM is linearly independent, zero otherwise. In particular, if M is a family of matrices known to be orthogonal to a given matrix m′, then N(M) is either zero (if dimspan M < 20) or proportional to m′ (if dimspan M = 20). Inlightofthislastcomment, andassumingthatarichenoughsetofsolutionsoftheelasticity system (1) gives rise to a family of matrices M of the form (12), with cardinality greater than 20 and spanning a hyperplane of S (R) at a given point x ∈ X, then for any given 20-uple 6 0 M′ ⊂ M, N(M′) is either zero or proportional to c(x0). Normalizing c˜(x0) = det(c(x0))−61c(x0) and enforcing the condition ˜c > 0 (this is because in the proof of Lemma 2.1, we notice that 11 8 the top-left 3×3 block of c matches that of c′ which is a symmetric positive definite matrix, so one must have c > 0, c > 0 and c > 0), then we have the equality 11 22 33 (±)M′N(M′) = (det(N(M′)))16c˜(x0), forevery 20-upleM′ ⊂ M,where(±) isthesignofthetop-leftentryofN(M′). Thisequation M′ is either trivial when M′ is linearly dependent, or reconstructs c˜(x ) otherwise. Condition 7 0 ensures that at least one subfamily M′ is linearly independent, so we sum the last equation over all subfamilies and arrive at the formula −1 ˜c(x0) =  (det(N(M′)))16 (±)M′N(M′). (14) M′⊂MX, #M′=20 M′⊂MX, #M′=20   Remark 4.1. The expliciteness of formula (14) is interesting in its own right and makes the stability of the problem straighforward to assess. On the other hand, the computation of 20×20 determinants can be expensive, and methods constructing a normal to a hyperplane without implementing (14) might reveal more practical. For a 20-uple M, an example of a potentially faster method for finding a scalar multiple of N(M) is to form a 20×21 matrix whose rows are the elements of M, and to find a vector in the nullspace of that matrix via Gaussian elimination. Now that the anisotropic structure C (or equivalently, c˜) is reconstructed, we now explain how the multiplicative scalar τ can be reconstructed via a standard transport equation. e 4.2 Reconstruction of the scalar factor τ Wenowswitchbackto4-tensornotation. PluggingthedecompositionC = τC intotheelasticity equation, where C is assumed to be known from the previous step, we obtain, for each elasticity e solution considered, e (C :ǫ)∇logτ = −div (C :ǫ), (15) e e where(C :ǫ) ∈ S (R)isnotnecessarilyinvertible. However, undertheassumptionthat{ǫ(j)}6 3 j=1 is a basis of S (R), by virtue of (5), so is {C : ǫ(j)}6 . Therefore, the identity tensor I e 3 j=1 3 decomposes into this basis by means of some functions µ (x) such that j e 6 I = µ (x) C : ǫ(j). 3 j Xj=1 e If we denote D := (C : ǫ(i)) : (C : ǫ(j)) for 1 ≤ i,j ≤ 6, the S (R)-valued function D = {D } ij 6 ij is known and invertible (as the Grammian matrix of the basis {C : ǫ(j)}6 ), and the entries of e e j=1 e 9 its inverse are denoted by Dij. In this case, the functions µ take the explicit form j 6 6 µ (x) = Dij I : (C :ǫ(i)) = Dij tr (C : ǫ(i)). j 3 Xi=1 Xi=1 e e Now taking a linear combination of (15) weighted by the functions µ , we obtain j 6 ∇logτ = − µ (x) div (C :ǫ(j)). (16) j Xj=1 e The right hand side of this equation is completely known and τ can therefore be reconstructed via either (i) integration of (16), or, after taking divergence, (ii) solving a Poisson equation on the domain Ω with known Neumann boundaryconditions. As discussed in [7], the approach (ii) maybemorerobusttonoisethantheformer,astheintegration ofordinarydifferentialequations in the presence of noisy measurements may strongly depend on the choice of integration path. Moreover, taking divergence of (16) has the advantage of naturally removing the curl part of noisy right-hand sides in (16). 4.3 Reconstruction of the tensor div C Now thatbothC andτ arereconstructed, weexplain howtogain stability onthereconstruction of the third-order tensor div C = ∂ C e ⊗e ⊗e . Note that this tensor is symmetric in the i ijkl j k l e pair of indices (k,l). Now the elasticity equation can be rewritten as (∂ C )ǫ +C ∂ ǫ = 0, 1 ≤ j ≤ 3, i ijkl kl ijkl i kl or in contracted notation, (div C) : ǫ = −C ∂ ǫ , 1 ≤ j ≤ 3. (17) j·· ijkl i kl For every 1 ≤ j ≤ 3, (div C) can beseen as an S (R)-valued function, andas such, usingagain j·· 3 the assumption that {ǫ(j)}6 is a basis of S (R), upon defining E = ǫ(i) : ǫ(j) for 1≤ i,j ≤ 6, j=1 3 ij the S (R)-valued function E = {E } is uniformly invertible over Ω and we denote by Eij the 6 ij components of its inverse. With such notation, Cramer’s rule in S (R) then reads 3 (div C) = Epq((div C) : ǫ(p))ǫ(q), 1≤ j ≤ 3, j·· j·· and combining this with (17), we obtain a reconstruction formula for the tensor div C (div C) = −Epq(C ∂ ǫ(p)) ǫ(q), 1≤ j ≤ 3. (18) j·· ijkl i kl 10