Stability in the linearized problem of quantitative elastography

The goal of quantitative elastography is to identify biomechanical parameters from interior displacement data, which are provided by other modalities, such as ultrasound or magnetic resonance imaging. In this paper, we analyze the stability of several linearized problems in quantitative elastography. Our method is based on the theory of redundant systems of linear partial differential equations. We analyze the ellipticity properties of the corresponding PDE systems augmented with the interior displacement data; we explicitly characterize the kernel of the forward operators and show injectivity for particular linearizations. Stability criteria can then be deduced. Our results show stability of shear modulus, pressure and density; they indicate that singular strain fields should be avoided, and show how additional measurements can help in ensuring stability.


Introduction
Elastography is a medical imaging technology; its current applications range from detection of cancer in the breast and in the prostate, liver cirrhosis and characterization of artherosclerotic plaque in hardened coronary vessels [19,42,54,4,53,56,15].
Elastography is based on the fact that tissue has high contrast in biomechanical quantities and the health state of organs is reflected in the elastic properties of tissue [22,2]. The most important of these is the shear modulus µ, which is the dominant factor in the propagation of shear waves in tissue; shear wave speed in tissue can change up to a factor of 4 with disease [47].
Elastography is performed by coupling with various established imaging techniques, such as ultrasound [27,41], MRI [38,29] or OCT [51,40] -What is common to these elastography techniques is that they provide interior data of the displacement u| Ω of the tissue on the imaging domain Ω. According to the specific excitation used, u can be space or space-time-dependent (both cases are treated in this work).
In some applications, knowledge of the displacement u already gives qualitative diagnostic information (see, e.g., [57] for a dermatological application). More accurate information is provided by quantitative information of the underlying biomechanical parameters. For this one needs to formulate the elasticity problem as a model; for elastography, various models based on linear elasticity, viscoelasticity or hyperelasticity have been considered [19].
To recover the material parameters, an inverse problem based on an elasticity model has to be solved. Given the displacement u, the mathematical problems in quantitative elastography are the recovery of parameters such as the Lamé parameter λ, the shear modulus µ, the density ρ, or recovery of the shear wave speed µ ρ [19,13]. These problems are non-linear inverse problems. Much effort in elastography research has been concerned with developing adequate numerical inversion schemes for all kind of experimental varieties of elastography (see, e.g., [30,31,32,6,39] and the reviews in [19,49]); among the proposed algorithms, optimization procedures, which need the linearized inverse problem, form an important class.
Mathematical results about uniqueness and stability of the inverse problems in elastography have been gathered recently, mostly for the simplest models in linear elasticity. In [20,34], it was proved, that -given one piece of dynamic interior information u subject to restrictions such as u = 0 -, one can uniquely recover material parameters. Other uniqueness results have been reported in [14] for two and more measurements, and in [46] for a hyperelastic model. The recent paper [10], showed unique reconstruction of λ and µ given two sets of exact measurements subject to some non-trivial conditions on the displacements. The stability of the nonlinear problem has been studied in [10] using ODE-based and variational tools.
Elastography can be seen as part of coupled-physics imaging methods [7]. This body of literature centers on novel imaging methods involving more than one physical modality; emphasis is laid upon the quantitative imaging problems for mechanical, optical or electrical parameters that have high diagnostic contrast. For reviews on the typical problems in coupled-physics imaging, see [5,8,23,55,7,24].
In many coupled-physics problems, high resolution in the reconstructions was observed. To explain this phenomenon in a unified manner, a general strategy is to linearize the corresponding nonlinear quantitative problems. Then the stability properties of the linearized problem have been treated with tools of linear PDE theory [26,9,12]. In [26], problems in conductivity imaging and quantitative photoacoustic imaging were treated. The linearized forward operators were studied with pseudodifferential theory and shown to be Fredholm, i.e., with stable inversion up to a finite-dimensional kernel. In [9], a general framework was proposed for treating linearizations of quantitative problems with interior data; in this method, the parameters are considered as variables and the information is recast into a single redundant PDE system, to which theory such as [50] is applied. This method was applied to the power density problem in [9] and to the problem of acousto-optic imaging in [12].
In this article, we treat the linearized problem in quantitative elastography using the general coupled-physics approach in [9], using [50]. This is the first time that this technique is applied to a problem in elasticity imaging. We treat stability and explicitly characterize the finite-dimensional kernel, and show injectivity in several versions of the linearized inverse problems.
The structure of the paper is as follows: in Section 2, we review the elasticity model which we are using, and in Section 3, we review necessary background of linear PDE theory. In Section 4, we apply this theory to the elasticity equation. We investigate how several kinds of linearizations perform analytically, derive the stability results and the characterization of the kernel. In Section 5, we discuss the stability conditions with respect to the literature. The appendix contains topological lemmas needed for the investigation of the kernel.

.1 Experiments and interior information in elastography
The general principle in elastography [19] is to • perturb the tissue using a suitable mechanical source • determine the internal tissue displacement using an ultrasound, magnetic resonance or optical displacement estimation method • infer the mechanical properties from the interior information, using a mechanical model Note that in elastography, there are two forms of excitation: an elastic deformation from the mechanical source, and the excitation from the ground modality. Also, the reconstruction procedure involves two steps: the recovery of the mechanical displacement u(x)| Ω resp. u(x, t)| Ω from measurements on the boundary, and the recovery of the mechanical parameters and properties from u| Ω . In this last step of quantitative reconstruction (which is treated in this article), the mechanical displacement u| Ω is also referred to as interior information.

