Universality of the angled shear wave identity in soft viscous solids

Mechanical stress within biological tissue can indicate an anomaly, or can be vital of its function, such as stresses in arteries. Measuring these stresses in tissue is challenging due to the complex, and often unknown, nature of the material properties. Recently, a method called the angled shear wave identity was proposed to predict the stress by measuring the speed of two small amplitude shear waves. The method does not require prior knowledge of the material's constitutive law, making it ideal for complex biological tissues. We extend this method, and the underlying identity, to include viscous dissipation, which can be significant for biological tissues. To generalise the identity, we consider soft viscoelastic solids described by a generalised Newtonian viscous stress, and account for transverse isotropy, a feature that is common in muscle tissue, for instance. We then derive the dispersion relationship for small-amplitude shear waves superimposed on a large static deformation. Similarly to the elastic case, the identity is recovered when the stress in the material is coaxial with the transverse anisotropy. A key result in this paper is that to predict the stress in a viscous material one would need to measure the wave attenuation as well as the wave speed. The case of viscoelastic materials with memory is also discussed.


Introduction
The function of most materials, and durability, is intimately linked to the level of mechanical stress within the material.For instance, in civil engineering, structures are designed so that internal forces remain below the elastic limit (that is, the yield stress), while also counter-acting the effects of gravity [1].Similarly, mechanical stress plays a major role in biomechanics, such as in the study of bone fracture [2].
Despite its importance, measuring the internal stresses, in a non-destructive way, has proven very challenging.Early attempts to determine residual stress in metals using ultrasonic waves rely on the acousto-elasticity theory, which connects changes in the wave speed to a level of applied stress and some of the material parameters [3].This approach can be adapted to modern imaging techniques based on ultrasound, see for instance Bied and Gennisson [4].However, using these techniques to measure the stress has proven difficult, as they require prior knowledge of the material parameters [5].
To develop a way to measure the stress, without needing to know the material constants, we can use universal relations [6].In mechanics, universal relations are mathematical formulas that connect a material's stress to deformation and motion, without involving the material's parameters explicitly.Instead, these relations hold for all materials within a certain class [7].For this reason, universal relations are useful for the design of experimental procedures aimed at testing a large range of materials.
In terms of wave propagation, there exists universal relations that connect the wave speeds to the stress in the Figure 1: Sketch of the angled shear wave identity in a rectangular cuboid maintained in a state of static stress (sectional view for constant ).The selected shear waves propagate along two directions with in-plane polarisation.
material, regardless of the material parameters.Examples include Ericksen's identities providing the speed of compression and shear waves in elastic solids [6], and an identity by Man and Lu [8], which links the stress to the speed of four different shear waves.In both cases, the superimposed waves are assumed non dispersive, i.e., their speed does not depend on frequency.Therefore, these identities do not apply to many soft solids which are inherently viscoelastic, hence dispersive [9,10].Named 'angled shear wave identity', one simple universal relation that connects stress to wave speeds takes the form [5] where the mass density  is a constant.This identity connects the difference σ1 − σ3 between the normal stresses along  and  to the speed ,  ′ of only two shear waves, see Figure 1.These waves are both polarised in the same plane (  ,   ) and propagate at specific angles  and  ′ = ∕2 ± .Here too, the superimposed waves were assumed non dispersive.In Zhang et al. [11], the method was applied to the measurement of stress in hydrogel and in skeletal muscle, without the knowledge of the constitutive parameters of the materials.In the supplementary material, the authors study the applicability of the method to transversely isotropic elastic media such as muscle tissue [12], whose anisotropy derives from the orientation of its fibres.Zhang et al. [11] prove that the universal relation still holds as long as the fibres are aligned with one of the principal stresses.Similarly Mukherjee et al. [13] demonstrated that the identity (1) holds for compressible elastic solids, with texture anisotropy due to fibres or other structural components, if the texture anisotropy is aligned with the initial stress, and if both the texture anisotropy and stress are small (in comparison to the Lamé constants of the material).
The main aim of this work is to extend this angled shear wave identity to a larger class of materials that are dispersive and incompressible.Examples of materials include soft solids such as hydrogels and biological tissues.The theory of non-linear viscoelasticity offers a highly suitable modelling framework to achieve this, cf.Destrade et al. [14].In this framework, harmonic plane waves no longer propagate with constant waveform and speed.Instead, due to dispersion, waves are attenuated as they propagate, and their speed depends on the frequency [15].Since the above universal relations were obtained in the elastic case, it remains unclear how they apply to viscoelastic materials.
In the present letter, we consider a general model of soft viscoelastic solid, and the particular case of Kelvin-Voigt solids with Newtonian viscosity (Section 2).These theories account for incompressibility, nonlinear elasticity, viscosity, and transverse isotropy.To analyse the properties of mechanical waves that propagate in such a material subjected to a homogeneous pre-stress, we then proceed with an acousto-elastic linearisation, or small-on-large (Section 3).The results are used to prove that the angled shear wave identity (1) still holds in viscoelastic solids, as long as the phase velocity is replaced by a complex wave speed which accounts for attenuation (Section 4).To reach this result, we made the following assumptions: 1. the material is incompressible, viscoelastic of differential type, and transversely isotropic; 2. there exists a dissipation potential from which the viscous stress contribution is derived; 3. the material is subjected to tri-axial stress, which is aligned with the directions of the transverse anisotropy.
Generalisation to isotropic materials with memory is addressed in the Appendix A.

