This post is largely inspired by a preprint of Capoferri and Vassiliev.
Capoferri and Vassiliev considered a Lorentzian analogue of elasticity theory: given $(M,g)$ a semi-Riemannian manifold, and let $\phi: M\to M$ be a diffeomorphism, one can define the analogue of the Cauchy strain tensor from classical, geometric formulation of elasticity theory (where in the static case concerns $g$ being a Riemannian metric) by the formula (given in abstract index notation) \begin{equation} S^a{}_b = g^{ac} (\phi^\star g)_{cb} - \delta^a{}_b. \end{equation} Here $\phi$ is considered as the body deformation, and the strain compares the pull-back metric (a measure of the intrinsic geometry of the deformed body) with the original metric. Following ideas of elasticity theory, after noting that $S$ can be regarded as a linear operator on the tangent spaces $T_xM$, they decide to study Lagrangian field theories where the Lagrangian density takes the form \begin{equation} L ~\mathrm{dvol}_g \end{equation} where the function $L$ is restricted to be a function of the algebraic invariants of the operator $S$. They further restrict to only considering volume-preserving diffeomorphisms $\phi$; though the actual physical motivation for this escapes me.
Doing computations with $\phi$ is, let's say, tricky at best. Even if we assume that $\phi$ is a small perturbation of the identity mapping, the computations are still difficult: unless one has a global chart that covers $M$, parametrizing $\phi$ is not simple. One way to do so, that Capoferri and Vassiliev alluded to, is to treat small perturbations of the identity as described by the inverse exponential map: that is, instead of considering $\phi$, one can consider the field $A(x) = \exp_x^{-1} \phi(x)$. Provided $\phi$ is sufficiently close to the identity and the geometry of $M$ is sufficiently well-controlled, we can assume at every $x$, the value $\phi(x)$ lies within the injectivity radius and $A(x)$ is well-defined as a section of the tangent bundle of $M$. (The reason to consider $A(x)$ is that this casts the data in a linear space, which makes analysis easier.)
The formula for $S$ however is very complicated when it comes to the description in $A$. While the authors did not explicitly describe it, based on geometry, we know that fixing $x\in M$ and $v\in T_x M$, the pushforward $\phi_\star v$ has a geometric meaning in terms of $A(x)$: if we let $\gamma: [0,1]\to M$ be the geodesic $t\mapsto \exp_x tA(x)$, so that $\gamma(0) = x$ and $\gamma(1) = \phi(x)$, we know that $\phi_\star v \in T_{\gamma(1)} M$ is the value $V(1)$, where $V$ is the Jacobi field along $\gamma$ with $V(0) = v$ and $V'(0) = \nabla_v A(x)$. I do not easily see a way where $S$ can be made much simpler for the analysis.
One can consider an alternative formulation of their theory, however.
Instead of insisting on studying macroscopic deformations, one can study their infinitesimal analogue, even nonlinearly. In this setting the descriptions localizes much better and formulas seem better within reach. Let me summarize this alternative here. First, we will take seriously the vector field $A$, instead of the mapping $\phi$. We think of the vector field $A$ as describing an infinitesimal deformation of the body. The analogue of the Cauchy strain will be the deformation tensor for $A$: \begin{equation} S^a{}_b = g^{ac} (\mathscr{L}_A g)_{cb} \end{equation} here we use $\mathscr{L}$ to denote Lie differentiation. One can expand the deformation tensor \begin{equation} S^a{}_b = g^{ac} ( \nabla_c A_b + \nabla_b A_c) \end{equation} in terms of the symmetrized metric-covariant derivative of the vector field $A$. One can further require that the theory is incompressible, so that the $A$ represents a volume-preserving infinitesimal deformation. As is well-known this means that $A$ must be divergence free \begin{equation} \nabla_c A^c = 0 \end{equation} and so in this context this means that the first algebraic invariant of $S$ automatically vanishes. (Remember that the algebraic invariants are the coefficients1 of the characteristic polynomial of a linear operator, and the first invariant corresponds to the trace.) A classic result of linear algebra says that the $j$-th algebraic invariant is proportional to the contraction \begin{equation} S^{[a_1}{}_{a_1} S^{a_2}{}_{a_2} \cdots S^{a_j]}{}_{a_j} \end{equation} where the bracket notation $[ a_1 a_2 \ldots a_j]$ in the upper indices indicates implicit signed permutations. We see that the second matrix invariant then becomes \begin{equation} \sigma_2 \propto (\sigma_1)^2 - (\nabla_c A_b + \nabla_b A_c)(\nabla^c A^b + \nabla^b A^c). \end{equation} Notice that the first term vanishes by the incompressibility condition. Similarly, we can compute, just for illustration, the third matrix invariant; simplifying with the vanishing trace we get \begin{equation} \sigma_3 \propto (\nabla^a A_b + \nabla_b A^a)(\nabla^b A_c + \nabla_c A^b)(\nabla^c A_a + \nabla_a A^c). \end{equation} Since we are mainly interested in the physically setting with 4 space-time dimensions, we complete this with looking at the final case of the fourth matrix invariant, which can be simplified to (using incompressibility) \begin{equation} \sigma_4 \propto -2 \mathrm{tr} (S^4) + (\mathrm{tr} S^2)^2 \end{equation} where $S^4$ is the linear operator formed by the four-fold composition of $S$.
Based on the above computation, we see that any Lagrangian density $L$ that depends only on the algebraic invariants can also be written as a function of $s_j$, $j = 2,3,4$, where $s_j = \mathrm{tr} S^j$. This allows us to analyze the hyperbolicity of such a theory. The Euler-Lagrange equation are given in divergence form easy enough: \begin{equation} 0 = \nabla_a \left( 2 \partial_2 L S^{ab} + 3 \partial_3 L (S^2)^{ab} + 4 \partial_4 L (S^3)^{ab} \right) \end{equation} where $\partial_j L = \partial L / \partial s_j$. A full analysis of the principle symbol of these equations is beyond the scope here today. But we can note that the principal part of the equation is given as some combination of \[ \nabla_a S^{ac} = \Box A^c + \mathrm{Ric}^{ac} A_a \] and terms of the form \[ S^{ab} \nabla_b S_{cd} \] and hence for sufficiently small $A$, the equations satisfied are quasilinear wave/Klein-Gordon equations (depending on the Ricci curvature of the background) and the dynamics are hyperbolic, provided $\partial_2 L \neq 0$.
- The computations of such invariants also featured in one of my papers concerning Lagrangian field theory. ^