A linear elasticity model for inhomogeneous linear isotropic media
There exist different variants of elastography using quasi-static, transient or time-harmonic mechanical excitations, but they can be described by common PDE models [42,19]. The elasticity models can be deduced from the equation of motion [28] ∇ · σ − ρu tt = F where F(x, t) is the excitation force density in N/m 3 . (We use the convention that letters printed in bold denote vectors in R 3 .) The mechanical displacement u(x, t) of the material point which is at position x at time t is measured in [u] = m. ρ(x) is the density in kg/m 3 . σ = (σ ij (x, t)) ij is the mechanical (or Cauchy) stress tensor with unit N/m 2 . Here, the divergence of a second-rank tensor A is computed column-wise: ∇ · A = ∇ · (a 1 , . . . , a n ) := (∇ · a 1 , . . . , ∇ · a n ).
The constitutive equation in elasticity (also termed Hooke's law) is with the stress tensor σ = σ ij (x, t) and the dimensionless strain tensor ε = ε kl (x, t). Here, * denotes the tensor multiplication. The material properties are incorporated into C = C ijkl (x), which is the rank-four stiffness tensor (with unit N/m 2 ). In most practical settings of elastography, one makes simplifications concerning the material parameters of the biological tissue [19,42]: Typically, the assumptions are that the material response is isotropic and linear.
Isotropy means that C is reduced to knowledge of two scalar quantites λ and µ, such that one has σ = 2µε + λ tr ε Id, instead of (2). Here, λ(x) is called the first Lamé parameter, and µ(x) is called the shear modulus or the second Lamé parameter. The physical units of λ and µ are the same as of σ and C, i.e. N/m 2 . Linear material behavior is encoded in the following representation of the strain ε: ε = 1 2 (∇u + ∇u ⊤ ).
Note that with this, the quantity tr ε in (3) is equal to ∇ · u. The equation of motion (1), with (3) and (4) is then augmented with boundary conditions and appropriate sources. There are different choices for initial and boundary conditions. One option is Existence, uniqueness and regularity properties for this model follow from the theory in [33] (for an earlier result for elastodynamic problems, see [17]). The initial and boundary values g and h are prescribed only for the complete the mathematical analysis. In practical physical experiments, the varying excitations enter in the source term F.
The standard problem of quantitative elastography is then to determine the material parameters λ(x), µ(x) and ρ(x) in the equation (5), given the interior information u(x, t)| Ω .

Adapting the quantities in the model for inversion
Models based on (3) and (5) are widely used in elastography for simulating the elastic behavior of tissue [42,19]. Nevertheless, the parameter which one reconstructs, is often only the shear modulus µ. Sometimes, one sets λ∇ · u = 0 or one assumes the incompressibility condition ∇ · u = 0. In these cases, λ does not occur in the model at all (see also the discussion in [29]).
We propose a different definition of quantities for the reconstruction. Precisely, we change the quantities and use the pressure p defined by With these quantities, it follows from (5) that The inverse problem is now to recover p(x, t), as well as µ(x) and ρ(x), given u(x, t)| Ω . It is this problem which we address in our stability analysis. Definition (6) has been used before, see e.g. [35,45]). Note that in tissue, one has that ∇ · u ≪ 1. Because of ill-posedness of differentiation, the quantity ∇ · u cannot be computed accurately from the data u in experiments. On the other hand, one has that λ ≫ 1. In numerical simulations, the pressure p turns out to be of order 1 and therefore should not be neglected [35].
Note that that p(x, t) is an elastic quantity, but not a material parameter: it depends on the particular displacement field induced by the excitation. Knowledge of p may or may not prove to be useful for diagnostic purposes. The reason, though, for introducing this quantity in the inversion model is that it numerically turned out to be useful. It was numerically observed that keeping p in the model improves the reconstruction of the shear modulus µ [35].
In our analysis, we will give a mathematical reason for using (7) instead of the first equation in 5. Before we come to that, we give the relevant background from PDE theory which we use in our work.