Constitutive theory
Let us introduce the deformation gradient tensor  = ∕ which links the current position of a point on a material , after a deformation, to its initial position  in the reference configuration, before the deformation.The material is assumed incompressible, which means that the volume does not change during any deformation.Thus, the determinant of  satisfies det  = 1, and the mass density  is constant [9].
According to our assumptions, the Cauchy stress tensor  can depend on the elastic deformation, deformation rate, and direction of anisotropy.Here, we decompose the stress as the sum of an elastic contribution  e and of a viscous one,  v , that is, where  is the identity tensor.The stress − is a reaction to enforce the incompressibility constraint.The conservation of angular momentum requires that  is symmetric, i.e., equal to its transpose  T .The expression of the elastic stress is deduced from the strain energy function  e as follows [9], where  e depends on  only through , and are the right and left Cauchy-Green deformation tensors.Due to incompressibility, the symmetric tensors ( 4) have unit determinant.
To account for anisotropic materials, we consider a unit vector  which defines the direction of anisotropy in the reference configuration, while  =   is a vector in the direction of anisotropy in the current configuration.Now by assuming that  e depends on only  and , we have that  e ∶=  e (, ), (5) where the term  =  ⊗  is used instead of  so that swapping  with − does not change the energy  e .Likewise, we define  =  ⊗ .
From the assumption (5), it can be deduced that  e depends on only ten scalar invariants of the tensors ,  [16,17].When further using det  = 1 and the fact that  is a unit vector, these invariants reduce to [12]  1 = tr(),  2 = tr( −1 ), The first two scalar functions   are the principal invariants of the deformation (4), whereas the latter ones are related to anisotropy.
To describe the stress due to viscosity ( v ), we need to introduce the rate of strain Ċ and the symmetric part of the Eulerian velocity gradient  = 1  2  −T Ċ −1 , where the overdot denotes the material time derivative.Due to incompressibility, det  = 1, and as a consequence [14], we have tr  = 0.
The viscous stress is deduced from the dissipation potential  v such that where  v can depend on the elastic deformation, deformation rate, and direction of anisotropy , in the form [18] The viscous stress (7) vanishes in equilibrium ( Ċ = ), where the total stress ( 2) is purely elastic.
From the form (8), it can be deduced that  v depends only on seventeen scalar invariants of the tensors , Ċ, , see for instance Merodio [19] and Anssari-Benam et al. [18], two of which can be eliminated by virtue of the incompressibility property [14].As pointed out by Coco and Saccomandi [20], it seems unpractical to consider such a general theory as it leads to lengthy mathematical expressions.
Although it is not absolutely needed here, we will present a much simpler theory for practical use, namely the transversely isotropic viscous model by Spencer [21].This constitutive model is based on an expression of the dissipation potential in terms of the tensors  and  = ⊗ only, where  = ∕‖‖ is a unit vector defining the current direction of anisotropy, in the deformed configuration (as opposed to the initial direction of anisotropy  =  −1 ).Thus, the symmetric tensor  = ∕ 4 is idempotent with unit trace.With this assumption, the dissipation potential  v depends on only ten scalar invariants of , , six of which can be eliminated due to the properties of  and incompressibility.
Next, a Kelvin-Voigt model can be derived, where the viscous stress  v is chosen to be linear with respect to .With this further simplification, the dissipation potential is now function of  v ∶=  v ( 2 ,  4 ,  5 ) only, with the invariants defined as Here, we have eliminated those invariants that are cubic in the strain rate, because they lead to non-Newtonian viscosity theories.
We can now deduce a reduced form for the stress tensors.As  e depends only on the invariants (6), and  v depends only on the invariants (9), we can calculate the stress tensors (3) and ( 7) by using the chain rule, where   = 2  e ∕  and   = 2  v ∕  .The isotropic theory is recovered for  4 ,  5 ,  4 ,  5 equal to zero.In general, the elastic moduli   (in Pa) are functions of all the selected invariants   in Eq. ( 6).Thus, we observe that the tensor  =  4  can be replaced with  in the expression of the elastic stress (10) 1 , provided that the coefficients  4 ,  5 are given a suitable redefinition.
Because we assumed that the stress is linear in , the viscosities   (in Pa s) are no longer arbitrary functions of the invariants   of Eq. ( 9).In fact,  2 and  5 must be constants, whereas  4 can be linear in .Thus, we set  4 = 2 4  4 with  4 constant.With these choices, Eq. ( 10) 2 recovers exactly the transversely isotropic Newtonian stress proposed by Spencer [21], Eq. (34) therein. 1  According to the principles of thermodynamics, the dissipation per unit volume of material must remain nonnegative to ensure physically admissible behaviour [22].
Here, the dissipation reads tr( v ) = tr . By substituting equation ( 9) into the above, and using that  v is linear in , it can be shown that tr( v ) is a homogeneous function of degree two with respect to Ċ, and that the above expression equals 2 v .Given that the tensor 2 is positive semi-definite, we have that the invariants  2 ,  5 are non-negative.Therefore, requiring that the viscosities  2 ,  4 ,  5 to be non-negative is sufficient to ensure thermodynamic consistency.

General theory
In what follows, we linearise the equations of motion with respect to a pre-deformed state, with the goal of studying the propagation of infinitesimal shear waves in a material with initial stress.
Let us consider a static pre-deformation  ↦ x with the deformation gradient F =  x∕, see notations in Berjamin and Destrade [23].In contrast with the updated direction of anisotropy m = F n, the initial direction of anisotropy n =  is unaffected by the pre-deformation.
The initial deformation creates a static initial stress T = − p + T e , where T e is deduced from the constitutive law (10).To describe how elastic waves couple to this stress, we consider a small displacement ũ =  − x, then follow the 'small-on-large' approach [14].For this purpose, we introduce the increment of the deformation gradient F =  F such that  = F + F , where  = ∇ ũ is the incremental displacement gradient (here, the gradient stands for ∕ x).
The equations of motion are then linearised with respect to ũ, in the absence of body forces.Since the deformation is homogeneous, we arrive at the linearised balance of momentum equation [14]  ü = ∇ ⋅ σ, σ = − p + σe + σv , (12) in the pre-deformed configuration, where σ is the incremental Cauchy stress.The incremental displacement ũ is divergence-free due to the incompressibility constraint.As we consider a homogenous pre-deformation we have that T e and ℂ do not depend on space, we then compute the divergence of σe to get σe , =    , , where indices after the comma indicate spatial differentiation.We note in passing that the coefficients  , vanish because of incompressibility, and that the last two indices of  , can be interchanged because the order of differentiation can be interchanged.The above expressions are equivalent to that found in the supplementary material of Zhang et al. [11] up to conventions and notations used therein.
Similarly, we compute the incremental viscous stress by linearisation of (7), using the absence of strain rates in the pre-deformed configuration ( ̄Ċ = ).This way, we arrive at an expression of the viscous stress increment of the form where  = [  ] is a stiffness tensor and  = [  ] is the viscosity tensor.
For the transversely isotropic Kelvin-Voigt model described earlier, detailed expressions of the stress increments could be derived directly from Eq. ( 10) by linearisation, see for instance Destrade et al. [14] for the methodology.In particular, we note that the coefficients of  are equal to zero in this special case.