A result from linear PDE theory
We first treat the background from the general theory of linear PDE systems.
Let Ω be a bounded domain in R n (smoothness requirements on Ω are specified later). We consider the redundant system of linear partial differential equations for m unknown functions u(x) = (u 1 (x), . . . u m (x)), comprising in total M equations. Here, L(x, ∂ ∂x ) is a matrix differential operator of dimension M × m, For each 1 ≤ i ≤ M , 1 ≤ j ≤ m and for each point . Redundancy of the system means that there are possibly more equations than unknowns: M ≥ m.
Similarly, B(x, ∂ ∂x ) has entries B kj (x, ∂ ∂x ) for 1 ≤ k ≤ Q, 1 ≤ j ≤ m, consisting of Q equations at the boundary. The operations are again polynomial in the second variable. -S(x) is a vector of length M , and ϕ(x) is a vector of length Q.
We now define the notions of ellipticity and the principal part of L and B, respectively, in the sense of Douglis and Nirenberg [18]. Definition 1. Let integers s i , t j ∈ Z be given for each row 1 ≤ i ≤ M and column 1 ≤ j ≤ m with the following property: For s i + t j ≥ 0, the order of L ij does not exceed s i + t j . For s i + t j < 0, one has L ij = 0. Furthermore, the numbers are normalized so that for all i one has s i ≤ 0. Such numbers s i , t j are called Douglis-Nirenberg numbers.
The principal part of L for this choice of numbers s i , t j is defined as the matrix operator L 0 whose entries L 0,ij are composed of those terms in L ij which are exactly of order s i + t j .
The principal part B 0 of B is composed of the entries B 0,ij , which are composed of those terms in B kj which are exactly of order σ k + t j . The numbers σ k , 1 ≤ k ≤ Q are computed as where b kj denotes the order of B kj . Real directions ξ = 0 with rank L 0 (x, ıξ) < m are called characteristic directions of L at x. (The complex unit is denoted by the symbol ı = √ −1.) The operator L(x, ∂ ∂x ) is said to be overdetermined elliptic in Ω if for all x ∈ Ω and for all real vectors ξ = 0 one has that for the M × m matrix L 0 (x, ıξ).
To illustrate in an example, we can consider the system where Ω is the unit circle in R 2 . With ν, we denote the unit normal on ∂Ω. If we choose numbers (t j ) 2 j=1 = (1, 3), (s i ) 3 i=1 = (−1, 0, 0), we have the principal symbols and and (σ i ) 2 i=1 = (−1, −1). Note that the principal symbols differ in this case. Nevertheless, with both choices of numbers, L is overdetermined elliptic, as there exists a non-vanishing subdeterminant of L 0 in both cases.
Next we define the condition of B covering L, or the Lopatinskii boundary condition [50].
Definition 2. Fix y ∈ ∂Ω, and let ν be the inward unit normal vector at y. Let ζ be any non-zero tangential vector to Ω at y. Consider the half-line {y + z ν, z > 0} and the following system of ordinary differential equations on it: Consider the vector space of all solutionsũ of (15) which satisfyũ(z) → 0 for z → ∞. If this vector space consists just of the trivial solutionũ(z) ≡ 0, then the Lopatinskii condition is said to be fulfilled for the pair (L, B) at y, or B covers the operator L at y.
A typical example of an overdetermined elliptic systems with Lopatinskii boundary conditions is on a domain Ω with normal component u · v| ∂Ω given on the boundary. Another example of an overdetermined elliptic system with Lopatinskii boundary condition is the system of time-harmonic Maxwell's equations, where u = (H, E) satisfies with the normal component of H as well as the tangential component of E given on the boundary (see [50, §2] for (16) and (17)) In the context of hybrid imaging, examples of overdetermined elliptic systems with Lopatinskii boundary conditions have been considered in [9,12].
For investigating the stability for linearized quantitative elastography, we are going to use the a-priori estimate in [50] for the solutions of system (8). This theory does not need smooth coefficients, but coefficients in the Sobolev spaces W α p (Ω) (for the usual definition, also for noninteger values of α, see [1]). In the setting of [50] with Douglis-Nirenberg numbers t j , s i , σ k , one has that the operator A with acts on the space where l ≥ 0, p > 1. Under suitable restrictions on the coefficients L ij and B kj (specified below in the conditions of the theorem), the operator A is bounded with range in In formulating the restrictions on the coefficients of L and B, we simplify the version of [50, Thm. 1.1] for the following result. Theorem 1. Let integers l ≥ 0, p > 1 be given. Let (S, ϕ), the data from (8) be in R(p, l) as defined in (20). Let Douglis-Nirenberg numbers s i and t j be given for L in (8), and let σ k be as in Definition 1. Let Ω be a bounded domain with boundary in C l+max tj . Assume furthermore that p(l − s i ) > n and p(l − σ k ) > n for all i and k. Let the coefficients of L ij be in W l−si p (Ω), and let the coefficients . Then the following statements are equivalent: (8) is overdetermined elliptic (see (11)) and the Lopatinskii covering condition (15) is fulfilled for (L, B) on ∂Ω.
2. There exists a left regularizer R for the operator A = L × B in (18), that is, we have

The following a-priori estimate holds
The assertion of the theorem gives a criterion for the existence of a left regularizer for the overdetermined redundant systems. For the case of boundary value problems for square systems with M = m, such an equivalence is established in the classical work of [3]. For square systems, one has the stronger statement that ellipticity and Lopatinskii condition are equivalent to the Fredholm property of a differential operator (which also needs the existence of a right regularizer). -The criterion for redundant systems with M ≥ m was established in [50], and investigates the stability estimate, and even gives a representation formula for the solution, provided it exists. The existence of a right regularizer Q with AQ = I − T (which would yield local existence) cannot be assured in general for overdetermined systems.
We will exploit this criterion for the linearized version of quantitative elastography.