Harmonic plane waves
We apply this theory to the case of harmonic plane wave solutions of the form where  > 0 is the angular frequency and  is a unit vector along the direction of propagation.Similar notations are introduced for the incremental pressure and the incremental stresses.If the wavenumber  is complex-valued, then the complex velocity  = ∕ is a complex number too.
In general, the complex velocity differs from the phase velocity   = ∕ℜ().More precisely, computation of the real and imaginary parts of 1∕ yields [15] where  = −ℑ() is the attenuation coefficient.If   and  are known at a given frequency , then the complex velocity  = (1∕) −1 can be retrieved based on the above formulas.In the absence of attenuation ( = 0), the complex velocity is equal to the phase velocity.The constraint of incompressiblity ∇ ⋅ ũ = 0 leads to the orthogonality property  ⋅ û = 0, whereas the linearised equation of motion becomes where defines the acoustic tensor  and the dynamic tensor .On scalar multiplication of ( 18) by , we obtain the harmonic amplitude of the incremental pressure, p.Then, substitution in ( 18) yields the eigenvalue problem whose solutions provide the dynamic modulus  2 .For non-trivial solutions, the polarisation vector û can be assumed unitary, without loss of generality.Then, the dynamic modulus can be retrieved directly from ( 18)-( 19) by means of a scalar product with û: Noting that the second term of ℂ from ( 14) is formally similar to  and  , see (15), we write where  = [  ] is defined by The minor symmetries of  follow from the symmetry of .
In general, the equality of mixed partials entails the major symmetry   =   for  provided that  has this symmetry too, a requirement that is not needed here.