Stability analysis 4.1 Setting and notation
For treating the hybrid imaging problem described in 2.1, we take the adapted forward model (5) + (6). We recast the equations of the forward problem for the displacement, and the interior information from the measurements, respectively, into a single system of partial differential equations: Here, we are given F (the excitation force), g and h, as well as the interior information K. We formally keep u and K distinct in the second equation, as u is treated as variable of the system, and K represents the data. -We aim at a quantitative estimate such as (23) with the measurement data in the inhomogeneity.
We consider the system (24), for the variables u, p, µ, ρ; therefore, we consider it as nonlinear, involving multiplication of the unknowns. In order to make it tractable for the analysis, we linearize this system to provide equations of the form (8). Several of these linearizations will be considered in the subsequent theory.
We first consider the parameter-to-solution operator which maps the choice of parameters (p, µ, ρ) to the displacement field u satisfying (24). Then we consider the linearization of V at a reference state (p, µ, ρ), By formal differentiation of (24), we find that, at the reference state u = V(p, µ, ρ), the increment δu satisfies the equations Note that F does not depend on the reference state, therefore no inhomogeneity appears in the first equation in (26).
We write and we introduce the operator (29) The operator L pµρ in (29) is the linearization of the redundant system (24) with respect to p, µ and ρ. Note that L pµρ is a matrix differential operator like L in (8).
Apart from L pµρ , we also introduce the operators corresponding to directional derivatives with respect to only one parameter: We also use the combination In the stability analysis, we will make a comparison of the properties of these operators.
In the model (24) which we started from, we incorporated the definition of the pressure in (6). An alternative is to use the original model (5) only, which involved the first Lamé parameter λ. This corresponds to re-substituting p = λ∇ · u in (24). In exact analogy to constructing L pµρ one can form the forward operator V λ . One then considers its linearization V ′ λ (λ, µ, ρ) : (δλ, δµ, δρ) → δu, and introduces the operator L λµρ . We particularly will consider Up to now, we defined several differential operators L of form (9). As in (18), we now combine them with boundary data and then form equations Aw = (S, 0) as in (8) resp. (21). The vector w changes according to the variables which are in the system, specified below. For the inhomogeneity S, we have We introduce the operator from which we get the following specializations.
The system has 6 interior equations for the unknowns (δρ, δu 1 , δu 2 , δu 3 ). All of these are linear differential systems which are redundant systems of form (8) with M ≥ m. Therefore we can apply the methodology of Section 3 to these.
Up to now, we only incorporated information from one imaging experiment in our operators. It is possible, though, to conduct more than one imaging experiment, and the consideration of multiple measurements in the inverse problem is typical in hybrid imaging. Therefore, we can use different excitations F i and possibly different functions g i , h i in (24) and obtain different versions of u i (x, t) and interior information K i (x, t). For each experiment, we also have a different variable p i in the system. While these quantities change with each excitation, the material parameters λ, µ and ρ remain the same.
For example, we write for the system corresponding to 2 experiments. Here, the operators L (i) In the inhomogeneity, we have the quantities S (i) = (0, 0, 0, δK µ in (41) shows the effect of adding one more experiment in the system: there are 6 more equations and 3 new variables: together 12 equations and 7 unknowns.
The shown procedure of addition addition experiments can be applied to any of the operators in (38), (39), (36). For each experiment we add, the inequality M ≥ m in the nomenclature of Section 3 is fulfilled and the system is redundant.
For reasons which are apparent later, we can augment the boundary operator with additional constraints and use Normally, we will use (28) as boundary operator. The conditions (42) will be used in special cases which are separately indicated.
Note that in this section, we have given the general form of the linearization operators for the dynamic case on a cylindrical domain Ω × T . Sometimes, we will consider these operators in the quasi-static case with using only the spatial domain Ω. This is specially indicated in each case.

Ellipticity
We want to apply the methodology of Section 3, and use the criterion in Theorem 1. Therefore, we have to we determine the ellipticity condition in Definition 1 for the operators L p , L µ , L ρ in (30), (31), (32). We first determine the principal symbol and possible characteristic directions for the operator L pµρ in (29), which is treated in Proposition 1. From this analysis we then draw some corollaries concerning the ellipticity of L p , L µ , L ρ , as well as ellipticity of L pµ . For the analysis of L pµρ , we choose Douglis-Nirenberg numbers (t j ) 6 j=1 = (1, 1, 0, 2, 2, 2) and corresponding to the variables (δp, δµ, δρ, δu 1 , δu 2 , δu 3 ) and numbers (s i ) 6 i=1 = (0, 0, 0, −2, −2, −2) corresponding to the six equations. If there are less variables in the system (as in L p , L µ , L ρ , L pµ ), then only the corresponding Douglis-Nirenberg numbers are used (e.g., for the analysis of L p , we have the numbers (1, 2, 2, 2) for the variables (δp, δu 1 , δu 2 , δu 3 ) in L p ). a) The principal symbol of L pµρ (in the dynamic case) is where ξ s := (ξ 1 , ξ 2 , ξ 3 ) for ξ ∈ R 4 . b) For every point (x, t), the operator L pµρ is not overdetermined elliptic.
c) In the quasi-static case (43), the principal symbol of L pµρ = L pµ is for ξ ∈ R 3 . d) For every point x consider L = L pµρ = L pµ in the quasi-static case (43). Then L is not overdetermined elliptic at x.
Proof. a) To compute the principal symbol (44), we refer to the definition of L pµρ in (29). We write the first three equations in L pµρ (δp, δµ, δρ, δu) = S as where we denote the columns of the (symmetric) strain as The meaning of ε(δu) i is analogous. The last three equations in L pµρ (δp, δµ, δρ, δu) = S are written as We now want to determine the entries L 0,ij of the principal symbol L 0 ((x, t), ıξ). Note that each of the columns in L 0 exactly corresponds to one of the variables (δp, δµ, δρ, δu 1 , δu 2 , δu 3 ). Starting with j = 0, we go through the variable list until j = 6. For each variable, we determine the term where the unknown appears in the equations (46) resp. in (47). Then we choose that component of the term which has order t j + s i . Substituting ıξ for ∂ ∂(x,t) in that component gives the corresponding entries (L 0,ij ) 6 i=1 in the j-th column of L 0 . For the first column corresponding to δp, the term ∂ i δp in (46) The second column corresponding to δµ, and the summand of highest order in the term 2 In the third column corresponding to δρ, no differentiation occurs, so we just have −(u i ) tt as entries in (L 0,i3 ) 3 i=1 . The last three columns correspond to the variables (δu 1 , δu 2 , δu 3 ). The relevant terms in (46) are We substitute iξ for differentiation in (48) and take the terms of highest order to find the entries of the columns (L 0 ) i,j , j = 4, 5, 6; these are the terms −µξ 2 i − µ|ξ| 2 +ρξ 2 4 in (L 0 ) i,i+3 , i = 1, 2, 3, and the term −µξ i ξ j in the entries (L 0 ) i,j , j = 4, 5, 6, j = i. -The last three equations (47) contain no derivatives and give rise to the identity matrix in the entries (L 0 ) ij , 4 ≤ i, j ≤ 6 of the principal symbol.
c) The calculations for the symbol (45) in the spatial case are exactly the same as for the spatio-temporal case. The only change in this case is that there is no temporal derivative. Therefore there is no variable δρ, and one column less than in (44), and ξ 4 can be set to zero everywhere. b) and d) We first observe that both in (44) as well in (45), three of the columns are clearly linearly independent. The first two columns, though, can be linearly dependent: Consider the symmetric strain ε = ε ⊤ . As the entries are real, there always exists an eigenvector v such that where * denotes matrix multiplication. Choosing (ξ 1 , ξ 2 , ξ 3 ) = v gives linear dependence of L 0 (ıξ) in the first two columns.
In conclusion, at each point (x, t) resp. x, there are choices of ξ such that the symbols in (44) resp. (45) do not have full rank. Therefore L pµρ is not overdetermined elliptic.
We now restrict the focus on linearizations in only one direction, which were introduced in Section 4.1. Then the corresponding principal symbol contains fewer columns and results on ellipticity can be obtained. Corollary 1. The operator L p in (30), considered in the stationary case, is elliptic everywhere.
Proof. Let L = L p . The principal symbol consists of the first and the three last columns of the matrix in (45): For ξ = 0, this symbol has maximal rank 4 everywhere, therefore L p is elliptic.
The operator L s µ , corresponding to s measurements (δu) k , 1 ≤ k ≤ s, is elliptic exactly at points where at least one of det(ε(u k (x, t)) = 0.
Proof. Let L = L µ (the case of one measurement). Then the principal symbol consists of the second and the three last columns of the matrix in (45): This symbol has rank 4 provided that the first column is non-degenerate, which is equivalent to the condition (50) of non-singular strain. Now let L = L In the first column of L s 0 , there are three new entries 2ıξ s · ε(u s ) 3 i=1 at position (L 0,i1 ) 3s i=3s−2 . If we have that for 1 ≤ k ≤ s, at least one of det(ε(u k (x, t)) = 0 is non-zero, then the principal symbol L s 0 has full rank: in the case k < s because of the induction assumption, and in the case k = s because of nondegeneracy in the first column due to the new measurement. is elliptic exactly at points (x, t) where at least one of (u k ) tt (x, t) = 0, 1 ≤ k ≤ s.
Proof. Let L = L ρ . The principal symbol consists of the last four columns of the matrix in (44): This symbol has rank 4 iff u tt is nonzero. The statement for multiple measurements is proved by induction, analogous to the case of L s µ .
Corollary 4. The operator L pµ is not elliptic at points (x, t) ∈ Ω × T in the spatio-temporal case, resp. points x ∈ Ω in the spatial case.
Proof. Let L = L pµ in the dynamic case. The principal symbol consists of the two first and the three last columns of the matrix in (44): Fix a point (x, t). As in the proof of Proposition 2c), choose an eigenvector v of ε and let ξ such that (ξ 1 , ξ 2 , ξ 3 ) = v. With this choice of ξ, the principal symbol (53) has two columns which are linearly dependent. Therefore, the operator L pµ is not overdetermined elliptic.
Remark 1. As starting point of our analysis, we have choosen the modified model (7) with the substitution p = λ∇ · u in (6). Of course, we could also analyze the system corresponding to the original model (5). Using linearization in direction of δλ, we have the operator L λ in (34). Now let L = L λ in the quasi-static case. Doing calculations as in Proposition 1, and restricting the focus on the first and last three columns as in Corollary 1, one finds the resulting principal symbol as In this case, ellipticity is harder to obtain, since it is tied up with non-vanishing of ∇ · u. This actually may be a condition which is harder to guarantee in experiments, see the discussion in Section 5.

Lopatinskii condition
We want to use the stability criterion in Theorem 1 for the systems A p , A µ and A ρ in (37), (38), (39). For this purpose, we are checking the covering condition in Definition 2 for the various differential operators L introduced in Section 4.1, together with the relevant boundary data where necessary.
Proposition 2. The systems L p and L µ in (30) and (31) satisfy the Lopatinskii condition with arbitrary boundary data.
Proof. For checking the Lopatinskii (or covering) condition in Definition 2, we have to consider the vector space of functions satisfying the system (15) and show that it is trivial. Let L = L p with principal symbol (49), and let y be a point on the boundary. Then the system of equations L 0 (y, ıζ + ν d dz )ũ(z) = 0 in (15) comprises 6 equations for 4 unknowns (ũ 1 (z),ũ 2 (z),ũ 3 (z),ũ 4 (z)).
In the entries (L 0 ) ij , 4 ≤ i ≤ 6, 2 ≤ j ≤ 4 of (49), there is a 3-by-3 identity matrix. The three last equations of (15) therefore mean that The first three equations in (15) then reduce to Any functions which are in the vector space considered in Definition 2 have to satisfy (56). The only possible solutions of (56) consist of functions of form e ıλz , where the parameter λ = ζi νi is a real number; neither of these functions tends to 0 for z → ∞. In this argument, we did not use any boundary constraint. Therefore the vector space to be considered in Definition 2 is trivial for every y ∈ ∂Ω, and the Lopatinskii covering condition is always satisfied. Now let L = L µ with principal symbol (51). Then, for y on the boundary, consider the system L 0 (y, ıζ + ν d dz )ũ(z) = 0 in (15). Similarly to (55), the last three components vanish, and the system reduces to Here, we use fixed vectors g j = ε(u(y)) j , j = 1, 2, 3 and ζ s = (ζ 1 , ζ 2 , ζ 3 ).
Suppose that the relation g j · ν = 0 holds for all j. Then we have solutions of (57) of form e ı g i · ζ s g i · ν z . The numbers gi · ζ s gi · ν are real, so neither of these functions tends to 0.
Assume, on the contrary, that all entries in A vanish: Incorporating the information (66)-(68), the nine equations in (70) can be written in matrix form as Here, * denotes matrix multiplication.
Using the definition of g j in (62), it follows that ε(u(y)) * ν(y) = κν(y), as we have that ε = ε ⊤ . But, by assumption, ν(y) must not be an eigenvector of the strain ε(u(y)), so this cannot happen. Therefore (70) is wrong and A = 0. Consequently, the equations (63)-(65) are always nontrivial in our case. As these are linear ordinary differential equations of second order, the basis of solutions consists of functions of form e λz and z e λz .
We now claim that there can be no solutions of form z e λz to the differential equations (59)-(61).
Assume, on the contrary, that a solution of form z e λz exists. According to the theory of linear ODE, λ is a double zero of the three characteristic polyomials The discriminant has to vanish, so we have λ = (2,3); consequently we also have c k = 4λ 2 a k . Because of these proportionalities between a k , b k and c k , the matrix A in (69) has rank 1. Therefore we have also a proportionality between the rows in A, i.e. the relations for some constants γ, δ. Using (66)-(68), the equations in (73) can be written in matrix form as  which we interpret, as in (71), as linear system for the unknown variables {g 1 · ζ s , g 2 · ζ s , g 3 · ζ s , g 1 · ν, g 2 · ν, g 3 · ν}. Now the matrix in (74) has rank 6, so the equations reduce to g j · ζ s = 0 j = 1, 2, 3 and g j · ν = 0 j = 1, 2, 3.
As in (72), it follows from the second equation in (75) that ν is an eigenvector of the strain ε(u(y)) corresponding to the eigenvalue 0. But this has been ruled out by hypothesis.
So the only solutions of (59)-(61) are of form e λz . The boundary condition (42) then leads to the unique solution ofũ 1 (z) = 0. The same chain of arguments can be invoked to obtainũ 2 (z) = 0. Therefore the Lopatinskii condition is satisfied at y.

Stability
In Sections 4.2 and 4.3, we collected results about the systems A p , A µ and A ρ defined in (37), (38), (39); this now allows to derive the following information about the stability properties of these operators, and about reconstruction of parameters in the linearized inverse problem of quantitative elastography.
We choose numbers p = 2, l with l > 0, pl > n, and a bounded and connected domain Ω with C l+2 -boundary.
In the following statements, we suppose that we have a reference state (p, µ, ρ) with µ ∈ H l+1 (Ω), ρ ∈ H l (Ω) and for the reference displacement u = V(p, µ, ρ) in (25) we require u ∈ H l+2 (Ω), u tt ∈ H l (Ω). The existence of such a displacement field u can be ensured by [33,Thm.8.1] if we require that the material parameter are in a Hölder space.
In the following theorem, we have recourse to the notion of an operator having a left regularizer, as in Theorem 1, see also [21]. Other names for this notion are that A is a left semi-Fredholm operator [16, Ch. XI, §2].
Theorem 2. The operator A (r) µ for r ≥ 1 measurements, described in (41) and considered in the quasi-static case, has a left regularizer on Ω precisely when, for all x ∈ Ω, det(ε(u k (x))) = 0 for at least one measurement 1 ≤ k ≤ r.
Then the stability estimate holds with a finite-dimensional kernel K 1 .
Proof. We first treat the case of A µ .
Existence of the left regularizer in (78) implies that ran(A µ ) is closed [16, Ch.XI, Thm. 2.3(ii)]. Therefore ran(A µ ) is a Hilbert space. We apply the open mapping theorem [16, Ch.III, Thm. 12.1] to find that the inverse A −1 µ : ran(A µ ) → D(p, l)/K is continuous. Therefore, we have the existence of a real number C such that for all (δµ, δu) ∈ D(p, l), the estimate holds. By (26), we have δK = δu in S. We set K = K 1 ×K 2 with K 1 ⊂ H l+1 (Ω). Then from (79), we obtain the estimate (77) for r = 1. The case of A (r) µ follows straight-forwardly by induction.
Remark 2. Theorem 2 gives the criterion (76), which means that at each point, at least one of the measured elastic displacement fields has non-singular strain.
The requirement of such qualitative conditions for the solutions is typical for the coupled-physics literature. In fact, the condition (76) for r = 2 is a generalization of the invertibility condition for the nonlinear reconstruction problem in elastography, which was found in the research of [10], namely Here, we have for k = 1, 2 that t k := tr(ε(u k )) = ∇ · u k and that ε(u k ) D := ε(u k ) − t k 3 Id is the deviatoric part of the strain. It can be verified by simple calculation that violation of (76) leads to violation of (80).
It is unknown whether (80) can be ensured with two vector fields for every distribution of material parameters. For several parameter classes, existence of boundary conditions ensuring (80) can be justified, see the discussion and examples in [10,Sec.3.3]. As (76) is a consequence of (80), the special argumentation for (80) can also be invoked for arguing for the premise of Theorem 2 in our case.
We now give an explicit characterization of the kernel of the operator A (1) µ , which will be exploited in Corollary 5 to show injectivity of A (2) µ , the operator corresponding to two measurements. with fixed p ∈ Ω. Here, the vector field a(x) is uniquely determined by where u is a reference state for which (76) holds.
Proof. In the proof of the statement, we derive a representation for δµ on a connected set and then infer that the representation is valid on Ω by a topological argument. Suppose, to begin with, that (δµ, δu) ∈ D(p, l) is in the kernel of A µ = L µ × B.
If δµ ≡ 0, then the assertion is trivially satisfied. Otherwise, there exists a point p ∈ A. In this case, consider the connected component V of p in the topology of A ⊂ Ω, that is Suppose now x ∈ V . We analyze the 6 equations From the last three equations, we immediately get δu| V = 0. From the first three equations, we then get that ∇ · (δµ ε(u)) = 0 on V for the element δµ, and hence Evaluating (84) at the point x ∈ V ⊂ A and dividing by δµ(x) (which, by (83), is non-zero) shows that with a determined by (82). Actually, the conditions (76) and (82) can be used to define a(x) uniquely for x ∈ Ω. To see this, set b i := ε(u) i . By (76), the vectors b i form a basis of R 3 . Then write Here, the vectors b ′ i are obtained by Gram-Schmidt orthogonalization, i.e., b ′ 1 := b 1 , and The scalar products in (86) can be represented as In (82), the scalar products a · b i are specified for i = 1, 2, 3. Now using (88), we see by a simple induction argument that (82) determines a · b ′ i in (86); thus we can determine a uniquely on Ω.
Furthermore, observe that u ∈ (H l+2 (Ω)) 3 implies b i ∈ (H l+1 (Ω)) 3 and Using this, as well as (86), (88), the inequality pl > n and the Sobolev embedding theorem [1,Thm.5.4.C], we obtain that a is continuous on Ω. Now consider the vector field a and calculate the path integral from p to x to find x p a(y)dy From this identity, we have the representation for the values δµ on the set V ⊂ A ⊂ Ω. Note that the function on the right hand side of (89) is continuous and defined on the whole domain Ω. We now claim that actually, we have such that the representation formula (89) holds for x ∈ Ω. Assume, on the contrary, that V Ω. We then also have that A Ω so Ω would not be connected). Therefore, the assumptions of Lemma 2 are satisfied. Consequently, there exists a point q ∈ ∂V \ A, where ∂V is the boundary of V in Ω. As q ∈ A, we have, by (83), that δµ(q) = 0.
As q ∈ ∂V , there exists a sequence v n ∈ V with v n → q in Ω.
By the representation (89), together with (92), we have On the other hand, continuity of δµ, equation (91) and (92) imply that As (93) and (94) contradict each other, we infer that the assumption V Ω is wrong. Therefore, as asserted in (90), V = Ω holds.
such that together with the boundary data we have Suppose that δµ = 0. Then there exists a point p ∈ Ω with δµ(p) = 0. As in the proof of Theorem 3, where we derived the representation formula (89) for x ∈ Ω, there exists a certain vector field a such that This implies that δµ(x) > 0 for all x ∈ Ω. Therefore, the condition ess inf Ω µ = ess inf Ω µ 0 > 0 in [33, (2.2)] is satisfied. The uniqueness result [33, Thm.5.2] then implies that, from (95), we have that u 1 = u 2 . But this is contradiction to our assumption. Therefore, we have ker(A In the subsequent part of the section we give the stability criteria for the operators A ρ , A p and A λ .
One has the stability estimate Proof. We first treat A ρ . The case of A (r) ρ follows by induction. The stability criterion in Theorem 1 is established for domains with C l+max tj boundary. Upon careful checking of the proof [50, §6], the only place where this assumption enters is the existence of a partition of unity. Now our domain is Ω×[0, T ], The construction of a partition of unity easily generalizes to cylindrical domains Ω × [0, T ], where Ω has C l+max tj boundary. Therefore, we can apply Theorem 1 to the problems with cylindrical domains.
The ellipticity condition has been assured in Corollary 3, and the Lopatinskii condition is satisfied according to Proposition 3. The assumptions in these results give the requirement u tt (x, t) W = 0. With that, the equivalent conditions of Theorem 1 are fulfilled and we apply the result as in the proof of Theorem 2.
Theorem 5. The operator A p in (37), considered in the stationary case, has a left regularizer on Ω, and we have the estimate Here, the kernel K 3 consists of the (one-dimensional) space of constant functions on Ω.
Note that, with the same method, but using Remark 1, one obtains a conditional stability result for the operator A λ : Theorem 6. The operator A (r) λ corresponding to r measurements, considered in the stationary case, has a left regularizer on Ω precisely when, for all x ∈ Ω, Then the stability estimate holds for a finite-dimensional kernel K 4 . In the research for coupled-physics conductivity problems, ellipticity has been investigated theoretically and numerically, and found to yield optimal stability estimates, avoid blurring effects, accurate reconstruction of edges, and absence of propagation of singularities [26,23,25,9,36,11].
Note that failure of ellipticity in our cases entails non-existence of a left regularizer non-existence of a left regularizer is equivalent to either dim ker(A) = ∞ or the range of A not being closed for the particular Sobolev spaces involved [16,XI,Thm. 2.3]. This does not mean that necessarily, the linearized problem will be unstable for all data in any function space. For example, consider the case of Corollary 2: at a point x, there might be just one direction ξ for which ellipticity does not hold. Then one can form the conjecture that reconstruction can still be stable if there is no edge along this direction (see the related discussion in [26, 6(ii)]). We plan to address this in future work.
2. The ellipticity conditions for A λ and A ρ seem to be natural. Concerning λ, literature actually often assumes the incompressibility condition ∇ · u = 0 on the whole of Ω [45,6,29,14]. In this case, of course, the measurement data are not dependent on λ, so this parameter cannot be reconstructed then. -But in the compressible case, where ∇ · u = 0 on the whole of Ω, there still might be single points x at which ∇ · u(x) = 0. Notice that, as stated in Remark 1, the ellipticity analysis along the lines of this article then entails that at such points x, ellipticity is lost and every direction is a characteristic. -Concerning the particular data one has, it might then be better to reconstruct the pressure p = λ∇ · u, with the operator A p being always elliptic.
Similarly for ρ: If u tt = 0 on the whole of Ω, the parameter ρ does not appear in the model, so it cannot be reconstructed from the measurements. If, on the other hand, u tt (x) = 0 only for particular points x, the analysis says that ellipticity is lost at these points x, and every direction is a characteristic for A ρ there.
3. The ellipticity condition for reconstruction of µ turned out to be the nonsingular strain condition in (76), which is a generalization of the condition in (80). Apart from this characterization, points of singular strain have been found in experiments, namely at the intersection of nodal lines or surfaces in early experiments of elastography using eigenmodes (see [44,43,52]). Empirically, it was observed that these patterns could be avoided by choosing multi-frequency excitation functions F [44].

Conclusion
We have applied a general method of linear PDE to linearized problems in quantitative elastography in R 3 , with interior data given. We analyzed ellipticity conditions of the PDE problem augmented with the interior data. We deduced simple criteria for the stability of the linearization. This analysis revealed stable reconstruction of the shear modulus µ and the hydrostatic pressure p = λ∇ · u, but pointed to a difficulty of reconstruction of λ. For the reconstruction of µ and ρ, the kernel in the linearization was shown to be trivial for choice of two measurements. The results give a mathematical explanation which biomechanical parameters can be stably reconstructed from interior measurement data u.
Among all subsets of A which are connected and contain p, the component V is maximal. Therefore p ∈ V ∪ U 1 = V , or equivalently U 1 ⊂ V . This shows that V is open in Ω. Proof. We use Lemma 1 and prove the statement in two steps: first, we find a point q ∈ ∂V \ V ; second, we show that q ∈ A.
Claim 1: There exists a point q ∈ ∂V \ V . We have that V ⊂ A Ω. Therefore, there exists an element y ∈ Ω \ V.
Consider the mapping f : V → R v → |y − v|.
Observe that V ⊂ Ω is closed and bounded, hence a compact set; observe also that f is continuous. Therefore, a minimum exists, that is: We now show that, actually, the point q ∈ V = V ∪ ∂V is not contained in V . Once this is shown, Claim 1 is proven. Assume, on the contrary, that q ∈ V . According to Lemma 1, we then would have an ε, such that U 2 := {z : |q − z| < ε} ⊂ V.
Without loss of generality, we can assume ε < 2. Now, using the element y from (102), define the point Then calculate This would contradict (103). -Therefore, q ∈ V .
Claim 2: The point q in (103) does not belong to A. We prove this claim indirectly. Assume that q ∈ A. (104) Recall that, according to Claim 1, q ∈ ∂V , where ∂V is the boundary of V in Ω.
Hence there exists a sequence v n ∈ V with v n → q in the topology of Ω. We assert that v n → q in the topology of A.
To see this, choose an open set U 3 ⊂ A with q ∈ U 3 . Because A is open in Ω, U 3 is open in Ω as well. Now the elements v n converge to q in Ω; therefore, there exists an N , such that for all n ≥ N : v n ∈ U 3 ; hence we have (105). The set V , which is the connected component of the point p, is closed in the topology of A [37, Thm. 23.4]. But a closed set contains all its limit points. Therefore, with (105), we would have that the limit of the sequence v n ∈ V lies in V , so q ∈ V . But this is a contradiction to Claim 1. -Therefore, contrary to (104), we have q ∈ A.