The angled shear wave identity
We now move on to the derivation of the angled shear wave identity.For this purpose, we assume that the material is subjected to a static triaxial pre-deformation which is aligned with coordinate axes so that whose determinant equals one, in agreement with the incompressibility property.This way, the unit vectors   ,   ,   in the reference configuration are equal to their counterpart   ,   ,   in the pre-deformed configuration.A picture is provided in Figure 1.Following similar steps shown in Zhang et al. [11] for the purely elastic case, we see that the angled shear wave identity only holds when the direction of anisotropy  is aligned with one of the principal directions of the predeformation   ,   , or   .This way, the vector m = λ   is aligned with  =   (no summation over repeated indices), where we have used the notation   =   for  = 1, etc. Furthermore, the tensor M = λ2    ⊗   is diagonal.Thus, the pre-stress T = diag[ σ1 , σ2 , σ3 ] is determined from the principal stresses σ = − p + σe  , where the coefficients σe  are the diagonal entries of T e , see Eq. ( 10) 1 , and there is no viscous stress as the pre-deformation is static.Along with these observations, it is worth mentioning that the tensors F , T and M are mutually commuting (i.e., their eigenspaces coincide), an alignment condition introduced in the supplementary material of Zhang et al. [11].
We now choose a direction of propagation  = cos    + sin    for the small amplitude shear wave (16).From the orthogonality property  ⋅ û = 0, we deduce that nontrivial wave solutions are either polarised along û =   (i.e., orthogonally to the propagation plane), or along û = − sin    +cos    , which belongs to the propagation plane.Similarly to earlier works, we now restrict the study to the latter case (in-plane polarisation), but a similar methodology could be used for other wave polarisations.
In a similar fashion to the elastic case [11], explicit calculation of the coefficients   shows that β vanishes.In fact, the coefficients of  involve tensor derivatives of the form  2  • ∕  evaluated at equilibrium ( = C, Ċ = ), with  ,  equal to  or Ċ.In the most general case, the scalar functions  • defined in ( 5)-( 8) depend on fifteen invariants of the form tr(  Ċ   ) with suitable integer exponents [19].Evaluation of the above partial derivatives by means of the chain rule then leads to β = 0, due to the fact that F and  are diagonal.
Finally, by taking β = 0, and rewriting in terms of the total stress components σ = − p+ σe  , equation (25) becomes whose coefficients do not depend on .Hence, for any frequency , we recover the angled shear wave identity (1) by subtraction of the dynamic moduli  2 and  ′2 measured along two angles  and  ′ = ∕2 ± , respectively (see Figure 1).This exact formula provides direct access to the stress difference σ1 − σ3 based on measurements of the phase velocity and wave attenuation from which  is deduced.
Illustration.Typically, in most ultrasonic experiments, it is simpler to measure the wave speed than it is to measure the attenuation.Here we illustrate what type of error would be expected when only using the wave speed in a viscoelastic solid to predict the stress when using the identity (1).
Using Eq. ( 17), we can write the squared phase velocity in the form where we have used a series expansions for small dispersion .Therefore, if only the wave speed   was measured and used instead of  in (1), then the predicted stress would be where  ′  ,  ′ are deduced from ( 29)-(30) with  replaced by  ′ .In contrast with the exact formula (1), the above expression includes an error of the order of ( 2 ) 2 .
This type of error might explain the mismatch found by Zhang et al. [11], according to whom "the stress is underestimated when the viscoelasticity comes into play".In the low frequency range of Figure S7B therein, the estimated stress decreases with increasing values of the frequency, which is coherent with Eq. (31) and ,  ′ ≥ 0. The above formula might also explain the frequency variations in Figure 6 of Delory et al. [24].

Concluding remarks
We have proved that the angled shear wave identity can be generalised, and still holds, for transversely isotropic soft solids with viscosity when replacing the phase velocity with the complex wave velocity.Therefore, if both the wavespeed and attenuation can be measured, then the identity (1) provides direct access to the initial stress, without extensive prior knowledge of the mechanical behaviour.However, for the identity to hold, the principal directions of the stress need to be aligned with the principal directions of the fibres or texture anisotropy.
The method we used to prove the identity in this paper involved a pre-deformation to create the initially stressed state.However we note that this proof could have been done for an arbitrary initial stress, without reference to how the stress was created, by following the approach shown in Mukherjee et al. [13] and Gower et al. [25], see also Ogden [26] for details.
Measuring the complex wave velocity is equivalent to measuring both the wave speed and wave attenuation.Such a measurement could be done with an imaging method, as is used for shear wave elastography [11].If the attenuation can not be measured, then using only the wave speed in place of the complex wave speed can lead to an approximate prediction of the stress.In the specific example considered here, the corresponding error is a quadratic order of the dispersion parameter.
By following similar steps, we believe the identity (1) could be generalised to more complex constitutive models, including compressible ones [5,13].A first step would be to apply the present method to orthotropic media with two orthogonal preferred directions.In a second step, frequencydependent responses that go beyond differential models could be investigated, such as models of fractional viscous stress or Maxwell rheologies based on an integral formulation.In the Appendix A, we show that the result still holds for some isotropic theories of this kind [23,27].
], and   denotes the Kronecker delta.In this formula, the overbar indicates that the tensor derivatives are evaluated in the pre-deformed equilibrium state where  = F .