\section*{Introduction}

Attempts to formulate \emph{Relativistic Elasticity} in the \emph{General Relativity} framework go back to 1916 with the pioneering work of Nordström~\cite{Nor1916}, in Dutch. Since then, several authors have first aimed at proposing constitutive equations for Relativistic fluids~\cite{Tau1954,Lic1955,CQ1972,MTW1973} and, then, at modeling Relativistic continuous media, most often at the astrophysics scale~\cite{Sou1958,Syn1959,deW1962,Ray1963,Sou1964,Ben1965,Old1970,Lam1989,KM1992,KM1997,BS2003,EBT2006,Wer2006,GHE2011,Bro2021}, for instance for the modeling of the solid crust of neutron stars, but also at a local scale~\cite{Mau1978a,Mau1978b,Mau1978c,PR2013,PR2013,PRA2015,NWP2022}, for mechanical engineering applications.

This has led Lichnerowicz to define \emph{pure matter}~\cite{Lic1955}, synonymous of dust, and Souriau to define \emph{perfect matter}~\cite{Sou1958,Sou1960,Sou1964}, as a continuous medium which can be described independently from electromagnetic phenomena. In the present work, we follow Souriau and model perfect matter with the \emph{Gauge Theory} mindset~\cite{Ble1981}. More precisely, we focus on Relativistic hyperelastic continuous media. We do not consider the coupling with electromagnetism, nor with temperature.

The work of Souriau (1958, in French), seems to have been unnoticed by the scientific community. It is prior to the works of Synge (1959), of DeWitt (1962) and of Rayner (1963) (all three criticized in the later papers by Bennoun~\cite{Ben1965} and Carter and Quintana~\cite{CQ1972}). As we shall see, Souriau did in fact formulate the correct framework to describe Relativistic Hyperelasticity, first in his long 1958 paper~\cite{Sou1958}, then in his 1964 book~\cite{Sou1964} (in French still). The modern geometric picture of the General Relativity framework for elastic media is, up to details that we shall discuss on the go, derived in~\cite{Sou1958,Sou1964}, and later in~\cite{CQ1972,KM1992,KM1997,BS2003}.

We stick to the chronology introduced by Souriau of the mathematical concepts, the idea being as often to make the least hypotheses as possible. This is why the introduction of a \emph{foliation by spacelike hypersurfaces} ---which is not assumed \emph{a priori}--- is only addressed in \autoref{sec:Universe-foliation}, and why the problems of the \emph{definition of time} and of the formulation of Relativistic Hyperelasticity in a \emph{spacetime} is addressed only in \autoref{sec:hyperelasticity-spacetime}.
We finally apply Künzle's methodology \cite{Kuen1976}, to mathematically recover Classical Galilean Hyperelasticity with gravity, as the Newton--Cartan infinite light speed limit of the described General Relativity formulation. Our calculations generalize the ones for relativistic fluids to relativistic solids.

We write this paper mainly with the mechanics ---not the astrophysics--- point of view. We seek for the adequacy with the geometric formulation on the body $\mathcal{B}$ of three-dimensional Hyperelasticity (the so-called \emph{intrinsic Lagrangian formulation}, developed by Noll~\cite{Nol1972,Nol1978} and Rougée~\cite{Rou1991a,Rou2006}, see also~\cite{KD2021,KD2021a}). The notations are chosen to be compatible with both the ones used classically in Continuum Mechanics of solid materials~\cite{LC1985} and the ones considered in~\cite{KD2021,KD2021a}.

\subsection*{Outline}

The article is organized as follows. In \autoref{sec:matter-fields}, we introduce the basic concepts of \emph{matter field} $\Psi$, of \emph{body World tube} $\mathcal{W}$, of \emph{mass measure $\mu$}, and of \emph{current of matter} $\mathbf{P}$, which are the starting point of the theory of Relativistic Continuum Mechanics. The normalization of the later allows to define the \emph{rest mass density} $\rho_{r}$ and the unit timelike vector $\mathbf{U}$. The section ends with the definition of the \emph{conformation} $\mathbf{H}$, the cornerstone of Souriau's general covariant formulation of Relativistic Hyperelasticity. Matter conservation is formulated in \autoref{sec:mass-conservation}. Definitions of relativistic strain tensors are provided in \autoref{sec:Conformation-strains}. Relativistic Hyperelasticity is formulated in \autoref{sec:Lagrangian-formulation}, in which the Lagrangian formulation of General Relativity is recalled and applied to this specific constitutive modeling. The \emph{stress-energy tensor} is introduced in \autoref{sec:stress-energy-tensor} and both four-dimensional and three-dimensional relativistic \emph{stress tensors} are defined. It is only in \autoref{sec:Universe-foliation} that a \emph{spacetime structure} is considered, allowing to better connect the preceding general geometric framework with Classical Continuum Mechanics, and to define the \emph{generalized Lorentz factor} $\gamma$ (\autoref{sec:Matter-field-spacetime}). This factor accounts for the distortion between the unit vector $\mathbf{U}$ (the matter) and the unit normal $\mathbf{N}$ to the spacelike hypersurfaces $\Omega_{t}$ (the observer), and allows for the geometric definition of the \emph{relativistic mass density} $\rho$. A first formulation of Relativistic Hyperelasticity in a static spacetime, including the generalization of the \emph{Cauchy stress tensor,} is derived in \autoref{sec:hyperelasticity-spacetime}. This framework is detailed in \autoref{sec:hyperelasticity-Schwarzschild} for the particular case of the \emph{Schwarzschild metric}. The Galilean (Newton--Cartan) infinite light speed limit of the theory is discussed in \autoref{sec:Gallilean}.

\subsection*{Notations}

Given a linear operator $L\colon E \to F$ between two finite dimensional vector spaces, we denote by $L^{\star} \colon F^{\star} \to E^{\star}$, $\beta \mapsto \beta \circ L$, its \emph{transpose}. If moreover, the vector space $E$ is equipped with an inner product $q_{E}$, and $F$, with an inner product $q_{F}$, we can define its \emph{adjoint}, defined implicitly by the relation $\overline{L} \colon F \to E$, $\langle Lv,w \rangle_{F} = \langle v,\overline{L}w \rangle_{E}$ for all $v \in E$ and $w \in F$. The relation between $L^{\star}$ and $\overline{L}$ is thus written as $\overline{L} = q_{E}^{-1} L^{\star} q_{F}$.
We denote by $\Sym^{k}E$ the set of totally symmetric tensors of order $k$ on $E$ and by $\Lambda^{k}E$ the set of alternate tensors of order $k$ on $E$.

Now, let $\mM$ be a differential manifold of dimension $n$, we denote by $\Omega^{k}(\mM)$, the set of \emph{differential $k$-forms} on $\mM$, that is smooth sections of the vector bundle $\Lambda^{k}T^{\star}\mM$. The contraction $i_{X} \alpha$, of components $(i_{X} \alpha)_{ \mu_{1} \dotsc \mu_{k-1}}=X^{\nu} \alpha_{\nu \mu_{1} \dotsc \mu_{k-1}}$, denotes the \emph{interior product} of a vector field $X\in \mathrm{Vect}(\mM)$ with a $k$-form $\alpha \in \Omega^{k}(\mM)$,
\begin{equation*}
  i_{X} \alpha:=\alpha(X, \cdot, \dotsc, \cdot) \in  \Omega^{k-1}(\mM).
\end{equation*}
If moreover, $\mM$ is endowed with a Riemannian or pseudo-Riemannian metric $g$ (and $\mM$ is orientable), we will denote by $\vol_{g}\in \Omega^{n}(\mM)$ the (pseudo-)Riemannian \emph{volume form} associated with $g$.

Given a 1-form $\alpha \in \Omega^{1}(\mM)$, the notation $\alpha^{\sharp}:=g^{-1}\alpha\in \mathrm{Vect}(\mM)$ stands for $\alpha^{\mu}=g^{\mu\nu} \alpha_{\nu}$. Conversely, given a vector field $X\in \mathrm{Vect}(\mM)$, $X^{\flat}:=gX\in \Omega^{1}(\mM)$ stands for $X_{\mu}=g_{\mu\nu} X^{\nu}$. When local coordinates are involved on a 4-dimensional manifold, the Greek subscripts or superscripts $\mu, \nu, \rho\dots$ range from 0 to 3, while the roman ones $i, j, k\dots$, or $I, J, K, \dots$ range from 1 to 3.

The light speed will be denoted by $c$ and we refer to the Galilean three-dimensional Continuum Mechanics of solids~\cite{TN1965,LC1985,MH1994} simply as \emph{Classical Continuum Mechanics}.

\section{Matter field, current of matter and conformation}
\label{sec:matter-fields}

The Universe is assumed to be a four-dimensional orientable manifold $\mM$, endowed with an hyperbolic metric $g$, of signature
$(-,+,+,+)$. Its pseudo-Riemannian volume form is denoted by
\begin{equation*}
  \vol_{g}\in \Omega^{4}(\mM).
\end{equation*}
In the present work we limit our study to a (non electromagnetic) continuous particles assembly, the so-called \emph{perfect matter}~\cite{Lic1955}. Its modeling adopted by Souriau in~\cite{Sou1958,Sou1960,Sou1964} is inspired by Gauge theory~\cite{Ble1981}, where \emph{matter fields} are described by sections of an \emph{associated bundle}, \textit{i.e.}, some vector bundle constructed using a linear representation of the \emph{structural group} of the considered Gauge theory on some given vector space. The specificity and relative simplicity of the present description of \emph{perfect matter} is, however, that we assume this linear representation, and thus the vector bundle, to be trivial. More precisely, we let $V$ be a three-dimensional vector space (taken as $\RR^{3}$ in~\cite{Sou1958}). A \emph{perfect matter field} (called the particles labelling in~\cite{Sou1958}, and the \emph{projection}, noted $\mathcal P$, in~\cite{CQ1972}) is then a smooth vector valued function
\begin{equation*}
  \Psi: \mM \to V.
\end{equation*}

\begin{rem}
  The notation $\Psi$ for the matter field is on purpose chosen similar to the one for the wave function in Quantum Mechanics.
\end{rem}

Matter is then described by the set of all the material points constitutive of the continuous medium under study in the Universe
(for example a mechanical structure). Their labels constitute a set $\mathcal{B}\subset V$, assumed to be (in general) a three-dimensional compact orientable manifold with boundary and called the \emph{body}. It is further assumed that $\Psi$ is a \emph{submersion} on $\mathcal{W}=\Psi^{-1} (\mathcal{B})$: the linear tangent map $T\Psi: T\mathcal{W} \to TV$ is of rank 3 at each point of $\mathcal{W}$. Thus, $\mathcal{W}$ is fibered by the particles World lines $\Psi^{-1}(\mathbf{X})$, $\mathbf{X}\in \mathcal{B}$, and is called for this reason the \emph{body's World tube}.

\begin{figure}[ht]\label{fig:fiberedtimelines}
  \centering
  \includegraphics[width=14cm]{Fig1.pdf}
  \caption{The World tube $\mathcal{W}=\Psi^{-1} (\mathcal{B})$ fibered by the particles World lines $\Psi^{-1}(\mathbf{X})$.}
\end{figure}

The body $\mathcal{B}$ is endowed with a volume form $\mu\in \Omega^{3}(\mathcal{B})$, the \emph{mass measure}, which carries the information about the distribution of matter present in $\mathcal{W}$~\cite{KM1992}. This interpretation is connected with the three-dimensional Classical Continuum Mechanics theory, in which the abstract manifold $\mathcal{B}$, equipped with the mass measure $\mu$, is in fact the \emph{body} introduced by Truesdell and Noll~\cite{TN1965,Nol1972,Nol1974,Nol1978}.

As we seek for a full consistency with the geometric framework of Classical Continuum Mechanics~\cite{Nol1978,Rou1991a,Rou2006,KD2021}, we have to emphasize a slight difference with previous works in astrophysics concerning the choice of the volume form on $\mathcal{B}$.
In~\cite{Sou1958}, $V$ is equipped with the canonical 3-form $\vol_{\bq}=\mathrm{d} X^{1} \wedge \mathrm{d} X^{2} \wedge \mathrm{d} X^{3}$ on $\RR^{3}$. In~\cite{CQ1972,BS2003,GHE2011}, the body $\mathcal{B}$ is equipped with a volume form $\Omega\in \Omega^{3}(\mathcal{B})$ which represents the number density of conserved idealized particles (meant to be identified with the baryon number density in~\cite{CQ1972}). The three of them are, of course, proportional to each other on $\mathcal{B}$. As pointed out by Carter and Quintana, the Relativistic Hyperelasticity theory does not depend on the particular choice of a volume form on $\mathcal{B}$~\cite{CQ1972}. Our choice, here, of a volume form $\mu$, interpreted as a ``mass measure'' allows us to recover the mass densities encountered in Classical Continuum Mechanics, and to assimilate the integral
\begin{equation*}
  m=\int_{\mathcal{B}} \mu,
\end{equation*}
as the total mass of the continuous medium/mechanical structure under study.

\begin{rem}\label{rem:F-CCM}
  It is worth mentioning that one takes here a point of view reverse to the one of Classical Continuum Mechanics of solids~\cite{TN1965,LC1985,Rou1991a,MH1994,Rou2006}, in which a configuration is an embedding $p\colon\mathcal{B} \to \espace$ of the body $\mathcal{B}$ into the three-dimensional space $\espace$, endowed with the Euclidean metric $\bq$. In the present formalism, the main concept is a mapping $\Psi\colon \mM \to V$ from the Universe $\mM$ to the space of labels $V$. A key difference is that, in Classical Continuum Mechanics, $p$ and its tangent map, the so-called \emph{deformation gradient}
  \begin{equation*}
    \mathbf{F}=Tp \colon T\mathcal{B} \to T\espace
  \end{equation*}
  are invertible, whereas here, the matter field $\Psi$ and its tangent map $T\Psi$ are not.
\end{rem}

The pullback by $\Psi$ of the mass measure $\mu$ on the body $\mathcal{B}$
\begin{equation*}
  \omega:=\Psi^{*} \mu=(\mu\circ \Psi)(T\Psi \cdot, T\Psi \cdot, T\Psi \cdot)
\end{equation*}
is a 3-form defined on the four-dimensional World tube $\mathcal{W}=\Psi^{-1}(\mathcal{B})$. Since $T_{m}\Psi$ is assumed to be of rank $3$ at each point of $\mathcal{W}$, there exists a \emph{nowhere vanishing} vector field $\mathbf{P}$ on $\mathcal{W}$, such that
\begin{equation}\label{eq:def-P}
  \omega=i_{\mathbf{P}} \vol_{g},
\end{equation}
where $i_{\mathbf{P}}$ is the interior product (or contraction) of $\omega$ by $\mathbf{P}$. This vector field $\mathbf{P}$ is the \emph{current of matter} (it was called \textit{vecteur courant de matière} in~\cite{Sou1958}).

\begin{rem}
  In 3D Classical Continuum Mechanics, the pushforward of the mass measure $\mu$ by the embedding $p\colon \mathcal{B} \to \espace$ \cite{KD2021}, when expressed using the 3D volume form $\vol_{q}$, is represented by a scalar density $\rho$ (indeed, $p_{*} \mu=\rho \, \vol_{q}$). In 4D, the pullback of the mass measure $\mu$ by the matter field $\Psi$, when expressed using the 4D volume form $\vol_{g}$, is represented by the quadrivector~$\mathbf{P}$ (indeed, $\Psi^{*} \mu=i_{\mathbf{P}} \vol_{g}$).
\end{rem}

At each point $m\in \mathcal{W}$, the tangent vector $\mathbf{P}(m)$ spans the one-dimensional subspace $\ker T_{m}\Psi$.
Indeed, the equality
\begin{equation*}
  \mu_{\Psi(m)}(T_{m}\Psi.\mathbf{P}, T_{m}\Psi.\xi_{1}, T_{m}\Psi.\xi_{2}) = \vol_{g_{m}}(\mathbf{P}(m),\mathbf{P}(m),\xi_{1},\xi_{2})=0,
  \quad
  \forall  \xi_{1},\xi_{2} \in T_{m}\mathcal{W},
\end{equation*}
implies
\begin{equation}\label{eq:TPsi-P-null}
  T_{m}\Psi.\mathbf{P}(m)=0,
\end{equation}
since $T_{m}\Psi$ is surjective. To describe perfect matter, Souriau assumes furthermore that $\mathbf{P}$ is timelike, \emph{i.e.} that
\begin{equation*}
  \norm{\mathbf{P}}_{g}^{2}=g(\mathbf{P}, \mathbf{P})<0
\end{equation*}
on the World tube $\mathcal{W}$ (we refer to~\cite{Sou1964} for the other cases, light for instance). Observe that $\mathbf{P}$ defines a time orientation on $\mathcal{W}$.

It will be proved as essential to define a unit timelike vector field $\mathbf{U}$ collinear to $\mathbf{P}$, and to write
\begin{equation}\label{eq:def-U}
  \mathbf{P}=\rho_{r} \mathbf{U},
  \quad
  \text{with}
  \quad
  \norm{\mathbf{U}}^{2}_{g}=-1.
\end{equation}
The function
\begin{equation}\label{eq:def-rho}
  \rho_{r} :=\sqrt{-\norm{\mathbf{P}}_{g}^{2}} \, ,
\end{equation}
defined on the World tube $\mathcal{W}$, is then interpreted as the \emph{rest mass density} \cite{Sou1958,Sou1964,KM1992}.

\begin{rem}\label{eq:cP}
  In Special Relativity, the vector field $c\mathbf{P}$ is the \emph{four-momentum quadrivector}.
\end{rem}

We will finish this section by defining the \emph{conformation}, a fundamental concept introduced by Souriau in 1958. It is the cornerstone of the formulation of Relativistic Hyperelasticity at large scale, in particular, for the modeling of neutron stars with a solid crust. In recent works, it is sometimes referred to as \emph{strain}, but since \emph{this term} has a slightly different meaning in Classical Continuum Mechanics, we prefer to keep the initial name given by Souriau. The conformation is defined as the vector-valued function~\cite{Sou1958}
\begin{equation}\label{eq:conformation}
  \mathbf{H}: \mM \to \Sym^{2}V, \qquad m \mapsto \mathbf{H}(m):= (T_{m}\Psi)\,g_{m}^{-1} (T_{m}\Psi)^{\star},
\end{equation}
where $\Sym^{2}V$ is the six-dimensional vector space of symmetric contravariant second-order tensors on $V$. In simpler words, $\mathbf{H}$ is a function from the four-dimensional manifold $\mM$ to the vector space of $3 \times 3$ matrices. The hypothesis we made that $\Psi$ is a submersion on $\mathcal{W}$, together with the hypothesis that $\ker T\Psi$ is generated by the timelike vector field $\mathbf{U}$ implies that $\mathbf{H}(m)$ is positive definite for all $m \in \mathcal{W}$.

\begin{rem}\label{rem:Psim1}
  Since the mapping $\Psi: \mathcal{W} \to \mathcal{B}$ is not invertible, the conformation cannot be considered, \textit{stricto sensu}, as the pushforward of $g^{-1}$, which is not defined. It is thus not, strictly speaking, a co-metric on $\mathcal{B}$, but a vector-valued function of $m \in \mathcal{W}$ with value a symmetric second-order contravariant tensor in $V$. Note that, if we forget that $\mathcal{B}$ is a domain in the vector space $V$ but consider that it is a manifold, then $\mathbf{H}$ is interpreted as a tensor field along $\Psi$ with values in $\Sym^{2} T\mathcal{B}$, in other words, it is a section of the pullback bundle $\Psi^{*}(\Sym^{2} T\mathcal{B})$.
\end{rem}

\section{Conservation of matter}
\label{sec:mass-conservation}

Since the exterior derivative of the mass measure $\mu$ on the 3-dimensional manifold $\mathcal{B}$ vanishes, $\mathrm{d} \mu=0$, we get the following conservation law.

\begin{lem}[Souriau, 1958] \label{lem:divP}
  We have the following conservation law on the World tube $\mathcal{W}$
  \begin{equation*}
    \dive^{g} \mathbf{P}=0.
  \end{equation*}
\end{lem}

\begin{proof}
  Let $\Lie_{\mathbf{P}}$ be the Lie derivative with respect to $\mathbf{P}$. Then,
  \begin{equation*}
    \Lie_{\mathbf{P}} \vol_{g} = (\dive^{g} \mathbf{P})\, \vol_{g},
  \end{equation*}
  but, using Cartan magic formula,
  \begin{equation*}
    \Lie_{\mathbf{P}} \vol_{g} = \mathrm{d} \, i_{\mathbf{P}} \vol_{g}=\mathrm{d}\Psi^{*} \mu=\Psi^{*} \mathrm{d}\mu=0.
  \end{equation*}
\end{proof}

\begin{rem}
  In Special Relativity, the equation
  \begin{equation*}
    \dive^{g} \mathbf{P} = \dive^{g} (\rho_{r} \mathbf{U})=0
  \end{equation*}
  recasts as the usual continuity equation of Classical Fluid Dynamics \cite{Eck1940,Sou1958,Sou1964}, and is interpreted as the Relativistic mass conservation. It will be shown, furthermore, in \autoref{sec:Gallilean}, that $ \dive^{g} c \mathbf{P}=0$ converges towards the classical continuity equation under subsequent hypothesis.
\end{rem}

If the body $\mathcal{B}$ is endowed with a Riemannian metric ${\bm{\gamma}}_{0}$, the rest mass density $\rho_{r}$ can be related to the conformation $\mathbf{H}$, as demonstrated by Souriau in~\cite{Sou1958}, where he chose ${\bm{\gamma}}_{0}:=\bq$, the canonical Euclidean metric on $V=\RR^{3}$ (see also~\cite{KM1992}). The notation ${\bm{\gamma}}_{0}$ is chosen for consistency with the intrinsic geometric framework of three-dimensional Hyperelasticity~\cite{Nol1978,Rou1991a,Rou2006,KD2021}, and ${\bm{\gamma}}_{0}$ is not necessarily equal to $\bq$, as discussed by several authors~\cite{BS2003,GHE2011} (see \autoref{sec:reference-metric} for a discussion about different choices for ${\bm{\gamma}}_{0}$).

\begin{lem}[Souriau, 1958] \label{lem:PCM_GR}
  Let ${\bm{\gamma}}_{0}$ be a fixed Riemannian metric on the body $\mathcal{B}$. Then, the rest mass density $\rho_{r}$ can be expressed as
  \begin{equation}\label{eq:PCM_GR}
    \rho_{r} = (\rho_{{\bm{\gamma}}_{0}} \circ \Psi) \sqrt{\det\left[ \mathbf{H}({\bm{\gamma}}_{0} \circ \Psi)\right]},
  \end{equation}
  where $\Psi$ is the matter field, $\mathbf{H}$ is the conformation, and
  \begin{equation*}
    \rho_{{\bm{\gamma}}_{0}} = \frac{\mu}{\vol_{{\bm{\gamma}}_{0}}}.
  \end{equation*}
\end{lem}

\begin{rem}\label{rem:PCM_GR}
  The function $\rho_{{\bm{\gamma}}_{0}}$ is defined on the body $\mathcal{B}$, and interpreted as the mass density with respect to the Riemannian volume form $\vol_{{\bm{\gamma}}_{0}}$. It is very important to note, for subsequent applications, that $\rho_{{\bm{\gamma}}_{0}}$ is independent of the metric $g$ on the Universe $\mM$. Moreover, one can check that the right hand-side of~\eqref{eq:PCM_GR} does not depend on ${\bm{\gamma}}_{0}$, as expected. Indeed if we substitute ${\bm{\gamma}}_{1}$ to ${\bm{\gamma}}_{0}$ in~\eqref{eq:PCM_GR}, one has $\rho_{{\bm{\gamma}}_{1}} = \sqrt{\det({\bm{\gamma}}_{1}^{-1}{\bm{\gamma}}_{0})}\,\rho_{{\bm{\gamma}}_{0}}$, whereas $\det(\mathbf{H}{\bm{\gamma}}_{1}) = \det({\bm{\gamma}}_{0}^{-1}{\bm{\gamma}}_{1}) \det( \mathbf{H}{\bm{\gamma}}_{0})$, and thus
  \begin{equation*}
    (\rho_{{\bm{\gamma}}_{0}} \circ \Psi) \sqrt{\det\left[ \mathbf{H}({\bm{\gamma}}_{0} \circ \Psi)\right]} = (\rho_{{\bm{\gamma}}_{1}} \circ \Psi) \sqrt{\det\left[ \mathbf{H}({\bm{\gamma}}_{1} \circ \Psi)\right]}.
  \end{equation*}
\end{rem}

\begin{proof}
  Note first that there exists a function $\rho_{{\bm{\gamma}}_{0}}$ (a mass density) defined on $\mathcal{B}$ such that
  \begin{equation*}
    \mu = \rho_{{\bm{\gamma}}_{0}} \vol_{{\bm{\gamma}}_{0}}.
  \end{equation*}
  Thus, by~\eqref{eq:def-P}--\eqref{eq:def-U}, we get
  \begin{equation*}
    \omega = \rho_{r} \, i_{\mathbf{U}}\vol_{g} = \Psi^{*}\mu = (\rho_{{\bm{\gamma}}_{0}}\circ\Psi) \Psi^{*}\vol_{{\bm{\gamma}}_{0}}.
  \end{equation*}
  Let $m \in \mathcal{W}$ and $(\bm{e}_{0} = \mathbf{U}(m),\bm{e}_{i})$ be a direct orthonormal basis of $T_{m}\mM$. Then, we have
  \begin{equation*}
    \vol_{g_{m}}(\mathbf{U}(m),\bm{e}_{1},\bm{e}_{2},\bm{e}_{3}) = 1
  \end{equation*}
  and
  \begin{equation*}
    \rho_{r}(m) = \rho_{r}(m)\vol_{g_{m}}(\mathbf{U}(m),\bm{e}_{1},\bm{e}_{2},\bm{e}_{3}) = \rho_{{\bm{\gamma}}_{0}}(\Psi(m)) \, \vol_{{\bm{\gamma}}_{0}(\Psi(m))}(T_{m}\Psi \bm{e}_{1}, T_{m}\Psi \bm{e}_{2},T_{m}\Psi \bm{e}_{3}).
  \end{equation*}
  Observe now that the vector space $T_{\Psi(m)}\mathcal{B}$ is endowed with two Euclidean structures; the first one, defined by ${\bm{\gamma}}_{0}(\Psi(m))$ and the second one, defined by $\mathbf{H}(m)^{-1}$. Besides, the restriction of $T_{m}\Psi$ to the three-dimensional subspace $\mathbf{U}(m)^{\bot}$ (the orthogonal complement of $\mathbf{U}(m)$ in $T_{m}\mM$) is a linear isomorphism and
  \begin{equation*}
    T_{m}\Psi\colon (\mathbf{U}(m)^{\bot},g_{m}) \to (T_{\Psi(m)}\mathcal{B},\mathbf{H}(m)^{-1})
  \end{equation*}
  is an isometry, by the very definition of the conformation $\mathbf{H}$. Hence, $(T_{m}\Psi \bm{e}_{1}, T_{m}\Psi \bm{e}_{2},T_{m}\Psi \bm{e}_{3})$ is a direct orthonormal basis of the Euclidean space $(T_{\Psi(m)}\mathcal{B},\mathbf{H}(m)^{-1})$ and thus
  \begin{equation*}
    \vol_{\mathbf{H}(m)^{-1}}(T_{m}\Psi \bm{e}_{1}, T_{m}\Psi \bm{e}_{2},T_{m}\Psi \bm{e}_{3}) = 1.
  \end{equation*}
  We have therefore
  \begin{align*}
    \vol_{{\bm{\gamma}}_{0}(\Psi(m))}(T_{m}\Psi \bm{e}_{1}, T_{m}\Psi \bm{e}_{2},T_{m}\Psi \bm{e}_{3})
     & = \sqrt{\det\left[\mathbf{H}(m){\bm{\gamma}}_{0}(\Psi(m))\right]}\, \vol_{\mathbf{H}(m)^{-1}}(T_{m}\Psi \bm{e}_{1}, T_{m}\Psi \bm{e}_{2},T_{m}\Psi \bm{e}_{3}) \\
     & =  \sqrt{\det\left[\mathbf{H}(m){\bm{\gamma}}_{0}(\Psi(m))\right]},
  \end{align*}
  and thus
  \begin{equation*}
    \rho_{r}(m) = \rho_{{\bm{\gamma}}_{0}}(\Psi(m)) \, \sqrt{\det\left[\mathbf{H}(m){\bm{\gamma}}_{0}(\Psi(m))\right]}.
  \end{equation*}
\end{proof}

\begin{rem}\label{rem:PCM_Omega0}
  In Classical Continuum Mechanics, the body is often identified with a reference configuration $\Omega_{0}$ embedded in $\RR^{3}$ and endowed with the euclidean metric ${\bm{\gamma}}_{0}=q$. Two mass densities, $\rho_{0}$ on $\Omega_{0}$ and $\rho$ on the deformed configuration $\Omega$ (also embedded in $\RR^{3}$), are usually defined. A classical expression of mass balance is formulated on $\Omega_{0}$ as
  \begin{equation}\label{eq:PCM_Omega0}
    \rho_{0} = (\rho \circ \phi) \sqrt{\det(\bq^{-1}\mathbf{C})},
  \end{equation}
  where $\phi \colon \Omega_{0}\to \Omega$ is the \emph{deformation}, $\mathbf{C}:=\phi^{*}\bq$ is the \emph{right Cauchy--Green tensor} (defined on $\Omega_{0}$ as the pullback by the deformation $\phi$ of the Euclidean metric $\bq$). The formal comparison of~\eqref{eq:PCM_GR}, recast as
  \begin{equation*}
    \rho_{{\bm{\gamma}}_{0}}\circ \Psi=\rho_{r}  \sqrt{\det\left[ ({\bm{\gamma}}_{0}^{-1} \circ \Psi)\mathbf{H}^{-1}\right]},
  \end{equation*}
  with ~\eqref{eq:PCM_Omega0}, shows that~\eqref{eq:PCM_GR} can be interpreted as a Relativistic generalization of the mass conservation law for Galilean deformable solids. It also shows that the inverse of the (contravariant) conformation $\mathbf{H}$ plays the role of the (covariant) right Cauchy--Green tensor $\mathbf{C}$.
\end{rem}

\section{Conformation and strains}
\label{sec:Conformation-strains}

The existence of the unit timelike vector field $\mathbf{U}$ on the World tube $\mathcal{W}$ allows to perform the related orthogonal decompositions of the metric $g$ and co-metric $g^{-1}$ (see~\autoref{sec:Orth-Decomp-2nd}),
\begin{equation}\label{eq:def-h}
  g = h-\mathbf{U}^\flat \otimes \mathbf{U}^{\flat},
  \qquad
  g^{-1} = h^{\sharp}-\mathbf{U} \otimes \mathbf{U},
  \qquad \text{on $\mathcal{W}$},
\end{equation}
where the tensor fields $h$ (noted $E$ in \cite{KM1997}) and $h^{\sharp}=g^{-1} h g^{-1}$, the spatial part of $g$ and $g^{-1}$ respectively, are uniquely defined by the conditions
\begin{equation}\label{eq:hU}
  h\mathbf{U} = 0, \quad \text{and} \quad h = g \quad \text{on} \quad  \mathbf{U}^{\perp},
\end{equation}
where $\mathbf{U}^{\perp}$ is the three-dimensional (necessarily spacelike) orthogonal subbundle to $\mathbf{U}$. Both $h$ and $h^{\sharp}$ have signature $(0,+,+,+)$. These orthogonal decompositions are highlighted at the beginning of most works on Relativistic Fluids or Solids \cite{Eck1940,Lic1955,CQ1972,KM1997}. We point out, however, that Souriau did not need to perform them to derive the general covariant formulation of Relativistic Hyperelasticity \cite{Sou1958,Sou1964}. There are two reasons for it. First, the four-dimensional symmetric second-order tensors $h$ and $h^{\sharp}$ are strongly related to the conformation
\begin{equation*}
  \mathbf{H} = (T\Psi)\,g^{-1} (T\Psi)^{\star},
\end{equation*}
by lemma~\ref{lem:h-H}. Secondly, $h$ and $h^{\sharp}$ do not appear naturally in the derivation of a general covariant formulation of Relativistic Hyperelasticity, contrary to the conformation $\mathbf{H}$ (see theorem~\ref{thm:LPsiH}).

\begin{lem}\label{lem:h-H}
  On the World tube $\mathcal{W}$, we have
  \begin{equation}\label{eq:Hhsharp}
    \mathbf{H} = (T\Psi)\,h^{\sharp} (T\Psi)^{\star},
    \quad \text{and}\quad
    h=(T\Psi)^{\star}\, \mathbf{H}^{-1} T\Psi.
  \end{equation}
  where $h=g+\mathbf{U}^\flat \otimes \mathbf{U}^{\flat}$ and $h^{\sharp}=g^{-1} h g^{-1}$.
\end{lem}

\begin{proof}
  First, since $T\Psi \mathbf{U}=0$, the conformation, when restricted to the World tube $\mathcal{W}$, recasts as
  \begin{equation*}
    \mathbf{H} = (T\Psi)\left(h^{\sharp} - \mathbf{U} \otimes \mathbf{U} \right)(T\Psi)^{\star}=(T\Psi)\,h^{\sharp} (T\Psi)^{\star}.
  \end{equation*}
  Then, to prove the second equality, remark that the statement is pointwise. Therefore, we can use an orthonormal basis $(\bm{e}_{\mu})$ of $T_{m}\mM$ with $\bm{e}_{0} = \mathbf{U}(m)$. In this basis, $h_{m}$ is represented by the $4 \times 4$ matrix
  \begin{equation*}
    \begin{pmatrix}
      0 & 0     \\
      0 & I_{3}
    \end{pmatrix},
  \end{equation*}
  where $I_{3}$ is the $3 \times 3$ identity matrix. Now, respectively to this basis and the canonical basis of $\RR^{3}$, the linear map $T_{m}\Psi\colon T_{m}\mM \to \RR^{3}$ is represented by the matrix
  \begin{equation*}
    \begin{pmatrix}
      0 & M \\
    \end{pmatrix},
  \end{equation*}
  where $M$ is a $3 \times 3$ invertible matrix, and its transpose $(T_{m}\Psi)^{\star}$ by the matrix
  \begin{equation*}
    \begin{pmatrix}
      0         \\
      M^{\star} \\
    \end{pmatrix}.
  \end{equation*}
  Thus, we have
  \begin{equation*}
    \mathbf{H}(m) = (T_{m}\Psi)\,g_{m}^{-1} (T_{m}\Psi)^{\star} =
    \begin{pmatrix}
      0 & M \\
    \end{pmatrix}
    \begin{pmatrix}
      -1 & 0     \\
      0  & I_{3}
    \end{pmatrix}
    \begin{pmatrix}
      0         \\
      M^{\star} \\
    \end{pmatrix}
    = MM^{\star},
  \end{equation*}
  and
  \begin{equation*}
    (T_{m}\Psi)^{\star}\, \mathbf{H}(m)^{-1} T_{m}\Psi =
    \begin{pmatrix}
      0         \\
      M^{\star} \\
    \end{pmatrix}
    M^{-\star}M^{-1}
    \begin{pmatrix}
      0 & M \\
    \end{pmatrix}
    =
    \begin{pmatrix}
      0     \\
      I_{3} \\
    \end{pmatrix}
    \begin{pmatrix}
      0 & I_{3} \\
    \end{pmatrix}
    = h_{m}.
  \end{equation*}
\end{proof}

The definition of a strain in (hyper)elasticity is usually obtained by comparing two metrics. If the body $\mathcal{B}$ is endowed with a fixed Riemannian metric ${\bm{\gamma}}_{0}$, it can be used to define a strain tensor in Relativistic Hyperelasticity. A first possibility~\cite{CQ1972,Mau1978b} is to introduce the pullback by $\Psi$ of ${\bm{\gamma}}_{0}$,  given by
\begin{equation}\label{eq:h0}
  h_{0} := \Psi^{*}{\bm{\gamma}}_{0}=(T\Psi)^{\star} ({\bm{\gamma}}_{0}\circ \Psi) T\Psi,
\end{equation}
and called a \emph{frozen metric} in~\cite{KM1992,KM1997} (these authors note it $h$ rather than $h_{0}$).
It is defined on the World tube $\mathcal{W}$ and is of signature $(0, +, +, +)$. Conversely, given a quadratic form $h_{0}$ on $\mathcal{W}$ with signature $(0, +,+,+)$, the question of when it can be realized as the pullback by $\Psi$ of a fixed Riemannian metric ${\bm{\gamma}}_{0}$ on the body, has been investigated by Kijowski and Magli~(see \autoref{sec:reference-metric}).

A possible generalization of the \emph{Euler-Almansi strain tensor}~\cite{CQ1972,Mau1978b} is then obtained as the four-dimensional symmetric covariant tensor field,
\begin{equation}\label{eq:strains-h}
  \mathbf{e} := \frac{1}{2} (h-h_{0}).
\end{equation}
Note that $\mathbf{e}=0$ for $h=h_{0}$ and that $\mathbf{e}$ is degenerate since $\mathbf{e} \mathbf{U}=0$.

\begin{rem}
  As observed by Carter and Quintana \cite{CQ1972}, since the linear tangent map $T\Psi$ plays a role similar to that of the inverse of the tangent map $\mathbf{F}=Tp$ in Classical Continuum Mechanics (see remark~\ref{rem:F-CCM}), the frozen metric $h_{0}$ plays a role similar to that of the inverse, sometimes called the \emph{finger deformation tensor}, of the \emph{left Cauchy--Green tensor} $\bb:=\mathbf{F}{\bm{\gamma}}_{0}^{-1} \mathbf{F}^{\star}$.
\end{rem}

Other choices for strain tensors similar to the ones of Classical Continuum Mechanics can be made, for instance the following ones which are simpler and probably more relevant,
\begin{equation}\label{eq:strains-H}
  \mathfrak E:= \frac{1}{2} \left(\mathbf{H}^{-1}-\mathbf{H}_{0}^{-1}\right)
  \quad \text{or} \quad
  \widehat{\mathfrak E} := -\frac{1}{2}\log \big(\mathbf{H}\,\mathbf{H}_{0}^{-1}\big),
\end{equation}
where
\begin{equation} \label{eq:H0}
  \mathbf{H}_{0}:={\bm{\gamma}}_{0}^{-1} \circ \Psi.
\end{equation}
The first one generalizes the \emph{Green--Lagrange strain}, whereas the second one generalizes the \emph{logarithmic strain} introduced by Becker~\cite{Bec1893} and Hencky~\cite{Hen1928} (see~\cite{MMEN2018}). They both vanish when $\mathbf{H}^{-1}=\mathbf{H}_{0}^{-1}={\bm{\gamma}}_{0}\circ \Psi$. These strain tensors are three-dimensional second-order tensors. Like the conformation, they are not tensor fields on $\mathcal{B}$ but vector valued functions defined on the World tube $\mathcal{W}$ with values in $\Sym^{2}V$.

Note that $\mathbf{H}_{0}$ is related to $h_{0}$ by
\begin{equation*}
  h_{0}=(T\Psi)^{\star}\,\mathbf{H}_{0} \, (T\Psi)
\end{equation*}
and that, by lemma~\ref{lem:h-H} and definitions~\eqref{eq:h0}--\eqref{eq:strains-h}, $\mathfrak E$ is connected to $\mathbf{e}$, defined by~\eqref{eq:strains-h}, by
\begin{equation*}
  \mathbf{e}=(T\Psi)^{\star}\, \mathfrak E\, (T\Psi)
  \quad
  \text{on $\mathcal{W}$}.
\end{equation*}

\section{Lagrangian formulation}
\label{sec:Lagrangian-formulation}

In~\cite{Sou1958,Sou1964}, Souriau has proposed a clear and detailed formulation of Hyperelasticity in the framework of General Relativity. He called this formulation \emph{Variational Relativity} (which is the title of~\cite{Sou1958}). His approach consists in writing Lagrangians (\textit{i.e.} functionals depending on tensorial fields) and looking for critical points of them (Principle of Least, or Stationary, Action). This formulation is inspired by Gauge Theory~\cite{Ble1981}, which is the main framework of Fields Theory and Quantum Mechanics and can also be used to formulate General Relativity using variational principles (see Palatini's Method~\cite{Pal1919,FFR1982}).

The starting point is the \emph{Hilbert-Einstein functional}
\begin{equation}\label{eq:Hilbert-Einstein}
  \mathcal{H}(g) = \int (a R_{g} + b)\, \vol_{g},
\end{equation}
defined formally on the set of all Lorentzian metrics on the Universe $\mM$. Here, the two constants $a$ and $b$ are related to the Einstein constant $\kappa$ (depending on the Newton constant $G$) and the cosmological constant $\Lambda$ by
\begin{equation*}
  \kappa = \frac{8\pi G}{c^{4}} = \frac{1}{2a}, \qquad \Lambda = -\frac{b}{2a}.
\end{equation*}
As derived first by Hilbert~\cite{Hil1924}, the $L^{2}$-gradient of $\mathcal{H}$ (for Ebin's metric~\cite{Ebi1968}) is the symmetric second order covariant tensor field
\begin{equation}\label{eq:H-gradient}
  \grad \mathcal{H} = a\, \Ric_{g}-\frac{1}{2}(a\, R_{g}+b)g = \frac{1}{2\kappa}\left( \Ein_{g} + \Lambda g\right),
\end{equation}
where $\Ric_{g}$ is the Ricci tensor of the metric $g$, $R_{g} = \tr(g^{-1}\Ric_{g})$ is the scalar curvature, and $\Ein_{g}$ is the \emph{Einstein tensor}, defined by
\begin{equation}\label{eq:Einstein-tensor}
  \Ein_{g} := \Ric_{g}-\frac{1}{2}R_{g}\, g.
\end{equation}

The critical points of $\mathcal{H}$ are the solutions of Einstein's equation in the vacuum (with cosmological constant)
\begin{equation*}
  \Ein_{g} + \Lambda g = 0.
\end{equation*}

To introduce the effects of matter in this framework, a second functional $\mathcal{L}^{\text{matter}}(g, \Psi)$, depending on the metric $g$ and the matter field $\Psi$, is added to $\mathcal{H}$ to build a new Lagrangian
\begin{equation*}
  \mathcal{L}(g,\Psi) = \mathcal{H}(g) + \mathcal{L}^{\text{matter}}(g, \Psi).
\end{equation*}
Following Souriau~\cite{Sou1958,Sou1964}, for Relativistic continua, one assumes that the Lagrangian for perfect matter $\mathcal{L}^{\text{matter}}(g, \Psi)$ depends only on the $0$-jet of the metric $g$ and of the $1$-jet of the matter field $\Psi$. In other words, it takes the form
\begin{equation}\label{eq:LagrangianFirstJets}
  \mathcal{L}^{\text{matter}}(g, \Psi) = \int L_{0}\left(g_{\mu\nu}, \Psi^{I}, \frac{\partial \Psi^{I}}{\partial x^{\mu}} \right) \vol_{g},
\end{equation}
where
\begin{equation*}
  L_{0} \colon ({\bm{\gamma}},\bm{v},{\bm{\tau}}) \to L_{0}({\bm{\gamma}},\bm{v},{\bm{\tau}})
\end{equation*}
is a smooth scalar function which has for arguments a quadratic form ${\bm{\gamma}}$ (of signature $(-,+,+,+)$) on $\RR^{4}$, a vector $\bm{v} \in \RR^{3}$ and matrix ${\bm{\tau}}$ with 3 raws and 4 columns. The function $L_{0}$ is called the \emph{Lagrangian density} of the functional $\mathcal{L}^{\text{matter}}$, and its evaluation on the fields $(g,\Psi)$, that is $L_{0}\left(g_{\mu\nu}, \Psi^{I}, \frac{\partial \Psi^{I}}{\partial x^{\mu}} \right)$, will be denoted as $L_{0}(g, \Psi, T\Psi)$. Its evaluation at a point $m\in \mM$ is then noted $L_{0}(g_{m}, \Psi(m), T_{m}\Psi)$.

\begin{rem}
  The Lagrangian density $L_{0}$ is noted $p$ and called the \emph{presence} in~\cite{Sou1958}. It is noted $\epsilon$ and called the \emph{rest frame energy density} in~\cite{KM1992}. It is noted $\rho$ or $\sigma$ in \cite{BS2003}.
\end{rem}

In order to avoid unnecessary analytical difficulties and since, in practice, we do not require that Lagrangian densities are integrable over the whole manifold $\mM$, usually not compact, Lagrangian densities are integrated only over relatively
compact domains $U$ (and furthermore contained in a local chart). Therefore, we shall write
\begin{equation*}
  \mathcal{L}_{U}(g, \Psi) = \int_{U} L(g_{m}, \Psi(m), T_{m}\Psi)\, \vol_{g},
\end{equation*}
to emphasize the dependence on $U$. When we just want to express that a Lagrangian $\mathcal{L}$ is defined by the Lagrangian density
$L$, we simply write
\begin{equation*}
  \mathcal{L}(g, \Psi) = \int L(g_{m}, \Psi(m), T_{m}\Psi) \, \vol_{g},
\end{equation*}
omitting the domain of integration.

The main postulate of General Relativity is precisely that \emph{Physical laws must be independent of the choice of coordinates}. This principle is known as \emph{General Covariance}, or invariance by coordinates change, or invariance by (local) diffeomorphisms. Let us describe this principle in more precise terms and formulate its consequences. Let $\varphi \colon U \to \widetilde{U}$ be a diffeomorphism between two open sets $U$ and $\widetilde{U}$. Then, the Lagrangian $\mathcal{L}$ is invariant by $\varphi$ if
\begin{equation}\label{eq:invmL}
  \mathcal{L}_{U}(\varphi^{*} g, \varphi^{*}\Psi) = \mathcal{L}_{\widetilde{U}}(g, \Psi),
\end{equation}
for every Lorentzian metrics $g$ on $\mM$, and vector valued functions $\Psi:\mM\to V$. Here, the action of a (local) diffeomorphism $\varphi$ on these field variables is defined by
\begin{equation*}
  \varphi^{*} g = (T\varphi)^{\star}\,(g \circ \varphi) (T\varphi), \quad \text{and} \quad \varphi^{*} \Psi=\Psi \circ \varphi .
\end{equation*}
If the invariance~\eqref{eq:invmL} holds for every local diffeomorphism $\varphi \colon U \to \widetilde{U}$, then, $\mathcal{L}$ is said to be \emph{general covariant}.

\begin{rem}\label{rem:H-invariance}
  It is well-known that the Hilbert-Einstein functional $\mathcal{H}$ is general covariant. Indeed,
  \begin{equation*}
    \mathcal{H}_{U}(\varphi^{*} g) = \mathcal{H}_{\widetilde{U}}(g),
  \end{equation*}
  for every diffeomorphism $\varphi\colon U \to \widetilde{U}$, by virtue of the \emph{change of variables formula}
  \begin{equation*}
    \int_{U} \varphi^{*} \omega = \int_{\widetilde{U}} \omega,
  \end{equation*}
  and because
  \begin{equation*}
    (a R_{\varphi^{*} g} + b)\, \vol_{\varphi^{*} g} = \varphi^{*}[(a R_{g} + b)\, \vol_{g}].
  \end{equation*}
  As shown by Noether \cite{Noe1918,Kos2011}, a direct consequence of this invariance is the fundamental property~\cite{Ein1915,Wey1917}
  \begin{equation*}
    \dive^{g}(\Ein_{g} + \Lambda g) = \dive^{g} \Ein_{g} = 0.
  \end{equation*}
\end{rem}

\begin{lem}
  If the Lagrangian
  \begin{equation*}
    \mathcal{L}^{\text{matter}}(g, \Psi) = \int L_{0}(g_{m}, \Psi(m), T_{m}\Psi) \, \vol_{g}
  \end{equation*}
  is general covariant, then, its Lagrangian density satisfies
  \begin{equation}\label{eq:Lagrangian-density-covariance}
    L_{0}(\mathbf{A}^{\star} {\bm{\gamma}} \,\mathbf{A}, \bm{v}, {\bm{\tau}}\mathbf{A}) = L_{0}({\bm{\gamma}}, \bm{v}, {\bm{\tau}}),
    \quad
    \forall \mathbf{A} \in \GL(4).
  \end{equation}
\end{lem}

\begin{proof}
  Let $\varphi \colon U \to \widetilde{U}$ be a diffeomorphism between two open sets $U$ and $\widetilde{U}$ and set
  \begin{equation*}
    f(\mathbf{y}) := L_{0}(A(\mathbf{y})^{\star}g_{\mathbf{y}} A(\mathbf{y}),\Psi(\mathbf{y}), T_{\mathbf{y}}\Psi A(\mathbf{y}))
  \end{equation*}
  for $\mathbf{y} \in \widetilde{U}$, where $A(\mathbf{y}) = T_{\varphi^{-1}(\mathbf{y})}\varphi$. Then, $T_{m}\varphi = A(\varphi(m))$, for $m \in U$ and
  \begin{equation*}
    L_{0}((\varphi^{*}g)_{m}, (\varphi^{*}\Psi)(m), T_{m}(\varphi^{*}\Psi)) = f(\varphi(m)).
  \end{equation*}
  Therefore
  \begin{equation*}
    \mathcal{L}_{U}(\varphi^{*} g, \varphi^{*}\Psi) = \int_{U} f(\varphi(m)) \varphi^{*} \vol_{g} = \int_{\widetilde{U}} f(\mathbf{y}) \vol_{g},
  \end{equation*}
  by the change of variable formula, and the general covariance property leads to
  \begin{equation*}
    f(\mathbf{y}) = L_{0}(A(\mathbf{y})^{\star}g_{\mathbf{y}} A(\mathbf{y}),\Psi(\mathbf{y}), T_{\mathbf{y}}\Psi A(\mathbf{y})) = L_{0}(g_{\mathbf{y}},\Psi(\mathbf{y}), T_{\mathbf{y}}\Psi), \quad \forall \mathbf{y}.
  \end{equation*}
  Hence, the Lagrangian density is subject to the following invariance
  \begin{equation*}
    L_{0}(\mathbf{A}^{\star} {\bm{\gamma}} \mathbf{A}, \bm{v}, {\bm{\tau}} \mathbf{A}) = L_{0}({\bm{\gamma}}, \bm{v}, {\bm{\tau}}),
    \quad
    \forall \mathbf{A} \in \GL(4).
  \end{equation*}
\end{proof}

\begin{rem}
  Since the Lie derivative is the infinitesimal version of the pullback, meaning that
  \begin{equation*}
    \Lie_{X}\mathbf{T} := [\partial_{s}\varphi(s)^{*}\mathbf{T}]_{s=0}
  \end{equation*}
  for every tensor field $\mathbf{T}$ and every path of (local) diffeomorphisms $\varphi(s)$ with
  \begin{equation*}
    \varphi(0)=\mathrm{id}, \quad \text{and} \quad [\partial_{s}\varphi(s)]_{s=0}=X,
  \end{equation*}
  there is also an almost\footnote{indeed equivalent to covariance by diffeomorphisms isotopic to the identity.} equivalent \emph{infinitesimal formulation of general covariance}~\cite{Noe1918}, which is used by several authors (such as in \cite{Wer2006}). For instance, in the present case, the general covariance of the matter Lagrangian $\mathcal{L}^{\text{matter}}$
  \begin{equation*}
    \mathcal{L}^{\text{matter}}_{U}(\varphi^{*} g, \varphi^{*}\Psi) = \mathcal{L}^{\text{matter}}_{\widetilde{U}}(g, \Psi),
  \end{equation*}
  for every local diffeomorphism $\varphi \colon U \to \widetilde{U}$ leads to
  \begin{equation*}
    \frac{\delta\mathcal{L}^{\text{matter}}_{U}}{\delta g} . \Lie_{X}g + \frac{\delta\mathcal{L}^{\text{matter}}_{U}}{\delta \Psi}.\Lie_{X}\Psi = 0.
  \end{equation*}
  Therefore, its Lagrangian density must satisfy (see~\cite{Sop2008})
  \begin{equation*}
    \frac{\partial L_{0}}{\partial {\bm{\gamma}}}:(\Lie_{X}g)_{m} + \frac{\partial L_{0}}{\partial \bm{v}}\cdot (\Lie_{X}\Psi)(m) + \frac{\partial L_{0}}{\partial {\bm{\tau}}}:(T_{m}\Lie_{X}\Psi) = 0, \qquad \forall m.
  \end{equation*}
\end{rem}

The following result is essential for the formulation of Relativistic Hyperelasticity and exhibits the fundamental role played by the conformation. It must be compared to the fact that an elastic energy in Classical Continuum Mechanics, which is \emph{objective} (\textit{i.e.} satisfies the \emph{material frame indifference} principle~\cite{TN1965}) depends on the deformation $\varphi$ only through the right Cauchy--Green tensor $\mathbf{C} =\varphi^{*}q$.

\begin{thm}[Souriau (1958)]\label{thm:LPsiH}
  Suppose that the Lagrangian
  \begin{equation*}
    \mathcal{L}^{\text{matter}}(g, \Psi) = \int L_{0}(g_{m}, \Psi(m), T_{m}\Psi) \, \vol_{g}
  \end{equation*}
  is general covariant. Then, its Lagrangian density can be written as
  \begin{equation*}
    L_{0}(g, \Psi, T\Psi) = L(\Psi, \mathbf{H}),
  \end{equation*}
  for some function $L$, where $\mathbf{H}= (T\Psi)\,g^{-1} (T\Psi)^{\star}$ is the conformation.
\end{thm}

The proof provided below is simpler and shorter that the original one given by Souriau in~\cite{Sou1958}. The reason for it is that, in this paper, we consider only perfect matter, in which case the conformation $\mathbf{H}$ is positive definite at each point $m$ of the World tube.
This is not an hypothesis which is made in~\cite{Sou1958}.

\begin{proof}
  Consider a smooth Lagrangian density $L_{0}({\bm{\gamma}}, \bm{v}, {\bm{\tau}})$, where ${\bm{\gamma}}$ is a quadratic form of signature $(-,+,+,+)$ on $\RR^{4}$, $\bm{v} \in \RR^{3}$ and ${\bm{\tau}}\in \mathcal{L}(\RR^{4},\RR^{3})$ satisfies ${\bm{\tau}}\overline{{\bm{\tau}}} >0$. Suppose moreover that this Lagrangian density satisfies the following covariance property
  \begin{equation*}
    L_{0}(\mathbf{A}^{\star} {\bm{\gamma}} \mathbf{A}, \bm{v}, {\bm{\tau}} \mathbf{A}) = L_{0}({\bm{\gamma}}, \bm{v}, {\bm{\tau}}),  \qquad \forall \mathbf{A}\in \GL(4).
  \end{equation*}

  First, we can find $\mathbf{A}\in \GL(4)$ such that $\mathbf{A}^{\star} {\bm{\gamma}} \mathbf{A}=\eta$, where
  \begin{equation*}
    \eta =
    \begin{pmatrix}
      -1 & 0 \\
      0  & q
    \end{pmatrix}
  \end{equation*}
  is the canonical Lorentz inner product. Hence we get
  \begin{equation*}
    L_{0}({\bm{\gamma}}, \bm{v}, {\bm{\tau}}) = L_{0}(\eta, \bm{v}, {\bm{\tau}}_{1}), \quad \text{with} \quad {\bm{\tau}}_{1} = {\bm{\tau}} \mathbf{A}
    \quad \text{and} \quad \mathbf{A}^{\star} {\bm{\gamma}} \mathbf{A} = \eta.
  \end{equation*}

  Now, we introduce the following change of variables ${\bm{\tau}}_{1} \mapsto (R, \mathbf{H})$, where
  \begin{equation*}
    \mathbf{H} = {\bm{\tau}}_{1}\eta^{-1} {\bm{\tau}}_{1}^{\star} = {\bm{\tau}}{\bm{\gamma}}^{-1}{\bm{\tau}}^{\star}, \qquad R = V^{-1}{\bm{\tau}}_{1},
  \end{equation*}
  and $V$ is the positive square root of the positive definite symmetric operator on $(\RR^{3},q)$
  \begin{equation*}
    {\bm{\tau}}_{1}\overline{{\bm{\tau}}_{1}} = \mathbf{H} q =  {\bm{\tau}}_{1}\eta^{-1}{\bm{\tau}}_{1}^{\star}q
  \end{equation*}
  with $q$, the canonical Euclidean metric on $\RR^{3}$.

  We can check that $R\overline{R} = R \eta^{-1}R^{\star}q = I_{3}$ is a condition which defines a submanifold of the vector space of linear mappings $\mathcal{L}(\RR^{4},\RR^{3})$, and that ${\bm{\tau}}_{1} \mapsto (R, \mathbf{H})$ is a diffeomorphism from the open set
  \begin{equation*}
    \set{{\bm{\tau}}_{1} \in \mathcal{L}(\RR^{4},\RR^{3});\; {\bm{\tau}}_{1}\overline{{\bm{\tau}}_{1}} > 0}
  \end{equation*}
  onto the manifold
  \begin{equation*}
    \set{R \in \mathcal{L}(\RR^{4},\RR^{3});\; R\overline{R} = I_{3}} \times \set{\mathbf{H} \in \Sym^{2}(\RR^{3});\; \mathbf{H} >0}.
  \end{equation*}
  Hence, we can find a smooth function $L_{1}(\bm{v}, R, \mathbf{H})$, such that
  \begin{equation*}
    L_{0}(\eta, \bm{v}, {\bm{\tau}}_{1}) = L_{1}(\bm{v}, R, \mathbf{H}),
  \end{equation*}
  with the property that
  \begin{equation*}
    L_{1}(\bm{v}, R, \mathbf{H}) = L_{0}(\eta, \bm{v}, {\bm{\tau}}_{1}) = L_{0}(\eta, \bm{v}, {\bm{\tau}}_{1}Q) = L_{1}(\bm{v}, RQ, \mathbf{H}),
  \end{equation*}
  for every Lorentz transformation $Q$. Next, we can find a Lorentz transformation $Q$ such that $RQ = R_{0}$ with
  \begin{equation*}
    R_{0} =
    \begin{pmatrix}
      0 & I_{3}
    \end{pmatrix},
  \end{equation*}
  because $R\overline{R} = R_{0}\overline{R_{0}} = I_{3}$. Therefore, we get finally
  \begin{equation*}
    L_{0}({\bm{\gamma}}, \bm{v}, {\bm{\tau}}) = L_{0}(\eta, \bm{v}, {\bm{\tau}}_{1}) =  L_{1}(\bm{v}, R, \mathbf{H}) = L_{1}(\bm{v}, R_{0}, \mathbf{H}),
  \end{equation*}
  and $L_{1}(\bm{v}, R_{0}, \mathbf{H})$ is a function $L(\bm{v}, \mathbf{H})$, which depends only on $\bm{v}$ and $\mathbf{H} = {\bm{\tau}}_{1}\eta^{-1} {\bm{\tau}}_{1}^{\star} = {\bm{\tau}}{\bm{\gamma}}^{-1}{\bm{\tau}}^{\star}$.
\end{proof}

The following splitting of the Lagrangian density has been introduced by Souriau~\cite{Sou1958,Sou1964} and DeWitt~\cite{deW1962}:
\begin{equation}\label{eq:LSouriau}
  L(\Psi, \mathbf{H}) = \rho_{r} c^{2}+E(\Psi, \mathbf{H})=\rho_{r} c^{2}+\rho_{r} e(\Psi, \mathbf{H}),
\end{equation}
where $\rho_{r}$ is the rest mass density, expressed as
\begin{equation*}
  \rho_{r}=\rho_{{\bm{\gamma}}_{0}}(\Psi) \sqrt{\det\left[ \mathbf{H}({\bm{\gamma}}_{0} \circ \Psi)\right]},
\end{equation*}
by lemma~\ref{lem:PCM_GR}, provided a fixed metric ${\bm{\gamma}}_{0}$ has been given on the body $\mathcal{B}$ and $\rho_{{\bm{\gamma}}_{0}}=\mu/\vol_{{\bm{\gamma}}_{0}}$. The contribution $\rho_{r} c^{2}$ alone ($E=0$) allows for the modeling of perfect (non electromagnetic) dust. The function $E$ (resp. $e$) is the \emph{internal energy density} (resp. the \emph{specific  internal energy}). It is representative of perfect fluids when its dependency on $\mathbf{H}$ is introduced only through the determinant $\det\left[ \mathbf{H}({\bm{\gamma}}_{0} \circ \Psi)\right]$. The additional dependency on $\Psi$ and $\mathbf{H}$ through the energy density $E$ is more generally representative of Relativistic hyperelastic solids.

We conclude this section by emphasizing that the present formulation of Relativistic Hyperelasticity does not require the definition of a time function (which is indeed not a necessity in astrophysics) and the associated assumption of a foliation of the World tube $\mathcal{W}$ by spacelike hypersurfaces. All we need is to endow the body $\mathcal{B}$ with a fixed metric ${\bm{\gamma}}_{0}$ as in~\cite{Sou1958,Sou1964,CQ1972}.

\section{The stress--energy tensor}
\label{sec:stress-energy-tensor}

The \emph{stress-energy tensor}, also called \emph{energy-momentum tensor} can be considered as a four-dimensional generalization of the stress tensor in Classical three-dimensional Continuum Mechanics. In General Relativity, it is the source of the curvature of the metric $g$ of the Universe. It is usually defined as the variational derivative of a Lagrangian with respect to the metric $g$ \cite{Hil1915,Hil1924,Noe1918,Kos2011,BGRS2016} and, for this reason, it is thus a symmetric contravariant second-order tensor field (or a \emph{tensor distribution} defined on symmetric second-order covariant tensor fields, in more general situations \cite{Sou1997a}).

In the present case, the Euler-Lagrange stationary equation $\delta \mathcal{L}=0$ for the Lagrangian
\begin{equation*}
  \mathcal{L}(g,\Psi)= \mathcal{H}(g) + \mathcal{L}^{\text{matter}}(g, \Psi),
\end{equation*}
leads in particular to the equation
\begin{equation*}
  \frac{\delta \mathcal{H}}{\delta g} + \frac{\delta \mathcal{L}^{\text{matter}}}{\delta g} = 0,
\end{equation*}
when only variations of the metric $g$ are considered. It recasts as the Einstein field equation
\begin{equation}\label{eq:Einstein}
  \Ein_{g}^{\sharp} + \Lambda g^{-1} = \kappa\mathbf{T},
\end{equation}
if $\Ein_{g}^{\sharp}=g^{-1}\Ein_{g} g^{-1}$ is the contravariant form of Einstein's tensor~\eqref{eq:Einstein-tensor}, and
\begin{equation*}
  \mathbf{T}:=-2\frac{\delta \mathcal{L}^{\text{matter}}}{\delta g}.
\end{equation*}
is the \emph{stress-energy tensor} (the source term in Einstein's equation), which is a symmetric contravariant
second-order tensor field on the Universe $\mM$.

\begin{rem}
  Because $\dive^{g}(\Ein_{g}^{\sharp} + \Lambda g^{-1})=0$ (see remark~\ref{rem:H-invariance}), the stress-energy tensor $\mathbf{T}$ satisfies the conservation law
  \begin{equation*}
    \dive^{g}\mathbf{T}=0.
  \end{equation*}
  As observed by Einstein himself~\cite{Ein1988}, ``$\dive \mathbf{T} = 0$, that's mechanics''. Indeed, this equation generalizes in 4D (and non flat Universe) the three-dimensional equilibrium equations of Classical Continuum Mechanics. When a spacetime structure is adopted, the Cauchy stress tensor is related to the spacelike components of $\mathbf{T}$ (see \autoref{sec:hyperelasticity-spacetime}).
\end{rem}

The following result provides a general expression for the stress-energy tensor of $\mathbf{T}$ in the case of Relativistic Hyperelasticity (see also \cite{KM1992}).

\begin{thm}[Souriau, 1958]\label{thm:relativistic-hyperelasticity}
  Consider the general covariant matter Lagrangian
  \begin{equation*}
    \mathcal{L}^{\text{matter}}(g, \Psi) = \int L\, \vol_{g},
    \qquad
    L=\rho_{r} c^{2}+E,
  \end{equation*}
  with
  \begin{equation*}
    \rho_{r} = \rho_{{\bm{\gamma}}_{0}}(\Psi) \sqrt{\det\left[ \mathbf{H}({\bm{\gamma}}_{0} \circ \Psi)\right]} \quad \text{and} \quad E = E(\Psi, \mathbf{H}),
  \end{equation*}
  and where ${\bm{\gamma}}_{0}$ is a fixed metric on the body $\mathcal{B}$. Then, its stress-energy tensor has the following expression
  \begin{equation}\label{eq:def-S}
    \mathbf{T} = -2\frac{\delta \mathcal{L}^{\text{matter}}}{\delta g} = \rho_{r} c^{2} \mathbf{U} \otimes \mathbf{U} - \mathbf{S},
  \end{equation}
  where
  \begin{equation*}
    \mathbf{S} := E\,g^{-1} -2 g^{-1} (T\Psi)^{\star}\frac{\partial E}{\partial \mathbf{H}} (T\Psi) g^{-1}.
  \end{equation*}
  is the (four-dimensional) \emph{relativistic stress tensor} on $\mathcal{W}$. Moreover, we have
  \begin{equation*}
    \mathbf{S} \cdot \mathbf{U}^{\flat}=E\, \mathbf{U},
    \quad
    \text{and}
    \quad
    \mathbf{T} \cdot \mathbf{U}^{\flat}=-L \mathbf{U}.
  \end{equation*}
\end{thm}

\begin{rem}[Bennoun, 1965]\label{rem:T}
  Since $ \mathbf{S} \cdot \mathbf{U}^{\flat}\neq 0$, the decomposition~\eqref{eq:def-S} is not an orthogonal decomposition relative to $\mathbf{U}$ (see \autoref{sec:Orth-Decomp-2nd}). Writing $E = \rho_{r} e$, with $e = e(\Psi, \mathbf{H})$, the specific internal energy, the stress-energy tensor naturally recasts, using its orthogonal decomposition relative to $\mathbf{U}$, as
  \begin{equation}\label{eq:def-Sigma}
    \mathbf{T} = L\,\mathbf{U} \otimes \mathbf{U} -{\bm{\Sigma}},
  \end{equation}
  where its spatial part
  \begin{equation}\label{eq:dedH}
    {\bm{\Sigma}}:=- 2\rho_{r}\, g^{-1} (T\Psi)^{\star}\frac{\partial  e}{\partial \mathbf{H}} (T\Psi) g^{-1},
  \end{equation}
  is such that
  \begin{equation*}
    {\bm{\Sigma}}=\mathbf{S} + E\, \mathbf{U} \otimes \mathbf{U}
    \quad \text{and} \quad {\bm{\Sigma}} \cdot \mathbf{U}^{\flat}=0,
  \end{equation*}
  can also be interpreted as a (four-dimensional) relativistic stress tensor.
\end{rem}

\begin{proof}
  Consider the variation $\delta_{g}\mathcal{L}^{\text{matter}}$ of the Lagrangian $\mathcal{L}^{\text{matter}}$ with respect to the metric $g$. Then, we have
  \begin{equation*}
    \delta_{g}\mathcal{L}^{\text{matter}} = \int (\delta_{g}L)\, \vol_{g} + L\, \delta_{g}\vol_{g},
  \end{equation*}
  with
  \begin{equation*}
    \delta_{g}L = \tr \left(\frac{\partial L}{\partial \mathbf{H}} \delta_{g}\mathbf{H}\right) \quad \text{and} \quad \delta_{g} \vol_{g} = \frac{1}{2} \tr(g^{-1}\delta g) \vol_{g}.
  \end{equation*}
  But
  \begin{equation*}
    \delta_{g}\mathbf{H} = - (T\Psi) g^{-1}\delta g\, g^{-1}(T\Psi)^{\star},
  \end{equation*}
  and hence
  \begin{equation*}
    \delta_{g}L  = - \tr \left(\frac{\partial L}{\partial \mathbf{H}} (T\Psi) g^{-1}\delta g\, g^{-1}(T\Psi)^{\star}\right)
    = - \tr \left(g^{-1}(T\Psi)^{\star}\frac{\partial L}{\partial \mathbf{H}} (T\Psi) g^{-1}\delta g \right).
  \end{equation*}
  We get thus
  \begin{equation*}
    \delta_{g}\mathcal{L}^{\text{matter}} = - \int \tr \left[ \left( g^{-1}(T\Psi)^{\star}\frac{\partial L}{\partial \mathbf{H}} (T\Psi) g^{-1} - \frac{1}{2} L g^{-1} \right) \delta g \right]\vol_{g},
  \end{equation*}
  and therefore
  \begin{equation*}
    \mathbf{T} = 2 g^{-1} (T\Psi)^{\star}\frac{\partial L}{\partial \mathbf{H}} T\Psi g^{-1}-  L\,g^{-1}.
  \end{equation*}
  Now, we have
  \begin{equation*}
    \frac{\partial \rho_{r}}{\partial \mathbf{H}} = \frac{1}{2} \rho_{r} \; \mathbf{H}^{-1},
  \end{equation*}
  and thus
  \begin{align*}
    \mathbf{T} & = \rho_{r}c^{2}  \, g^{-1} (T\Psi)^{\star}\mathbf{H}^{-1} (T\Psi) \, g^{-1}
    +2 g^{-1} (T\Psi)^{\star}\frac{\partial E}{\partial \mathbf{H}} (T\Psi) g^{-1}
    -  \rho_{r} c^{2} \,g^{-1}-  E\,g^{-1}
    \\
        & = \rho_{r} c^{2}\left[g^{-1}h g^{-1}-   \,g^{-1}\right]
    +2 g^{-1} (T\Psi)^{\star}\frac{\partial E}{\partial \mathbf{H}} (T\Psi) g^{-1}
    -  E\,g^{-1}
    \\
        & =\rho_{r} c^{2} \mathbf{U} \otimes \mathbf{U}+
    2 g^{-1} (T\Psi)^{\star}\frac{\partial E}{\partial \mathbf{H}} (T\Psi) g^{-1}
    -  E\,g^{-1},
  \end{align*}
  where the second equality is gained by lemma~\ref{lem:h-H}, and the third one by~\eqref{eq:def-h}. Finally, since $T\Psi.\mathbf{U}=0$ and $\mathbf{U}$ is a unit timelike vector field, we get $ \mathbf{S} \cdot \mathbf{U}^{\flat}=E\, \mathbf{U}$, and $\mathbf{T} \cdot \mathbf{U}^{\flat}=-L \mathbf{U}$.
\end{proof}

\begin{exam}[Relativistic perfect fluid]\label{rem:RF}
  The stress-energy tensor of a Relativistic perfect fluid,
  \begin{equation*}
    \mathbf{T} = \left(L  + P \right)\mathbf{U} \otimes \mathbf{U} + P\, g^{-1},
    \qquad
    L=\rho_{r} c^{2}+ E,
  \end{equation*}
  corresponds to an internal energy density of the form $E= \rho_{r}e(\rho_{r})$, where $P = \rho_{r}^{2}  e'(\rho_{r})$ is the pressure. Indeed, in that case, we have
  \begin{equation*}
    \frac{\partial e}{\partial \mathbf{H}} = e^{\prime}(\rho_{r})\frac{\partial \rho_{r}}{\partial \mathbf{H}} \quad \text{with} \quad  \frac{\partial \rho_{r}}{\partial \mathbf{H}} = \frac{1}{2} \rho_{r} \; \mathbf{H}^{-1},
  \end{equation*}
  and thus, by lemma~\ref{lem:h-H}, we get
  \begin{equation*}
    (T\Psi)^{\star}\frac{\partial  e}{\partial \mathbf{H}} (T\Psi) = \frac{1}{2}\rho_{r}e^{\prime}(\rho_{r}) (T\Psi)^{\star}\mathbf{H}^{-1}(T\Psi) = \frac{1}{2}\rho_{r}e^{\prime}(\rho_{r}) h.
  \end{equation*}
  Therefore
  \begin{equation*}
    {\bm{\Sigma}} = - 2\rho_{r}\, g^{-1} (T\Psi)^{\star}\frac{\partial  e}{\partial \mathbf{H}} (T\Psi) g^{-1} = - \rho_{r}^{2}  e'(\rho_{r})\, g^{-1} h g^{-1} =  - P (g^{-1} + \mathbf{U}\otimes\mathbf{U}),
  \end{equation*}
  where we have set $P := \rho_{r}^{2}  e'(\rho_{r})$, and we get
  \begin{equation*}
    \mathbf{S} = {\bm{\Sigma}} - E\,\mathbf{U} \otimes \mathbf{U} = -(E + P)\mathbf{U}\otimes\mathbf{U} - P \,g^{-1}.
  \end{equation*}
  The corresponding stress--energy tensor is thus given by
  \begin{equation*}
    \mathbf{T} = \rho_{r} c^{2} \mathbf{U} \otimes \mathbf{U} - \mathbf{S} = L \mathbf{U} \otimes \mathbf{U} - {\bm{\Sigma}} = (\rho_{r} c^{2} + E + P)\mathbf{U}\otimes\mathbf{U} + P \,g^{-1}.
  \end{equation*}
\end{exam}

Even if the full theory is four-dimensional, the orthogonal decomposition~\eqref{eq:def-Sigma} of $\mathbf{T}$ relative to $\mathbf{U}$, and the definition~\eqref{eq:dedH} (\emph{i.e.}, the Relativistic Hyperelasticity law) naturally introduce a three-dimensional symmetric stress tensor, either \emph{covariant},
\begin{equation}\label{eq:def-s-3D}
  \bs:=- 2\frac{\partial  e}{\partial \mathbf{H}} ,
\end{equation}
or, \emph{contravariant},
\begin{equation*}
  \mathbf{H}\,\bs\,\mathbf{H}=2\frac{\partial  e\;\,}{\partial \mathbf{H}^{-1}},
\end{equation*}
since the conformation is invertible (and contravariant). The stress tensors $\bs$ and $\mathbf{H}\,\bs\,\mathbf{H}$ are generalizations of the \emph{second Piola--Kirchhoff stress tensor} (expressed on a reference configuration $\Omega_{0}$ of Classical Continuum Mechanics) or more precisely here of the \emph{Rougée stress tensor} \cite{Rou1991a,Rou2006} (defined on the body $\mathcal{B}$, see~\autoref{sec:stress-on-the-body}).
These constitutive equations are the three-dimensional Relativistic Hyperelasticity laws. They do not depend on the further assumption of a foliation of the World tube $\mathcal{W}$, nor on the consideration of a spacetime.

The underlying question~\cite{EBT2006,Bro2021} is then how to properly import in General Relativity existing Classical Continuum Mechanics constitutive laws formulated on the body $\mathcal{B}$ \cite{Rou2006,KD2021} (or a reference configuration $\Omega_{0}$). Indeed, many three-dimensional expressions of energy densities
\begin{equation}\label{eq:w3D}
  w=w(\widehat{{\bm{\gamma}}}), \quad \widehat{{\bm{\gamma}}}:={\bm{\gamma}}_{0}^{-1}{\bm{\gamma}}, \quad \text{on $\mathcal{B}$}
  \qquad \Big( \emph{i.e.}, \;w=w(\widehat\mathbf{C}), \; \widehat\mathbf{C}:=q^{-1}\mathbf{C}, \; \text{when $\mathcal{B} \equiv \Omega_{0}$} \Big),
\end{equation}
are available in the Classical Continuum Mechanics literature \cite{Moo1940,Har1966,Ogd1972,AB1993,Sto2009,GMDC2011}.
They are local function of the mixed tensor $\widehat{{\bm{\gamma}}} = {\bm{\gamma}}_{0}^{-1}{\bm{\gamma}}$, defined using the reference metric ${\bm{\gamma}}_{0}$ on $\mathcal{B}$ (equivalently, of the mixed right Cauchy--Green tensor $\widehat\mathbf{C}$ on $\Omega_{0}$), meaning that
\begin{equation*}
  w(\mathbf{X})=w\left(\widehat{\bm{\gamma}}(\mathbf{X}) \right), \quad \mathbf{X} \in \mathcal{B}.
\end{equation*}

The use of such energy densities is then straightforward in the Relativistic framework, if one sets (using definition~\eqref{eq:H0})
\begin{equation}\label{eq:e-from-w}
  e(\Psi, \mathbf{H})=w\left(({\bm{\gamma}}_{0}^{-1}\circ \Psi\right)\mathbf{H}^{-1})=w\left(\mathbf{H}_{0}\,\mathbf{H}^{-1}\right).
\end{equation}
Using~\eqref{eq:e-from-w}, we get
\begin{equation}\label{eq:s-3D}
  \bs = -2\frac{\partial e}{\partial \mathbf{H}}=2\,\mathbf{H}^{-1} \mathbf{H}_{0}\frac{\partial  w}{\partial \widehat{{\bm{\gamma}}}}\mathbf{H}^{-1},
  \qquad
  \mathbf{H}\, \bs \, \mathbf{H}=2\, \mathbf{H}_{0}\frac{\partial  w}{\partial \widehat{{\bm{\gamma}}}},
\end{equation}
so that the stress-energy tensor $\mathbf{T}$ and the four-dimensional stress ${\bm{\Sigma}}$ recast finally  as
\begin{equation}\label{eqTs-3D}
  \mathbf{T}= L \mathbf{U} \otimes \mathbf{U}-{\bm{\Sigma}}, \qquad {\bm{\Sigma}}=\rho_{r}\, g^{-1} (T\Psi)^{\star}\, \bs\, (T\Psi) g^{-1}.
\end{equation}
We refer to \autoref{sec:stress-on-the-body}  for the full link ---which needs the consideration of a spacetime--- with stresses on the body $\mathcal{B}$.

\section{Universe's foliation by spacelike hypersurfaces}
\label{sec:Universe-foliation}

There is no Mechanics without the proper definition of \emph{time} and \emph{space}. To introduce these concepts in General Relativity, one usually starts by introducing a smooth submersion (a \emph{time function}) $\hat{t}$ on the Universe $\mM$ with a timelike gradient everywhere. Then, spacelike hypersurfaces are defined as
\begin{equation}\label{eq:Omegatini}
  \espace_{t} := \set{m \in \mM;\; \hat{t}(m) = t},
\end{equation}
and one expects the Universe to be foliated by these hypersurfaces~\cite{Lic1955,Gou2012}. The problem is that, in general, a global foliation of the Universe might not exist (see~\cite[Chapter 4]{Gou2012}). Anyway, if such a foliation exists or is given, one will say that the Universe has been endowed with a \emph{spacetime structure} or a \emph{$(3+1)$-structure} as defined in~\cite{Dar1927,Lic1939,FBCB1948,Lic1952,Fou1956,ADM1962,Yor1979,Gou2012}.

Fortunately, for our concerns, we do not have to address this problem \emph{globally}. In the present paper, we will simply admit that such a foliation exists on a \emph{local chart} which contains the body World tube $\mathcal{W}$, or a part of it. Indeed, in Mechanical Engineering, the spacetime domain occupied by a continuous medium/a structure, embedded for example in a laboratory, a building, a city, a country, a domain of space \ldots can be considered as included into such a local chart.

Moreover, since the presence of the studied matter in the laboratory does not affect (much) the Universe metric $g$ compared to the one of Earth (\emph{passive matter} assumption), we can choose, among the numerous spacetimes encountered in General Relativity and available in~\cite{HE1973,MTW1973,MG2009}, those describing solutions of Einstein equations in the vacuum. These spacetimes are usually described using a coordinate system $(x^{\mu})$, for which the time function is chosen as
\begin{equation*}
  \hat{t}=\frac{x^{0}}{c},
\end{equation*}
and we then define $\Omega_{t}$ as the intersection of the World tube $\mathcal{W}$ with the spacelike hypersurface $\espace_{t}$,
\begin{equation*}
  \Omega_{t} := \mathcal{W} \cap \espace_{t} = \mathcal{W} \cap \set{x^{0}=ct}.
\end{equation*}
The three-dimensional hypersurfaces $\Omega_{t}$ of the Universe play the same role as the configurations, parameterized by time $t$, of Classical Continuum Mechanics~\cite{TN1965,Nol1972,Nol1974,MH1994,Ber2012,Ste2015,For2022}, with the difference that the later are embedded in the three-dimensional Euclidean space $\espace$, not in the four-dimensional Universe $\mM$. This construction is illustrated in \autoref{fig:foliationOmegat}, where a second time, $t_{0}$, and the associated hypersurface $\Omega_{t_{0}}$ (a possible reference configuration) are represented.

\begin{figure}[ht]
  \centering
  \includegraphics[width=12cm]{Fig2.pdf}
  \caption{The foliation of the World tube $\mathcal{W}$ by spacelike hypersurfaces $\Omega_{t}$.}
  \label{fig:foliationOmegat}
\end{figure}

The canonical embedding of these submanifolds into the Universe $\mM$ is noted $j_{t}\colon \Omega_{t} \to \mM$ (rather than $\theta$ as in \cite[Chapter 7]{HE1973}). Then, for each $t$, the pullback of the matter field $\Psi$ by~$j_{t}$,
\begin{equation*}
  \Psi_{t} := j_{t}^{*} \Psi = \Psi \circ j_{t} \colon \Omega_{t} \to V,
\end{equation*}
is just the restriction of $\Psi$ to $\Omega_{t}$, and we have $\Psi_{t}(\Omega_{t})=\mathcal{B}$.

\begin{rem}\label{rem:TPsit-iso}
  We have not made, so far, the assumption that $\Psi_{t}$ is a diffeomorphism. However, the tangent map
  \begin{equation*}
    T_{\mathbf{x}} \Psi_{t} \colon T_{\mathbf{x}} \Omega_{t} \to T_{\Psi_{t}(\mathbf{x})}\mathcal{B}
  \end{equation*}
  is an isomorphism for each $\mathbf{x} \in \Omega_{t}$ and each $t$, since we have assumed that $\Psi$ is a submersion. We shall denote the inverse mapping of $T\Psi_{t}$ by $\mathbf{F}(t)$ and omit, when not necessary, the explicit dependence on time and write simply $\mathbf{F} = (T\Psi_{t})^{-1}$. If we make, furthermore, \emph{the stronger hypothesis} that $\Psi_{t}$ is a diffeomorphism we can set
  \begin{equation*}
    p(t)= \Psi_{t}^{-1}\colon \mathcal{B} \to \Omega_{t},
  \end{equation*}
  and then $\mathbf{F}=Tp$.
\end{rem}

The unit normal to the spacelike hypersurfaces $\Omega_{t}$, proportional to the gradient of the time function $\hat{t}$, is denoted by the quadrivector $\mathbf{N}$. We have two opposite choices to define such a unit vector and we define $\mathbf{N}$ as~\cite{Gou2012}
\begin{equation}\label{eq:defN}
  \mathbf{N} := -\frac{\grad\, \hat{t}}{\sqrt{-\norm{\grad\, \hat{t}}^{2}}},
\end{equation}
where the gradient and the norm are relative to the metric $g$. The minus sign is chosen so that the quadrivector $\mathbf{N}$ is \emph{future-oriented}, meaning that the value of $\hat{t}$ increases along the flow curves of $\mathbf{N}$.
Note that, at each point $m \in \mM$ where the spacetime structure is defined, we have the orthogonal decomposition
\begin{equation*}
  T_{m}\mM = \langle \mathbf{N}(m) \rangle \oplus \mathbf{N}(m)^{\bot},
\end{equation*}
where, for each $m\in \Omega_{t}$, the orthogonal subspace $\mathbf{N}(m)^{\bot} = T_{m}\Omega_{t}$ coincides with the tangent space at $m$ to the spacelike hypersurface $\Omega_{t}$.

An important special case, and the only one used in this paper for our description of Relativistic Continuum Mechanics of solids, is the one of a \emph{static spacetime}. Such a spacetime is induced by a \emph{static metric}, \emph{i.e.}, a metric $g$ for which there exists a timelike Killing vector field $X$ (\textit{i.e.} $\Lie_{X}g=0$), which is moreover the gradient of a time function $\hat{t}$. There exists then a coordinate system $(x^{\mu})$, with $x^{0}=ct$ the time coordinate, for which
\begin{equation*}
  \frac{\partial g_{\mu \nu}}{\partial x^{0}} = 0 \quad \text{and} \quad g_{0i}=0.
\end{equation*}
In that case, we get
\begin{equation*}
  \grad\, \hat{t} = \frac{1}{c} g^{00}\partial_{x^{0}},
\end{equation*}
and the unit normal $\mathbf{N}$ is written as
\begin{equation*}
  \mathbf{N} = \sqrt{-g^{00}}\partial_{x^{0}}.
\end{equation*}
Examples of static spacetimes are the Minkowski and the Schwarzschild~\cite{Sch1916} spacetimes.

\section{Matter field in a spacetime -- Generalized Lorentz factor}
\label{sec:Matter-field-spacetime}

Perfect matter in the Universe $\mM$ is represented by a matter field $\Psi$. This field generates a timelike quadrivector $\mathbf{P}$ on the World tube $\mathcal{W}$, as introduced in \autoref{sec:matter-fields}, and a unit timelike quadrivector $\mathbf{U}=\mathbf{P}/\rho_{r}$. Therefore, if a spacetime structure is introduced on the Universe as it has been explained in \autoref{sec:Universe-foliation}, we get a second unit timelike quadrivector $\mathbf{N}$, normal to the hypersurfaces $\Omega_{t}$, and in general not collinear to $\mathbf{U}$. By changing the sign of the time function $\hat{t}$ if necessary, we can assume, anyway, that both $\mathbf{U}$ and $\mathbf{N}$ define the same time orientation. This is characterized by the condition
\begin{equation*}
  \langle \mathbf{U}, \mathbf{N} \rangle_{g} < 0.
\end{equation*}
Thus, the quadrivector $\mathbf{U}$ can be written uniquely using the orthogonal decomposition
\begin{equation*}
  \mathbf{U} = \mathbf{U}^{N} + \mathbf{U}^{\top},
\end{equation*}
where $\mathbf{U}^{N} = -\langle \mathbf{U}, \mathbf{N} \rangle_{g}\,\mathbf{N}$ is the normal component of $\mathbf{U}$ and $\mathbf{U}^{\top}$ is the tangential component to $\Omega_{t}$. Introducing the function (see for instance~\cite{Gou2006})
\begin{equation}\label{eq:def-gamma}
  \gamma: = - \langle \mathbf{U}, \mathbf{N}\rangle_{g},
\end{equation}
one can write thus,
\begin{equation}\label{eq:U-orthogonal-decomposition}
  \mathbf{U}^{N} = \gamma \mathbf{N}, \quad \text{and} \quad \mathbf{U}^{\top} = \frac{\gamma}{c} \bm{u},
\end{equation}
where $c$ is the light speed, and $\bm{u} \in \mathbf{N}^{\bot} = T\Omega_{t}$.

\begin{rem}
  Since we deal with a foliation by hypersurfaces $\Omega_{t}$, rather than just one hypersurface, the tangential component $\mathbf{U}^{\top}$ of a vector field $\mathbf{U}$ defined on $\mM$ (or an open subset of $\mM$) can be simultaneously interpreted as a vector field defined on $\mM$ but tangential to $\Omega_{t}$ at each point $m$, or as a (time-dependent) vector field on $\Omega_{t}$ (when restricted to $\Omega_{t}$). In our notations, we do not distinguish between these two interpretations.
\end{rem}

The orthogonal decomposition $\gamma \mathbf{N} + \frac{\gamma}{c} \bm{u}$ of $\mathbf{U}$ and the relation $T\Psi. \mathbf{U}=0$ deduced from~\eqref{eq:TPsi-P-null} allow to express the three-dimensional velocity $\bm{u}$ on $\Omega_{t}$, as
\begin{equation}\label{eq:def-u}
  \bm{u} = -c\, \mathbf{F}\, T\Psi. \mathbf{N},
  \qquad
  \mathbf{F}^{-1} := T\Psi_{t},
\end{equation}
where $\mathbf{F}\colon T\mathcal{B} \to T \Omega_{t}$ is defined as the inverse of $T\Psi_{t}$, the restriction of $T\Psi$ to $T\Omega_{t}$, which is an invertible linear mapping (remark~\ref{rem:TPsit-iso}). In \autoref{sec:Gallilean}, the expression~\eqref{eq:def-u} will allow us to interpret the Galilean limit of $\bm{u}$ as the three-dimensional Eulerian velocity on $\Omega_{t}$.

Using the fact that $\norm{\mathbf{U}}^{2}_{g}=-1$, one gets furthermore that
\begin{equation*}
  \gamma = \frac{1}{\sqrt{1-\displaystyle\frac{\norm{\bm{u}}_{g}^{2}}{c^{2}}}},
\end{equation*}
and $\gamma\geq 1$ since $\bm{u}$ is spacelike. This function plays a fundamental role in General Relativity and its notation is not accidental. In the special case of the Minkowski spacetime, where $g=\eta$ is the Minkowski metric, one recovers the traditional Lorentz factor
\begin{equation*}
  \gamma = \frac{1}{\sqrt{1-\frac{\bm{u}^{2}}{c^{2}}}},
\end{equation*}
where $\bm{u}^{2}:=\norm{\bm{u}}_{q}^{2}$ is the square Euclidean norm of $\bm{u}$. For this reason, we shall call $\gamma = - \langle \mathbf{U}, \mathbf{N}\rangle_{g}$ the \emph{generalized Lorentz factor}.

\begin{rem}[Rest frame and observers]
  The concept of \emph{rest frame} is well-defined for a particle in Special Relativity. Its definition for distributed matter in general Relativity is much less clear. In this paper, we will adopt the following definition. Given a matter field $\Psi$, a \emph{rest frame} will be defined as a spacetime in which $\mathbf{U} = \mathbf{N}$, \emph{i.e.}, as a spacetime in which the generalized Lorentz factor is $\gamma=1$. For such a spacetime, we will get of course $\bm{u} = 0$ and the particles can be considered at rest in it. The corresponding time coordinate will thus be interpreted as the \emph{proper time} $\tau$. More generally and heuristically, we can interpret $\mathbf{U}$ as ``inducing a splitting of the tangent bundle $T\mathcal{W}$ for matter'' and $\mathbf{N}$ as ``inducing an integrable splitting or (3+1) spacetime for an observer''. The Lorentz factor $\gamma = -\langle \mathbf{U}, \mathbf{N}\rangle_{g}$ is then the ``angle'' between these two timelike directions.
\end{rem}

Finally, the normal component of the current of matter $\mathbf{P}=\rho_{r} \mathbf{U}$ (definition~\eqref{eq:def-U}) is then simply
\begin{equation*}
  \mathbf{P}^{N}=\gamma \rho_{r}\, \mathbf{N}= \rho\, \mathbf{N},
\end{equation*}
where $\mathbf{N}$ is the unit timelike normal to the hypersurfaces $\Omega_{t}$.
The function
\begin{equation}\label{eq:rho}
  \rho:= \gamma \rho_{r},
\end{equation}
defined on the World tube $\mathcal{W}$, is interpreted as the \emph{relativistic mass density}. A geometric interpretation of the restriction of $\rho$ to $\Omega_{t}$ is provided in \autoref{sec:3D-metrics}.

\section{Relativistic stress tensors and constitutive laws in a spacetime}
\label{sec:hyperelasticity-spacetime}

We assume here that the World tube $\mathcal{W}$ is foliated by three-dimensional hypersurfaces $\Omega_{t}$, with unit timelike normal $\mathbf{N}$. Then, it is possible to use the orthogonal decomposition of each tangent space relative to $\mathbf{N}$ to split each tensor field accordingly. These splittings are referred to as \emph{(3+1)-decompositions} in the General Relativity literature~\cite{ADM1962,Yor1979,Gou2012}. Explicit formulas for second-order tensors are provided in \autoref{sec:Orth-Decomp-2nd}. We follow here the calculations of Souriau~\cite{Sou1958,Sou1964} for the flat Minkowski spacetime and extend them to any spacetime, thanks to this (3+1)-decomposition. These calculations generalize the ones given for relativistic fluids in~\cite{Gou2012} to relativistic solids. In particular, the orthogonal decomposition of the stress-energy tensor allows us to introduce relativistic generalizations of the \emph{Cauchy stress tensor} as 3D tensors and to recast the 4D Relativistic Hyperelasticity law (theorem~\ref{thm:relativistic-hyperelasticity}) as a \emph{three-dimensional constitutive equation}, relating these 3D stress tensors to the conformation $\mathbf{H}$.

The orthogonal decompositions of $g$ and $g^{-1}$, relative to the unit timelike vector $\mathbf{N}$ (instead of $\mathbf{U}$ as in \autoref{sec:Conformation-strains}) are
\begin{equation}\label{eq:g3D}
  g = {g^{3D}} - \mathbf{N}^{\flat} \otimes \mathbf{N}^{\flat},
  \qquad
  g^{-1} = ({g^{3D}})^{\sharp} - \mathbf{N} \otimes \mathbf{N},
  \qquad
  {g^{3D}} \mathbf{N} = 0,
\end{equation}
where the degenerate quadratic form ${g^{3D}}$ and $({g^{3D}})^{\sharp}=g^{-1} {g^{3D}} g^{-1}$ are of signature $(0, +, +, +)$, by lemma~\ref{lem:h-signature}. This decomposition allows, in particular, to recast the conformation as
\begin{equation*}
  \mathbf{H} = T\Psi \, g^{-1} \, (T\Psi)^{\star} = T\Psi (({g^{3D}})^{\sharp} - \mathbf{N} \otimes \mathbf{N}) (T\Psi)^{\star},
\end{equation*}
with
\begin{equation*}
  T\Psi (\mathbf{N} \otimes \mathbf{N}) (T\Psi)^{\star} = T\Psi \left(\frac{1}{\gamma} \mathbf{U}-\frac{1}{c} \bm{u} \right) \otimes \left(\frac{1}{\gamma} \mathbf{U}-\frac{1}{c} \bm{u}\right) (T\Psi)^{\star} = \mathbf{F}^{-1} \left(\frac{1}{c^{2}} \bm{u} \otimes \bm{u}\right) \mathbf{F}^{-\star},
\end{equation*}
since $T\Psi. \mathbf{U}=0$, and thus
\begin{equation}\label{eq:conformation-in-spacetime}
  \mathbf{H} = \mathbf{F}^{-1}\left( ({g^{3D}})^{\sharp} -\frac{1}{c^{2}} \bm{u} \otimes \bm{u}\right) \mathbf{F}^{-\star}.
\end{equation}

When applied to the stress-energy tensor $\mathbf{T}$, the orthogonal decomposition~\eqref{eq:T-orthogonal-decomposition} leads to
\begin{equation}\label{eq:stress-energy-tensor-decomposition}
  \mathbf{T} = E_{\mathrm{tot}} \, \mathbf{N} \otimes \mathbf{N} + \frac{1}{c} \left(\mathbf{N} \otimes \bm{p} + \bm{p} \otimes \mathbf{N} \right) + \mathbf{t},
\end{equation}
and allows to define the physical components of $\mathbf{T}$ in the considered spacetime:
\begin{itemize}
  \item $E_{\mathrm{tot}}$ is the \emph{total energy density},
  \item $\bm{p}$ is the \emph{momentum density vector field},
  \item and the spatial part $\mathbf{t}$ of $\mathbf{T}$ is related to the \emph{stress field}.
\end{itemize}
The question asked by Souriau is then: How are these quantities related to a three-dimensional relativistic generalization of the \emph{Cauchy stress tensor} ${\bm{\sigma}}$? As shown below, the answer depends on the choice of the decomposition  of the stress-energy tensor (see theorem~\ref{thm:relativistic-hyperelasticity} and remark~\ref{rem:T}). Indeed, we have seen  that there are two possible splittings of it:
\begin{enumerate}
  \item $\mathbf{T} = \rho_{r}c^{2} \mathbf{U} \otimes \mathbf{U} - \mathbf{S}$ (considered by Souriau \cite{Sou1958,Sou1964} and Synge \cite{Syn1959}),
  \item or, $\mathbf{T} = L \mathbf{U} \otimes \mathbf{U} - {\bm{\Sigma}}$ (considered by Eckart \cite{Eck1940} and Bennoun \cite{Ben1965}).
\end{enumerate}
The orthogonal decompositions (relative to $\mathbf{N}$) of the two four-dimensional stresses $\mathbf{S}$ and ${\bm{\Sigma}}$ interestingly give rise to two possible ways to define a three-dimensional stress tensor ${\bm{\sigma}}$:
\begin{enumerate}
  \item either as the spatial part of $\mathbf{S}$,
  \item or, as the spatial part of ${\bm{\Sigma}}$.
\end{enumerate}

\subsection*{First choice: $\pmb \sigma$ is defined as the spatial part of $\mathbf{S}$}

Using the fact that $\mathbf{S}\cdot\mathbf{U}^{\flat}=E\, \mathbf{U}$ by theorem~\ref{thm:relativistic-hyperelasticity} with $\mathbf{U}^{\flat} = \gamma \mathbf{N}^{\flat} + \gamma \bm{u}^{\flat}/ c$ by~\eqref{eq:U-orthogonal-decomposition}, the orthogonal decomposition of $\mathbf{S}$ can be expressed as
\begin{equation}\label{eq:S-from-sigma}
  \mathbf{S} = \alpha \, \mathbf{N} \otimes \mathbf{N} + \mathbf{N} \otimes \bm{a} +\bm{a} \otimes \mathbf{N} + {\bm{\sigma}},
\end{equation}
where the spatial part of $\mathbf{S}$ has been set equal to ${\bm{\sigma}}$, and
\begin{equation*}
  \bm{a} = \frac{1}{c}\left({\bm{\sigma}} \cdot \bm{u}^{\flat} - E\bm{u}\right), \quad \text{and} \quad \alpha =  \frac{1}{c^{2}}\bm{u}^{\flat} \cdot {\bm{\sigma}} \cdot \bm{u}^{\flat} - E\left( 1 + \frac{\norm{\bm{u}}^{2}}{c^{2}}\right).
\end{equation*}
The associated orthogonal decomposition of the stress-energy tensor $\mathbf{T} = \rho_{r}c^{2} \mathbf{U} \otimes \mathbf{U} - \mathbf{S}$ is then
\begin{equation*}
  \begin{cases}
    E_{\mathrm{tot}} = \gamma \rho c^{2} + E\Big(1+\frac{1}{c^{2}}  \norm{\bm{u}}^{2}\Big) -\frac{1}{c^{2}} \bm{u}^{\flat} \cdot {\bm{\sigma}} \cdot \bm{u}^{\flat},
    \\
    \bm{p} =(\gamma \rho c^{2}+E)\bm{u} - {\bm{\sigma}} \cdot \bm{u}^{\flat},
    \\
    \mathbf{t} = \gamma\rho \, \bm{u} \otimes \bm{u} - {\bm{\sigma}},
  \end{cases}
\end{equation*}
where $\rho= \gamma \rho_{r}$ is the relativistic mass density. The three-dimensional stress tensor ${\bm{\sigma}}$, defined as the spatial part of
\begin{equation}\label{eq:Hyperelasticity-law-S}
  \mathbf{S} = - 2\rho_{r}\, g^{-1} (T\Psi)^{\star}\frac{\partial e}{\partial \mathbf{H}} (T\Psi) g^{-1}- E\, \mathbf{U} \otimes \mathbf{U},
\end{equation}
is thus
\begin{equation}\label{eq:Hyperelasticity-a}
  {\bm{\sigma}} = - \frac{2}{\gamma} \rho\, ({g^{3D}})^{\sharp} \mathbf{F}^{-\star}\frac{\partial e}{\partial \mathbf{H}} \mathbf{F}^{-1} ({g^{3D}})^{\sharp} - \frac{\gamma^{2}E}{c^{2}} \, \bm{u} \otimes \bm{u}.
\end{equation}
The later equation can be interpreted as a three-dimensional Hyperelasticity law. Introducing the generalized second Piola--Kirchhoff stress tensor~\eqref{eq:def-s-3D}, we get
\begin{equation*}
  {\bm{\sigma}} = \frac{1}{\gamma} \rho\, ({g^{3D}})^{\sharp} \mathbf{F}^{-\star}\,\bs\,\mathbf{F}^{-1} ({g^{3D}})^{\sharp} - \frac{\gamma^{2}E}{c^{2}} \, \bm{u} \otimes \bm{u},
  \qquad
  \bs=\displaystyle-2\frac{\partial e}{\partial \mathbf{H}}.
\end{equation*}

\subsection*{Second choice: $\pmb \sigma$ is defined as the spatial part of $\pmb\Sigma$}

Using this time the fact that ${\bm{\Sigma}} \cdot \mathbf{U}^{\flat}=0$ by remark~\ref{rem:T}, we get the orthogonal decomposition
\begin{equation}\label{eq:Sigma-from-sigma}
  {\bm{\Sigma}} =  \frac{1}{c^{2}}(\bm{u}^{\flat} \cdot {\bm{\sigma}} \cdot \bm{u}^{\flat}) \, \mathbf{N} \otimes \mathbf{N}
  + \frac{1}{c} \left(\mathbf{N} \otimes ( {\bm{\sigma}} \cdot \bm{u}^{\flat}) +( {\bm{\sigma}} \cdot \bm{u}^{\flat}) \otimes \mathbf{N} \right)
  + {\bm{\sigma}},
\end{equation}
where the spatial part of ${\bm{\Sigma}}$ has been set equal to ${\bm{\sigma}}$. The associated orthogonal decomposition of the stress-energy tensor $\mathbf{T} = L \mathbf{U} \otimes \mathbf{U} - {\bm{\Sigma}}$, with $L=\rho_{r}c^{2}  + E$, is now
\begin{equation*}
  \begin{cases}
    E_{\mathrm{tot}} =  \gamma^{2} L -\frac{1}{c^{2}} \bm{u}^{\flat} \cdot {\bm{\sigma}} \cdot \bm{u}^{\flat}
    ,
    \\
    \bm{p} = \gamma^{2} L\, \bm{u} - {\bm{\sigma}} \cdot \bm{u}^{\flat},
    \\
    \mathbf{t} = \frac{\gamma^{2}}{c^{2}} L\, \bm{u} \otimes \bm{u} - {\bm{\sigma}} ,
  \end{cases}
\end{equation*}
where, using $\gamma^{2}=1+\frac{\gamma^{2}}{c^{2}} \norm{\bm{u}}^{2}$,
\begin{equation*}
  \gamma^{2} L = \gamma \rho c^{2} + E\Big(1+\frac{\gamma^{2}}{c^{2}} \norm{\bm{u}}^{2}\Big) .
\end{equation*}
The three-dimensional stress tensor ${\bm{\sigma}}$, defined as the spatial part of
\begin{equation}\label{eq:Hyperelasticity-law-Sigma}
  {\bm{\Sigma}} = - 2\rho_{r}\, g^{-1} (T\Psi)^{\star}\frac{\partial e}{\partial \mathbf{H}} (T\Psi) g^{-1}.
\end{equation}
is then
\begin{equation}\label{eq:Hyperelasticity-b}
  {\bm{\sigma}} := - \frac{2}{\gamma} \rho\, ({g^{3D}})^{\sharp} \mathbf{F}^{-\star}\frac{\partial  e}{\partial \mathbf{H}} \mathbf{F}^{-1} ({g^{3D}})^{\sharp},
\end{equation}
a relation which can be interpreted as a three-dimensional Hyperelasticity law. Introducing~\eqref{eq:def-s-3D}, we end up with
\begin{equation}\label{eq:Hyperelasticity-b-bis}
  {\bm{\sigma}} :=  \frac{1}{\gamma} \rho\, ({g^{3D}})^{\sharp} \mathbf{F}^{-\star}\,\bs\, \mathbf{F}^{-1} ({g^{3D}})^{\sharp},
  \qquad
  \bs=\displaystyle-2\frac{\partial e}{\partial \mathbf{H}}.
\end{equation}

\bigskip

Conversely, once the three-dimensional generalized Cauchy stress tensor ${\bm{\sigma}}$ is given (through a three-dimensional Hyperelasticity law), the four-dimensional stress tensors $\mathbf{S}$ and ${\bm{\Sigma}}$ are then fully determined, either by~\eqref{eq:S-from-sigma} or by~\eqref{eq:Sigma-from-sigma}. Even if the full theory is four-dimensional, \emph{the Relativistic Hyperelasticity laws are three-dimensional}. Finally, observe also that the difference between~\eqref{eq:Hyperelasticity-a} and~\eqref{eq:Hyperelasticity-b} is purely relativistic: both of them converge to the same three-dimensional stress tensor ${\bm{\sigma}}$ at the Galilean limit $c\to 0$ if one assumes $\lim_{c\to \infty}(E/c^{2})=0$.

\begin{rem}\label{rem:Kirchhoff}
  The three-dimensional stress tensor,
  \begin{equation*}
    {\bm{\tau}}:=\frac{{\bm{\sigma}}}{\rho_{r}}=\gamma\frac{{\bm{\sigma}}}{\rho},
    \quad\text{such as} \quad
    {\bm{\tau}} =({g^{3D}})^{\sharp} \mathbf{F}^{-\star}\,\bs\, \mathbf{F}^{-1} ({g^{3D}})^{\sharp},
  \end{equation*}
  is a first Relativistic generalization of the Kirchhoff stress tensor ${\bm{\sigma}}/\rho$ of Classical Continuum Mechanics (see \autoref{sec:stress-on-the-body} for a second generalization relative to the Schwarzschild spacetime).
\end{rem}

\section{Relativistic Hyperelasticity in Schwarzschild spacetime}
\label{sec:hyperelasticity-Schwarzschild}

In~\cite{Sou1958,Sou1960,Sou1964}, Souriau has provided a full description of Relativistic Hyperelasticity in the Minkowski spacetime, making implicitly the \emph{passive matter hypothesis}, meaning that the matter field $\Psi$ under study is negligible as a source of the gravitation field. Minkowski spacetime is a flat static spacetime, with no gravitational source, which is described in some coordinate system $(x^{0}=ct, x^{i})$ by the constant metric
\begin{equation*}
  \eta = -(dx^{0})^{2} + \delta_{ij}dx^{i}dx^{j}.
\end{equation*}
Of course, this situation is not fully realistic. However, due to the fact that for any point $m_{0}$ of the Universe, it is always possible to find a chart around $m_{0}$ in which the Christoffel symbols vanish at $m_{0}$, we can assume that the Christoffel symbols are almost zero around this point. This situation exactly corresponds to a free fall (like inside an orbital station), it approximately corresponds to mechanical situations on Earth surface for which gravity can be neglected or not taken into account.

Our goal here is to extend Souriau's results on Relativistic Hyperelasticity by taking into account gravity. These results will be used in the next section to formulate classical Galilean Hyperelasticity with gravity (or Newton--Cartan theory of Continuum Mechanics~\cite{Car1923,Car1924,Car1925}). To do so, rather than using the Minkowski spacetime, we shall assume here that the continuous medium/the structure considered is embedded in the \emph{Schwarzschild spacetime}~\cite{Sch1916,MTW1973,MG2009} (see \autoref{sec:Schwarzschild-spacetime}). Therefore, we neglect the influence of the matter under study (passive matter hypothesis) as a source of the gravity field. The \emph{exterior Schwarzschild metric} is a static solution $g$ of Einstein equation $\Ein_{g} = 0$ in the vacuum (with vanishing cosmological constant $\Lambda=0$). It is representative of the gravity field around a spherical and nonrotating planet (or a star or a black hole) of mass $M$ and radius $r_{0}$, such as the Earth. In this model, the rotation of the celestial body as a potential source of the gravitation field has been neglected. An alternative would have been to choose the \emph{Kerr metric}~\cite{MG2009} rather than the Schwarzschild metric, another possible choice which we did not make. A practical consequence of this choice is that the frame deduced from the coordinate system in which is described the
Schwarzschild metric corresponds to one pointing at fixed stars (the Earth is assumed to be nonrotating).

In the so-called \emph{Cartesian isotropic coordinates} (see \autoref{sec:Schwarzschild-spacetime}), the Schwarzschild metric has for expression
\begin{equation*}
  g = -\left(\frac{\bar{r}-\bar{r}_{\mathrm{s}}}{\bar{r} + \bar{r}_{\mathrm{s}}}\right)^{2} c^{2} \mathrm{d} t^{2}
  + \left(1+\frac{\bar{r}_{\mathrm{s}}}{\bar{r}}\right)^{4} \delta_{ij} x^{i}x^{j},
  \qquad
  \bar{r}:=\sqrt{\delta_{ij} x^{i}x^{j}},
\end{equation*}
with $\bar{r}=0$ at the center of the planet and $\bar{r}\approx r_{0}$ on its surface. The \emph{reduced Schwarzschild radius}
\begin{equation*}
  \bar{r}_{\mathrm{s}} = \frac{1}{4} r_{\mathrm{s}} = \frac{GM}{2c^{2}}
\end{equation*}
depends on the gravitational constant $G$, the mass $M$ of the celestial body, and the light speed $c$. It is much smaller than the radius $r_{0}$ of the planet~\cite{MTW1973} ($\bar{r}_{\mathrm{s}}\approx 2$ mm for the Earth).

Introducing the \emph{lapse function}~\cite{ADM1962}
\begin{equation*}
  \mathcal{N}:=\frac{1}{\sqrt{-\norm{\grad^{g} x^{0}}^{2}}} = \frac{1}{\sqrt{-g^{00}}} = \sqrt{-g_{00}} =  \frac{1-\bar{r}_{\mathrm{s}}/\bar{r}}{1+\bar{r}_{\mathrm{s}}/\bar{r}},
\end{equation*}
the metric can be written as
\begin{equation}\label{eq:Schwarzschild-metric}
  g = - \mathcal{N}^{2} c^{2}  \mathrm{d} t^{2}+ {g^{3D}},
\end{equation}
where ${g^{3D}}$ is the spatial (conformal) metric,
\begin{equation}\label{eq:k}
  {g^{3D}} = k q,
  \qquad
  k = \left(1+\frac{\bar{r}_{\mathrm{s}}}{\bar{r}}\right)^{4}
  \qquad
  q = \delta_{ij} \mathrm{d} x^{i} \mathrm{d} x^{j}.
\end{equation}
The unit normal $\mathbf{N}$, defined by~\eqref{eq:defN}, to the three-dimensional spatial hypersurfaces $\Omega_{t}$ is simply
\begin{equation}\label{eq:Schwarzschild-normal}
  \mathbf{N} = \frac{1}{\mathcal{N}}\, \partial_{x^{0}}=\frac{1}{c\mathcal{N}}\, \partial_{t}.
\end{equation}

\begin{rem}
  The flat Minkowski spacetime (Special Relativity) is simply the limiting massless case $M=0$ and thus $\bar{r}_{\mathrm{s}} = r_{\mathrm{s}} = 0$. It is a special case of this more general framework.
\end{rem}

The Schwarzschild metric is not flat. The non-vanishing symmetric Christoffel symbols, in the Cartesian isotropic coordinate systems $(x^{\mu})$ and $(t, x^{i})$, can easily be recovered using the usual formula~\eqref{eq:Gammas} since ${g^{3D}}=k(\bar{r}) q$ is conformal (see also~\cite{MG2009}). They are written as
\begin{equation}\label{eq:Schwarzschild-Christoffels}
  \begin{split}
    \Gamma^{t}_{ti} & = \Gamma^{0}_{0i} = (\mathrm{d} \ln \mathcal{N})_{i} = 2 \frac{\bar{r}_{\mathrm{s}}}{\bar{r}^{2}} \frac{ 1}{\left(1-\bar{r}_{\mathrm{s}}^{2}/\bar{r}^{2}\right)} \delta_{ik} \frac{x^{k}}{\bar{r}},
    \\
    \Gamma^{i}_{tt} & = c^{2} \Gamma^{i}_{00} = c^{2}\frac{\mathcal{N}}{k} (\grad \mathcal{N})^{i} = 2c^{2} \frac{\bar{r}_{\mathrm{s}}}{\bar{r}^{2}} \frac{(1-\bar{r}_{\mathrm{s}}/\bar{r})}{(1+\bar{r}_{\mathrm{s}}/\bar{r})^{7}} \frac{x^{i}}{\bar{r}},
    \\
    \Gamma^{i}_{jj} & = -\frac{1}{2}(\grad \ln k)^{i} = 2 \frac{\bar{r}_{\mathrm{s}}}{\bar{r}^{2}} \frac{ 1}{\left(1+\bar{r}_{\mathrm{s}}/\bar{r}\right)} \frac{x^{i}}{\bar{r}}, \quad \text{for $i\ne j$},
    \\
    \Gamma^{j}_{ji} & = \frac{1}{2}(\mathrm{d} \ln k)_{i} = - 2 \frac{\bar{r}_{\mathrm{s}}}{\bar{r}^{2}} \frac{ 1}{\left(1+\bar{r}_{\mathrm{s}}/\bar{r}\right)}\delta_{ik}\frac{x^{k}}{\bar{r}}, \quad \text{for $i = j$ and $i\ne j$}.
  \end{split}
\end{equation}
with no sum on $j$, and where, for a function $f$, $\grad f=q^{-1} \mathrm{d} f$ is the gradient relative to the Euclidean metric $q=(\delta_{ij})$. The related divergence operators are detailed in~\autoref{sec:divergences-spacetime}.

In the following, we particularize the relations established in \autoref{sec:hyperelasticity-spacetime} for the special case of the Schwarzschild spacetime when expressed in Cartesian isotropic coordinates.

$\bullet$ The \emph{three-dimensional velocity} $\bm{u}$ defined by~\eqref{eq:def-u} is then given by
\begin{equation}\label{eq:u-s}
  \bm{u} = -\frac{1}{\mathcal{N}}\, \mathbf{F}\, \frac{\partial \Psi}{\partial t},
\end{equation}
where $\mathbf{F}=T\Psi_{t}^{-1}$ was introduced in~\eqref{eq:def-u}.

$\bullet$ The \emph{generalized Lorentz factor}~\eqref{eq:def-gamma} has then for expression
\begin{equation}\label{eq:gamma-s}
  \gamma = \frac{1}{\sqrt{1-k \frac{\bm{u}^{2}}{c^{2}}}},
\end{equation}
where $\bm{u}^{2}:=\norm{\bm{u}}_{q}^{2}$ is the Euclidean squared norm and $k=\left(1+\frac{\bar{r}_{\mathrm{s}}}{\bar{r}}\right)^{4}$.

$\bullet$ The \emph{conformation}~\eqref{eq:conformation-in-spacetime} reduces to
\begin{equation}\label{eq:conformation-s}
  \mathbf{H}=\mathbf{F}^{-1}\left( \frac{1}{k}\bq^{-1}-\frac{1}{c^{2}} \bm{u} \otimes \bm{u}\right) \mathbf{F}^{-\star}.
\end{equation}

$\bullet$ For the two choices of a \emph{three-dimensional stress} ${\bm{\sigma}}$ introduced in \autoref{sec:hyperelasticity-spacetime}, where $\rho=\gamma\rho_{r}$ is the relativistic mass density, we get, in the coordinate system $(t, x^{i})$,
\begin{equation}\label{eq:Ttxi}
  \mathbf{T} =\left(\begin{array}{cc}\displaystyle
      \frac{1}{c^{2}\mathcal{N}^{2}} E_{\text{tot}} & \displaystyle\frac{1}{c^{2} \mathcal{N}} \bm{p}^{\star}
      \\
      \displaystyle\frac{1}{c^{2} \mathcal{N}} \bm{p}  & \bs
    \end{array}\right)
\end{equation}

\begin{enumerate}
  \item When the stress tensor ${\bm{\sigma}}$ is defined as the spatial part of $\mathbf{S}$:
        \begin{equation}\label{eq:T1-s}
          \begin{cases}
            E_{\mathrm{tot}}=  \gamma \rho c^{2} + E\Big(1+k \frac{\bm{u}^{2}}{c^{2}} \Big) - \frac{1}{c^{2}} \bm{u}^{\flat} \cdot {\bm{\sigma}} \cdot \bm{u}^{\flat},
            \\
            \bm{p} = (\gamma \rho c^{2}+E)\bm{u}- {\bm{\sigma}} \cdot \bm{u}^{\flat},
            \\
            \mathbf{t} = \gamma\rho \, \bm{u} \otimes \bm{u}- {\bm{\sigma}},
            \\
            {\bm{\sigma}} = - \frac{2}{\gamma k^{2}} \rho\, q^{-1} \mathbf{F}^{-\star}\frac{\partial  e}{\partial \mathbf{H}} \mathbf{F}^{-1} q^{-1}
            - \frac{\gamma^{2} E}{c^{2}} \bm{u} \otimes \bm{u}.
          \end{cases}
        \end{equation}

  \item When the stress tensor ${\bm{\sigma}}$ is defined as the spatial part of ${\bm{\Sigma}}$:
        \begin{equation}\label{eq:T2-s}
          \begin{cases}
            E_{\mathrm{tot}}=  \gamma \rho c^{2} + E\Big(1+k \gamma^{2} \frac{\bm{u}^{2}}{c^{2}}\Big)-\frac{1}{c^{2}} \bm{u}^{\flat} \cdot {\bm{\sigma}} \cdot \bm{u}^{\flat}
            ,
            \\
            \bm{p} = \left(\gamma \rho c^{2} + E\Big(1+k \gamma^{2}\frac{\bm{u}^{2}}{c^{2}} \Big)\right) \bm{u} - {\bm{\sigma}} \cdot\bm{u}^{\flat},
            \\
            \mathbf{t} = \left(\gamma \rho + \frac{1}{c^{2}}E \Big(1+k\gamma^{2} \frac{\bm{u}^{2}}{c^{2}} \Big)\right) \bm{u} \otimes \bm{u}- {\bm{\sigma}},
            \\
            {\bm{\sigma}} = - \frac{2}{\gamma k^{2}} \rho\, q^{-1} \mathbf{F}^{-\star}\frac{\partial  e}{\partial \mathbf{H}} \mathbf{F}^{-1} q^{-1}.
          \end{cases}
        \end{equation}
\end{enumerate}

\begin{rem}
  The above expressions generalize Souriau's 1958--1964 results when gravity is taken into account. Note that in~\eqref{eq:T1-s}, in the special case $k=1$ (Special Relativity), there is a contribution $E\bm{u}^{2}/c^{2}$ in the expression of $E_{\mathrm{tot}}$. This term is however missing in~\cite[p. 153]{Sou1958} but this typo was corrected in~\cite[p. 376]{Sou1964}.
\end{rem}

The full description of Relativistic Hyperelasticity must be completed by writing balance laws for the four-momentum quadrivector and the stress-energy tensor. They are written as
\begin{align}
  \label{eq:divP=0}
  \dive^{g} \mathbf{P} & = 0,
  \\
  \label{eq:divT=0}
  \dive^{g} \mathbf{T} & = 0,
\end{align}
where $\mathbf{P}$ is the current of matter and $\mathbf{T}$ is the stress-energy tensor. Using explicit formulas for these divergence operators provided in \autoref{sec:divergences-spacetime}, where we use here the Cartesian isotropic coordinates system $(t, x^{i})$, we obtain the following equations.

$\bullet$ The conservation law~\eqref{eq:divP=0} for the \emph{current of matter}
\begin{equation}\label{eq:P-s}
  \mathbf{P} = \rho_{r} \mathbf{U} = \frac{\rho}{c\mathcal{N}} \partial_{t} + \frac{1}{c}\rho\bm{u},
\end{equation}
leads (after multiplication by $c$) and according to~\eqref{eq:divP-static}, to
\begin{equation}\label{eq:Schwarzschild-divP}
  \dive^{g} c\mathbf{P} = \frac{1}{\mathcal{N}} \frac{\partial \rho}{\partial t}+ \frac{\partial}{\partial x^{i}} (\rho u^{i})+  (\Gamma^{j}_{j i } +\Gamma^{t}_{t i } ) \rho u^{i}=0,
\end{equation}
or, in a more intrinsic form
\begin{equation*}
  \dive^{g} c\mathbf{P} = \frac{1}{\mathcal{N}} \frac{\partial \rho}{\partial t} + \dive^{{g^{3D}}} (\rho \bm{u}) + \rho \bm{u}\cdot \mathrm{d} \ln \mathcal{N} = 0.
\end{equation*}

$\bullet$ The conservation law~\eqref{eq:divT=0} for the \emph{stress-energy tensor} $\mathbf{T}$, expressed in components by~\eqref{eq:Ttxi},
has for time component, according to~\eqref{eq:divT0-static}
\begin{equation}\label{eq:Schwarzschild-divT-t}
  (\dive^{g} \mathbf{T})^{t} =  \frac{1}{c^{2}\mathcal{N}^{2}}\frac{\partial E_{\mathrm{tot}}}{\partial t} + \frac{1}{c^{2}\mathcal{N}} \frac{\partial p^{i}}{\partial x^{i}}
  + \left(  \Gamma^{j}_{j i} + 3\Gamma^{t}_{ti} - \frac{\partial \mathcal{\ln N}}{\partial x^{i}}\right) \frac{p^{i}}{c^{2}\mathcal{N}} = 0,
\end{equation}
or, in a more intrinsic form
\begin{equation*}
  (\dive^{g} \mathbf{T})^{t} = \frac{1}{c^{2}\mathcal{N}^{2}}\frac{\partial E_{\mathrm{tot}}}{\partial t} + \frac{1}{c^{2}\mathcal{N}} \dive^{{g^{3D}}} \bm{p}+ 2 \frac{\bm{p}}{c^{2}\mathcal{N}} \cdot \mathrm{d} \ln \mathcal{N} =0,
\end{equation*}
and for spatial component, according to~\eqref{eq:divTi-static}
\begin{equation}  \label{eq:Schwarzschild-divT-i}
  (\dive^{g} \mathbf{T})^{i} =  \frac{1}{c^{2}\mathcal{N}} \frac{\partial p^{i}}{\partial t} + \frac{\partial s^{i j}}{\partial x^{j}}
  + (\Gamma^{k}_{k j}+\Gamma^{0}_{0 j}) s^{ij}+ \Gamma^{i}_{jk}  s^{jk}+ \Gamma^{i}_{tt} \frac{E_{\mathrm{tot}}}{c^{2}\mathcal{N}^{2}} =0,
\end{equation}
or, in a more intrinsic form
\begin{equation*}
  (\dive^{g} \mathbf{T})^{\top} =  \frac{1}{c^{2}\mathcal{N}} \frac{\partial \bm{p}}{\partial t} + \dive^{{g^{3D}}} \bs
  +   \bs \cdot \mathrm{d} \ln \mathcal{N}
  + E_{\mathrm{tot}} \grad^{{g^{3D}}}\!\! \ln \mathcal{N} =0,
\end{equation*}
where $ \grad^{{g^{3D}}}\!\! \ln \mathcal{N}= \frac{1}{k} q^{-1} \mathrm{d} \ln \mathcal{N}$.

\section{The Galilean limit of Relativistic Hyperelasticity}
\label{sec:Gallilean}

Soon after Einstein's formulation of the theory of General Relativity (1915), Cartan introduced, in a series of papers~\cite{Car1923,Car1924,Car1925}, a general covariant formulation of Newtonian gravity, today called \emph{Newton-Cartan theory of gravitation}. It can be obtained as the limit of Lorentzian spacetimes whose light cones open up to hyperplanes at each tangent space~\cite{Kuen1976}. It allows to recast the equations of energy and momentum balance of Classical Continuum Mechanics in a four-dimensional general covariant form, similar to the relativistic equation $\dive^{g} \mathbf{T} = 0$, provided the specific internal energy and the energy flux are suitably interpreted. This 4D general covariant formalism~\cite{Tou1958,Kuen1972,Dix1975,DK1977,Hav1964,DBKP1985,dSax2020}, derived from General Relativity, is also useful to better understand the foundations of Classical Continuum Mechanics and avoid \emph{ad hoc} assumptions in its formulation.

A \emph{Galilean structure} on a four-dimensional manifold $\mM$ is a pair $(\got,\theta)$, where $\got$ is a symmetric second-order contravariant tensor of signature $(0,+,+,+)$ (the classical spatial metric) and $\theta$ is a one-form which spans the kernel of $\got$ (the clock). This means that $\got \theta = 0$ and that $\theta$ vanishes nowhere. The one-form $\theta$ defines a distribution of hyperplanes $E_{m} := \ker \theta_{m}$, each of them, carrying an Euclidean metric induced by $\got$. The Galilean structure is called integrable if $\theta$ is closed. In that case, $\theta$ defines a time function $\hat{t}$ satisfying $\mathrm{d}\hat{t} = \theta$ (at least locally) and a foliation by hypersurfaces, $\espace_{t} := \hat{t}^{-1}(t)$, tangent to the distribution $(E_{m})$, and which are moreover Riemannian manifolds.

A covariant derivative $\nabla$ on $\mM$ is said to be Galilean if it is symmetric (torsion-less) and satisfies moreover
\begin{equation*}
  \nabla \got = 0 \quad \text{and} \quad \nabla \theta = 0.
\end{equation*}
Such a covariant derivative exists only if the Galilean structure $(\got,\theta)$ is integrable (since $\nabla \theta = 0$ implies $\mathrm{d}\theta = 0$ for a symmetric covariant derivative). Note, however, that contrary to the canonical covariant derivative of a Riemannian (or pseudo-Riemannian) manifold, a Galilean covariant derivative is not uniquely defined.

In practice, a Galilean structure $(\got,\theta)$ is obtained as a limit of a one-parameter family of smooth Lorentz metrics ${\overset{\scriptscriptstyle\lambda}{g}}{}$, such that
\begin{equation*}
  {\overset{\scriptscriptstyle\lambda}{g}}{}^{-1} = \got + \lambda \bm{\kappa} + O(\lambda^{2}),
\end{equation*}
with $\got$ of signature $(0,+,+,+)$, and $\theta$ is a generator of the kernel of $\got$~\cite{Dau1964,Kuen1972}. Note that $\theta$ can be fixed uniquely up to a sign by the normalization $\bm{\kappa}(\theta,\theta) = -1$. It has been shown in~\cite{Kuen1976} that, provided that $\mathrm{d}\theta = 0$ (a condition which is always satisfied by static spacetimes such as the Minkowski or the Schwarzschild spacetimes), the one-parameter family of Riemannian covariant derivatives $\nabla^{\lambda}$ converges then to a symmetric covariant derivative $\nabla^{\text{\tiny NC}}$ which is compatible with the Galilean structure $(\got,\theta)$.

Recall moreover, that, on a Riemannian or pseudo-Riemannian manifold, the Riemann tensor $\mathbf{R}$ (defined in components by $R_{\alpha\beta\gamma\delta} = g_{\alpha\rho}{R^{\rho}}_{\beta\gamma\delta}$) has the additional symmetry
\begin{equation*}
  R_{\alpha\beta\gamma\delta} = R_{\gamma\delta\alpha\beta},
\end{equation*}
and that, uprising the first and third indices, we obtain the following identity
\begin{equation*}
  g^{\gamma\lambda}{R^{\alpha}}_{\beta\lambda\delta} = g^{\alpha\lambda}{R^{\gamma}}_{\delta\lambda\beta}.
\end{equation*}
Therefore, when a Galilean covariant derivative is obtained as a limit of (pseudo-)Riemannian covariant derivatives, it must satisfies the additional property
\begin{equation}\label{eq:Newtonnian-connexion}
  \got^{\gamma\lambda}{R^{\alpha}}_{\beta\lambda\delta} = \got^{\alpha\lambda}{R^{\gamma}}_{\delta\lambda\beta}
\end{equation}
and is then called a \emph{Newtonian covariant derivative}.

Applying this procedure to the Schwarzschild metric~\eqref{eq:Schwarzschild-metric},
\begin{equation*}
  {\overset{\scriptscriptstyle\lambda}{g}}{} = - {\mathcal{N}}^{2} c^{2}  \mathrm{d} t^{2}+ k \, q,
\end{equation*}
with $\lambda := 1/c^{2}$, we get
\begin{equation}\label{eq:Nc4}
  {\mathcal{N}}=1-\frac{1}{c^{2}} \frac{GM}{\bar{r}}+O(1/c^{4}),
  \qquad
  k=1+\frac{2}{c^{2}} \frac{GM}{\bar{r}}+O(1/c^{4}),
\end{equation}
and
\begin{equation}\label{eq:metric-expansion}
  {\overset{\scriptscriptstyle\lambda}{g}}{}^{-1} = q^{-1} - \frac{1}{c^{2}} \left( (\partial_{t})^{2} + 2\frac{GM}{\bar{r}}q^{-1}\right) + O(1/c^{4}),
\end{equation}
where $q^{-1} = (\partial_{x^{1}})^{2} + (\partial_{x^{2}})^{2} + (\partial_{x^{3}})^{2}$. Observe also that we have,
\begin{equation}\label{eq:volg-expansion}
  \vol_{{\overset{\scriptscriptstyle\lambda}{g}}{}} = c f\, \mathrm{d} t \wedge \vol_{q},
  \quad \text{where} \quad  f={\mathcal{N}} k^\frac{3}{2} =1+\frac{2}{c^{2}} \frac{GM}{\bar{r}}  + O(1/c^{4}),
\end{equation}
for the Riemannian volume form associated with the metric ${\overset{\scriptscriptstyle\lambda}{g}}{}$. Note that $\vol_{{\overset{\scriptscriptstyle\lambda}{g}}{}}$ diverges as $c \to \infty$.

We obtain thus the following Galilean structure,
\begin{equation*}
  \got = q^{-1}, \qquad \bm{\kappa} = - (\partial_{t})^{2} - 2\frac{GM}{\bar{r}}q^{-1}, \qquad \theta = \mathrm{d} t,
\end{equation*}
where the normalisation condition $\bm{\kappa}(\theta,\theta) = -1$ has been used. This structure is of course integrable, and the time function is the same, in either the relativistic context or the Galilean one. Therefore, the foliation by the hypersurfaces $\Omega_{t}$ is common to both structures. Note however that by~\eqref{eq:Schwarzschild-normal}, the relativistic normal $\mathbf{N}^{\lambda}$ to these hypersurfaces converges towards $0$ as $c \to \infty$.

An immediate consequence of~\eqref{eq:metric-expansion} is the fact that the \emph{conformation} ${\overset{\scriptscriptstyle\lambda}{\bH}}$ defined by~\eqref{eq:conformation} has a limit when $c\to \infty$. Indeed,
\begin{equation}\label{eq:HapproxCm1}
  {\bar{\bH}} := \lim_{c\to \infty}{\overset{\scriptscriptstyle\lambda}{\bH}} = \mathbf{F}^{-1}\bq^{-1} \,\mathbf{F}^{-\star}.
\end{equation}

\begin{rem}
  Observe the similarity between the (limit) conformation ${\bar{\bH}}$ and the inverse $\mathbf{C}^{-1}$ of right Cauchy--Green tensor $\mathbf{C}:=\mathbf{F}^{\star} \bq\, \mathbf{F}$ in Classical Continuum Mechanics. However, ${\bar{\bH}}$ is not exactly $\mathbf{C}^{-1}$ because it is a function from the World tube $\mathcal{W}$ to $\Sym^{2}T\mathcal{B}$ (see remark~\ref{rem:Psim1}), while $\mathbf{C}^{-1}$ is a tensor field on $\mathcal{B}$.
\end{rem}

Concerning the Riemannian covariant derivative $\nabla^{\lambda}$ of ${\overset{\scriptscriptstyle\lambda}{g}}{}$, the expansion of its non vanishing Christoffel symbols is easily deduced from~\eqref{eq:Schwarzschild-Christoffels} and recalling that $\bar{r}_{\mathrm{s}} = GM/2c^{2}$. We get
\begin{equation}\label{eq:Schwarzschild-Christoffels-expansions}
  \begin{split}
    \Gamma^{i}_{tt} & =  2c^{2} \frac{\bar{r}_{\mathrm{s}}}{\bar{r}^{2}} \frac{(1-\bar{r}_{\mathrm{s}}/\bar{r})}{(1+\bar{r}_{\mathrm{s}}/\bar{r})^{7}} \frac{x^{i}}{\bar{r}} = -\mathrm{g}^{i} -\frac{1}{c^{2}}\frac{4GM}{\bar{r}^{2}}\mathrm{g}^{i} + O(1/c^{4}),
    \\
    \Gamma^{t}_{ti} & =  2 \frac{\bar{r}_{\mathrm{s}}}{\bar{r}^{2}} \frac{ 1}{\left(1-\bar{r}_{\mathrm{s}}^{2}/\bar{r}^{2}\right)} \delta_{ik} \frac{x^{k}}{\bar{r}} = -\frac{1}{c^{2}}\delta_{ik}\mathrm{g}^{k}+O(1/c^{4}),
    \\
    \Gamma^{i}_{jj} & = 2 \frac{\bar{r}_{\mathrm{s}}}{\bar{r}^{2}} \frac{ 1}{\left(1+\bar{r}_{\mathrm{s}}/\bar{r}\right)} \frac{x^{i}}{\bar{r}} = - \frac{1}{c^{2}}\mathrm{g}^{i}+O(1/c^{4}), \quad \text{for $i\neq j$},
    \\
    \Gamma^{j}_{ji} & = - 2 \frac{\bar{r}_{\mathrm{s}}}{\bar{r}^{2}} \frac{ 1}{\left(1+\bar{r}_{\mathrm{s}}/\bar{r}\right)}\delta_{ik}\frac{x^{k}}{\bar{r}} = \frac{1}{c^{2}}\delta_{ik}\mathrm{g}^{k}+O(1/c^{4}), \quad \text{for $i = j$ and $i\ne j$}.
  \end{split}
\end{equation}
with no sum on $j$, and where $\mathbf{g}$ is the Newtonian (centripetal) gravity field,
\begin{equation}\label{eq:gxi}
  \mathbf{g} := -\frac{GM} {\bar{r}^{2}} \frac{x^{i}}{\bar{r}} \partial_{x^{i}},
\end{equation}
with $\frac{GM} {\bar{r}^{2}}\approx \frac{GM} {r_{0}^{2}} = 9.8\, \text{m/s}^{2}$ on Earth surface. We deduce therefore that the Christoffel symbols of the Newton--Cartan limit $\nabla^{\text{\tiny NC}}$ are all vanishing except
\begin{equation}\label{eq:lmitGammaSchw}
  \Gamma^{i}_{tt} = -\mathrm{g}^{i}.
\end{equation}

\begin{rem}[Weak gravity]\label{rem:gravity}
  The approximations~\eqref{eq:Schwarzschild-Christoffels-expansions} at large $c$ coincide in fact with those corresponding to so-called \emph{weak gravity} (see \cite{Der2006}), valid for the Earth, and more generally for a planet or a star (any object whose radius is smaller than its Schwarzschild radius being called a black hole).
\end{rem}

The divergence of a quadrivector $\mathbf{P}$, relative to the Newtonian covariant derivative $\nabla^{\text{\tiny NC}}$, is given by
\begin{equation}\label{eq:P-Galilean-divergence}
  \divz \mathbf{P} = \frac{\partial P^{t}}{\partial t} + \frac{\partial P^{j}}{\partial x^{j}}.
\end{equation}
It corresponds to the zero order terms in the expansions in $\lambda=1/c^{2}$ of the divergence $\divl \mathbf{P}$ relative to the Schwarzschild metric ${\overset{\scriptscriptstyle\lambda}{g}}{}$ given by~\eqref{eq:divP-static}, since we have
\begin{equation*}
  \Gamma^{j}_{j i } +\Gamma^{t}_{t i }=\frac{2}{c^{2}} \delta_{ik} \mathrm{g}^{k}+ O(1/c^{4}).
\end{equation*}
We get therefore
\begin{equation}\label{eq:divPc4}
  \divl \mathbf{P}=\divz \mathbf{P}+\frac{2}{c^{2}} \mathbf{g}\cdot \mathbf{P}^{\top}+O(1/c^{4}),
\end{equation}
where $\mathbf{P}^{\top}$ is the spatial part of the quadrivector $\mathbf{P}$.

The divergence of a symmetric second-order contravariant tensor $\mathbf{T}$ is given by
\begin{equation}\label{eq:T-Galilean-divergence}
  (\divz \mathbf{T})^{t} =  \frac{\partial T^{tt}}{\partial t} + \frac{\partial T^{tj}}{\partial x^{j}},
  \qquad (\divz \mathbf{T})^{i} =  \frac{\partial T^{it}}{\partial t} +  \frac{\partial T^{ij}}{\partial x^{j}} - \mathrm{g}^{i}T^{tt}.
\end{equation}
Indeed,  by~\eqref{eq:divT0-static}, and since by~\eqref{eq:Schwarzschild-Christoffels},
\begin{equation}\label{eq:sumGammai}
  \sum_{j=1}^{3}\Gamma^{j}_{ji}+3 \Gamma^{t}_{ti} =3\left(\mathrm{d} \ln(\mathcal{N}\sqrt{k})\right)_{i}= O(1/c^{4}),
\end{equation}
we obtain
\begin{equation}\label{eq:divTtc4}
  (\divl \mathbf{T})^{t} = (\divz \mathbf{T})^{t}+O(1/c^{4}),
\end{equation}
whereas~\eqref{eq:divTi-static-components}, combined with the fact that all involved Christoffel's symbols are $O(1/c^{2})$, but
\begin{equation*}
  \Gamma^{i}_{tt} = -\mathrm{g}^{i} + O(1/c^{2}),
\end{equation*}
ends up to
\begin{equation*}
  (\divl \mathbf{T})^{i} = \frac{\partial T^{it}}{\partial t} + \frac{\partial T^{ij}}{\partial x^{j}} - \mathrm{g}^{i}T^{tt}+ O(1/c^{2}) = (\divz \mathbf{T})^{i} + O(1/c^{2}).
\end{equation*}

The \emph{current of matter}~\eqref{eq:def-P} ${\overset{\scriptscriptstyle\lambda}{\bP}}$ for the metric ${\overset{\scriptscriptstyle\lambda}{g}}{}$ is defined implicitly by
\begin{equation*}
  i_{\mathbf{P}^{\lambda}} \vol_{{\overset{\scriptscriptstyle\lambda}{g}}{}}=\omega=\Psi^{*} \mu,
\end{equation*}
where the 3-form $\omega$ does not depend on the light speed $c$ (the matter field $\Psi$ and the mass measure $\mu$ do not depend on $c$, which is only introduced through the metrics). By~\eqref{eq:volg-expansion}, we get thus
\begin{equation*}
  cf\,  i_{\mathbf{P}^{\lambda}} \,\mathrm{d} t\wedge \vol_{q}= i_{cf \mathbf{P}^{\lambda}} \,\mathrm{d} t\wedge \vol_{q}=\omega=\Psi^{*} \mu,
\end{equation*}
from which it is seen that $cf {\overset{\scriptscriptstyle\lambda}{\bP}}$ is independent of $c$, and is therefore equal to its Newtonian limit $(c\mathbf{P})^{\scriptscriptstyle 0}$ defined by $i_{(c\mathbf{P})^{\scriptscriptstyle 0}} \,\mathrm{d} t\wedge \vol_{q}=\omega$ (since $f \to 1$ as $c\to \infty$).
Setting
\begin{equation}\label{eq:cPLcP0}
  c{\overset{\scriptscriptstyle\lambda}{\bP}} = {\overset{\scriptscriptstyle\lambda}{\rho}} \,\partial_{t} + {\overset{\scriptscriptstyle\lambda}{\rho}}\, {\overset{\scriptscriptstyle\lambda}{\uu}}{},
  \qquad
  (c\mathbf{P})^{\scriptscriptstyle 0} = {\bar{\rho}}\, \partial_{t} + {\bar{\rho}}\, {\bar{\uu}},
\end{equation}
defines the mass density ${\overset{\scriptscriptstyle\lambda}{\rho}}$ and the velocity ${\overset{\scriptscriptstyle\lambda}{\uu}}{}$, as well as their Newtonian limits ${\bar{\rho}}$ and ${\bar{\uu}}$. The equality $cf \mathbf{P}^{\lambda}=(c\mathbf{P})^{\scriptscriptstyle 0}$ leads to
\begin{equation*}
  {\overset{\scriptscriptstyle\lambda}{\rho}} = \frac{1}{f} {\bar{\rho}},\qquad {\overset{\scriptscriptstyle\lambda}{\uu}}{} = {\bar{\uu}}.
\end{equation*}
and allows us to expand the mass density as
\begin{equation*}
  {\overset{\scriptscriptstyle\lambda}{\rho}} = {\bar{\rho}}- \frac{2}{c^{2}} \frac{GM}{\bar{r}} {\bar{\rho}} + O(1/c^{4}).
\end{equation*}
We have therefore
\begin{equation*}
  \lim_{c\to\infty} {\overset{\scriptscriptstyle\lambda}{\rho}} = {\bar{\rho}},
  \qquad
  \lim_{c\to\infty} {\overset{\scriptscriptstyle\lambda}{\uu}}{} = {\bar{\uu}},
\end{equation*}
and
\begin{equation*}
  \lim_{c\to\infty} \gamma=\lim_{c\to\infty} \frac{1}{\sqrt{1-k \frac{\bm{u}^2}{c^{2}}}}=1,
\end{equation*}
since the \emph{generalized Lorentz factor}~\eqref{eq:gamma-s} has the classical expansion
\begin{equation*}
  \gamma = 1+\frac{1}{2} \frac{{\bar{\uu}}^2}{c^{2}} + O(1/c^{4}).
\end{equation*}

Using~\eqref{eq:P-s}, we see that $({\overset{\scriptscriptstyle\lambda}{\rho}},{\overset{\scriptscriptstyle\lambda}{\uu}}{})$ are connected to $(\rho,\bm{u})$ by
\begin{equation*}
  {\overset{\scriptscriptstyle\lambda}{\rho}}=\frac{\rho}{{\mathcal{N}}},
  \qquad
  {\overset{\scriptscriptstyle\lambda}{\uu}}{}={\mathcal{N}} \bm{u} ,
  \qquad
  {\overset{\scriptscriptstyle\lambda}{\rho}}\,{\overset{\scriptscriptstyle\lambda}{\uu}}{}=\rho \bm{u},
\end{equation*}
when the so-called relativistic mass density $\rho$ and velocity $\bm{u}$ are defined by the orthogonal decomposition
\begin{equation*}
  c{\overset{\scriptscriptstyle\lambda}{\bP}}=\rho c\, {\overset{\scriptscriptstyle\lambda}{\mathbf{N}}}+\rho \bm{u}.
\end{equation*}
By~\eqref{eq:u-s}, we deduce that the Newtonian limit of the three-dimensional velocity $\bm{u}={\overset{\scriptscriptstyle\lambda}{\uu}}{}/{\mathcal{N}}$ is
\begin{equation*}
  \lim_{c\to \infty} \bm{u} ={\bar{\uu}} =  - \mathbf{F} \frac{\partial \Psi}{\partial t}.
\end{equation*}

\begin{rem}
  In Classical Continuum Mechanics, the \emph{Eulerian velocity} is defined as the vector field on the deformed configuration given by
  $\partial_{t} p \circ p^{-1}$, where $p$ is the embedding of the body $\mathcal{B}$ into the Euclidean space, and where $\VV:=\partial_{t} p$ is the \emph{Lagrangian velocity}. If we assume, furthermore, that $\Psi_{t}$ is a diffeomorphism and we set $p=\Psi_{t}^{-1}$ (see remark \autoref{rem:TPsit-iso}), the vector field ${\bar{\uu}}$ recasts as
  \begin{equation*}
    {\bar{\uu}} = - \mathbf{F}\,\partial_{t} \Psi = \partial_{t} \Psi_{t}^{-1} \circ \Psi_{t} =
    \partial_{t} p \circ p^{-1},
  \end{equation*}
  and we recognize ${\bar{\uu}}$ as the Eulerian velocity of Classical Continuum Mechanics.
\end{rem}

The \emph{stress-energy tensor} ${\overset{\scriptscriptstyle\lambda}{\bT}}:=\mathbf{T}$ has for expression, in the coordinate system $(t, x^{i})$,
\begin{equation}\label{eq:TtxiN}
  {\overset{\scriptscriptstyle\lambda}{\bT}} =\begin{pmatrix}
    \displaystyle \frac{1}{c^{2} {\mathcal{N}}} {\overset{\scriptscriptstyle\lambda}{E}}_{\text{tot}} & \displaystyle \frac{1}{c^{2} {\mathcal{N}}} {\overset{\scriptscriptstyle\lambda}{\bp}}{}^{\star}
    \\
    \displaystyle\frac{1}{c^{2} {\mathcal{N}}} {\overset{\scriptscriptstyle\lambda}{\bp}}{}               & \displaystyle {\overset{\scriptscriptstyle\lambda}{\bs}}
  \end{pmatrix}
  \qquad\text{where}\qquad
  \begin{cases}
    {\overset{\scriptscriptstyle\lambda}{E}}_{\text{tot}}:=E_{\text{tot}}/\mathcal{N},
    \\
    {\overset{\scriptscriptstyle\lambda}{\bp}}{}:= \bm{p},
    \\
    {\overset{\scriptscriptstyle\lambda}{\bs}}:=\bs.
  \end{cases}
\end{equation}
To determine its limit, observe that the energy density ${\overset{\scriptscriptstyle\lambda}{E}}_{\text{tot}}$, the linear momentum ${\overset{\scriptscriptstyle\lambda}{\bp}}{}$ and the spatial part ${\overset{\scriptscriptstyle\lambda}{\bs}}$ behave as
\begin{align}\label{eq:Etotlambda}
  \frac{{\overset{\scriptscriptstyle\lambda}{E}}_{\text{tot}}}{c^{2}} &
  = {\overset{\scriptscriptstyle\lambda}{\rho}}+ \frac{1}{c^{2}}\left(E+\frac{1}{2} \rho \bm{u}^{2}\right) +  O(1/c^{4}),
  \\
  \label{eq:plambda}
  \frac{{\overset{\scriptscriptstyle\lambda}{\bp}}{}}{c^{2}}              & =   {\overset{\scriptscriptstyle\lambda}{\rho}} \,{\overset{\scriptscriptstyle\lambda}{\uu}}{} + \frac{1}{c^{2}} \left(\left(E+\frac{1}{2} \rho \bm{u}^{2}\right) \bm{u}-{\bm{\sigma}}\cdot\bm{u}^{\flat}\right) + O(1/c^{4}),
  \\
  \label{eq:slambda}
  {\overset{\scriptscriptstyle\lambda}{\bs}}                           & = {\overset{\scriptscriptstyle\lambda}{\rho}}\,  {\overset{\scriptscriptstyle\lambda}{\uu}}{}   \otimes {\overset{\scriptscriptstyle\lambda}{\uu}}{}   -{\bm{\sigma}}
  + O(1/c^{2}),
\end{align}
where we have used indifferently either~\eqref{eq:T1-s} or~\eqref{eq:T2-s}, and we have assumed that the internal energy density $E$ (function of $\mathbf{H}$) is $o(1/c^{2})$ (see also the discussion in \cite{Sou1960}). In that case, the quantities $\frac{1}{2}\rho \bm{u}^{2}$, $E$ and ${\bm{\sigma}}$ converge respectively to $\frac{1}{2}{\bar{\rho}}\, {\bar{\uu}}^{2}$, ${\bar{E}}$ and ${\bar{\bsigma}}$.
Therefore, in the coordinate system $(t, x^{i})$, the stress-energy tensor ${\overset{\scriptscriptstyle\lambda}{\bT}}$ converges  to the Newtonian limit
\begin{equation}\label{eq:TGal}
  {\bar{\bT}} = \lim_{c\to \infty}\mathbf{T}^{\lambda}=
  \begin{pmatrix}
    {\bar{\rho}}        & {\bar{\rho}}\, {\bar{\uu}}^{\star}              \\
    {\bar{\rho}}\, {\bar{\uu}} & {\bar{\rho}}\, {\bar{\uu}} \otimes {\bar{\uu}}- {\bar{\bsigma}}
  \end{pmatrix}.
\end{equation}

We now discuss the limits of the balance laws. First, by~\eqref{eq:divPc4} and~\eqref{eq:cPLcP0}, we get
\begin{equation*}
  \divl (c{\overset{\scriptscriptstyle\lambda}{\bP}}) = \frac{\partial {\overset{\scriptscriptstyle\lambda}{\rho}}}{\partial t}+\dive({\overset{\scriptscriptstyle\lambda}{\rho}}\, {\overset{\scriptscriptstyle\lambda}{\uu}}{})+\frac{2}{c^{2}} \mathbf{g}\cdot {\overset{\scriptscriptstyle\lambda}{\rho}}\, {\overset{\scriptscriptstyle\lambda}{\uu}}{}+O(1/c^{4})=0,
\end{equation*}
where $\dive$ is the canonical divergence in $\RR^{3}$, and which can be recast as
\begin{equation}\label{eq:divrholul}
  \frac{\partial {\overset{\scriptscriptstyle\lambda}{\rho}}}{\partial t}+\dive({\overset{\scriptscriptstyle\lambda}{\rho}}\, {\overset{\scriptscriptstyle\lambda}{\uu}}{})=-\frac{2}{c^{2}} \mathbf{g}\cdot \rho {\overset{\scriptscriptstyle\lambda}{\uu}}{}+O(1/c^{4}).
\end{equation}
It converges to
\begin{equation*}
  \divz (c\mathbf{P})^{\scriptscriptstyle 0} = \frac{\partial {\bar{\rho}}}{\partial t} +\dive({\bar{\rho}}\, {\bar{\uu}})
  = 0.
\end{equation*}
One recovers thus the usual expression of mass conservation in Classical Continuum Mechanics (omitting the bars),
\begin{equation}\label{eq:massconservationMMC}
  \frac{\partial \rho}{\partial t} + \dive(\rho \bm{u})=0,
\end{equation}
for the Euclidean metric $q$.

Then, the relativistic conservation law $\divl \mathbf{T}^{\lambda} = 0$ converges to the equation $\divz {\bar{\bT}}=0$, with (by~\eqref{eq:T-Galilean-divergence})
\begin{equation*}
  (\divz {\bar{\bT}})^{t} = \frac{\partial {\bar{\rho}}}{\partial t} + \frac{\partial}{\partial x^{j}}( {\bar{\rho}}\, {\bar{u}}^{j}),
  \quad \text{and} \quad
  (\divz {\bar{\bT}})^{i} =  \frac{\partial}{\partial t}({\bar{\rho}}\, {\bar{u}}^{i}) +  \frac{\partial}{\partial x^{j}}( {\bar{\rho}}\, {\bar{u}}^{i} {\bar{u}}^{j} - {\bar{\sigma}}^{ij}) - {\bar{\rho}}\, \mathrm{g}^{i},
\end{equation*}
so that
\begin{equation*}
  \frac{\partial {\bar{\rho}}}{\partial t} +\dive({\bar{\rho}}\, {\bar{\uu}})=0,
  \qquad
  \frac{\partial }{\partial t} ({\bar{\rho}}\,{\bar{\uu}})+\dive\left({\bar{\rho}}\, {\bar{\uu}} \otimes {\bar{\uu}}-{\bar{\bsigma}}\right)-{\bar{\rho}} \mathbf{g}=0.
\end{equation*}

The first equation is (again) recognized as the mass conservation~\eqref{eq:massconservationMMC} and the second one as the linear momentum balance of Classical Continuum Mechanics, with gravity $\mathbf{g}$ (omitting the bars),
\begin{equation*}
  \frac{\partial }{\partial t} (\rho\,\bm{u})+\dive\left(\rho\, \bm{u} \otimes \bm{u}-{\bm{\sigma}}\right)-\rho \mathbf{g}=0.
\end{equation*}
By using the mass conservation law, the later can be recast as the classical expression,
\begin{equation}\label{eq:LinearMomentumCCM}
  \rho \left( \frac{\partial \bm{u}}{\partial t}+ \nabla_{\bm{u}} \bm{u} \right)= \dive  {\bm{\sigma}} + \rho \mathbf{g},
\end{equation}
where $\nabla$ is the covariant derivative for the Euclidean metric $q$.

It is furthermore possible to recover the so-called \emph{local form of energy balance} of Classical Continuum Mechanics, as a term of order $O(1/c^{2})$ in the expansion of a combination of both $\divl (c{\overset{\scriptscriptstyle\lambda}{\bP}})=0$ and $(\divl \mathbf{T}^{\lambda})^{t} = 0$ (see for instance~\cite{Sou1958} for the case of the flat Minkowski spacetime or~\cite{Gou2012} for relativistic fluids in the case of weak gravity or \cite{Kuen1976} for general discussions about this balance law).

The balance law~\eqref{eq:Schwarzschild-divT-t} expresses the vanishing of the time component $(\dive^{g} \mathbf{T})^{t}=(\divl {\overset{\scriptscriptstyle\lambda}{\bT}})^{t}=0$  in which here $\mathbf{T}={\overset{\scriptscriptstyle\lambda}{\bT}}$ is given by~\eqref{eq:TtxiN}.
Since ${\overset{\scriptscriptstyle\lambda}{E}}_{\mathrm{tot}}=E_{\mathrm{tot}}/{\mathcal{N}}$ and ${\overset{\scriptscriptstyle\lambda}{\bp}}{}=\bm{p}$, it recasts as
\begin{equation*}
  (\dive^{g} \mathbf{T})^{t} = (\divl {\overset{\scriptscriptstyle\lambda}{\bT}})^{t} = \frac{1}{{\mathcal{N}}}\left[\frac{\partial }{\partial t} \Big( \frac{{\overset{\scriptscriptstyle\lambda}{E}}_{\mathrm{tot}}}{c^{2}}\Big) + \dive \Big( \frac{{\overset{\scriptscriptstyle\lambda}{\bp}}{}}{c^{2}}\Big) + \frac{1}{c^{2}}{\overset{\scriptscriptstyle\lambda}{\bp}}{}\cdot  \left(3\,\mathrm{d} \ln(\mathcal{N}\sqrt{k}) - \mathrm{d} \,\mathcal{\ln N}\right)\right] = 0,
\end{equation*}
where $\dive$ is the divergence relative to the three-dimensional Euclidean metric $q$. We have introduced the 1-form~\eqref{eq:sumGammai}
\begin{equation*}
  3\,\mathrm{d} \ln(\mathcal{N}\sqrt{k}) := \Big(\sum_{j=1}^{3} \Gamma^{j}_{j i} +3\Gamma^{t}_{ti}\Big)  \mathrm{d} x^{i}=O(1/c^{4}),
\end{equation*}
which is of order $O(1/c^{4})$,  ${\overset{\scriptscriptstyle\lambda}{\bp}}{}$ which is of order $O(c^{2})$ and
\begin{equation*}
  \mathrm{d} \,\mathcal{\ln N}=-\frac{1}{c^{2}} q\mathbf{g}+O(1/c^{4}).
\end{equation*}
By~\eqref{eq:Nc4}, we get
\begin{equation}\label{eq:divTlc4}
  {\mathcal{N}}(\divl {\overset{\scriptscriptstyle\lambda}{\bT}})^{t} = \frac{\partial }{\partial t} \Big( \frac{{\overset{\scriptscriptstyle\lambda}{E}}_{\mathrm{tot}}}{c^{2}}\Big) + \dive \Big( \frac{{\overset{\scriptscriptstyle\lambda}{\bp}}{}}{c^{2}}\Big) + \frac{1}{c^{2}}\Big(\frac{{\overset{\scriptscriptstyle\lambda}{\bp}}{}}{c^{2}}\cdot \mathbf{g} \Big)+ O(1/c^{4}) = 0.
\end{equation}
Subtracting~\eqref{eq:divrholul} from~\eqref{eq:divTlc4} and using~\eqref{eq:Etotlambda} and~\eqref{eq:plambda}, we get now
\begin{multline*}
  \frac{1}{c^{2}} \left\{ \frac{\partial}{\partial t}\Big[ \Big(E+\frac{1}{2} \rho \bm{u}^{2}\Big) \bm{u}-{\bm{\sigma}}\cdot\bm{u}^{\flat} \Big]
  + \dive \Big[\Big(E+\frac{1}{2} \rho \bm{u}^{2}\Big) \bm{u}-{\bm{\sigma}}\cdot\bm{u}^{\flat}\Big] + \rho\,\mathbf{g}\cdot \bm{u} \right\}
  \\
  - \frac{2}{c^{2}}\rho\,\mathbf{g}\cdot \bm{u} + O(1/c^{4}) = 0,
\end{multline*}
and thus
\begin{equation*}
  \frac{\partial}{\partial t} \Big( E+\frac{1}{2}\rho\bm{u}^{2}\Big)
  +\dive\left(\Big(E+\frac{1}{2}\rho\bm{u}^{2}\Big) \bm{u} -{\bm{\sigma}}\cdot\bm{u}^{\flat}\right)
  -\rho\, \mathbf{g}\cdot \bm{u}=O(1/c^{2}).
\end{equation*}
Passing to the limit $c \to \infty$, we obtain therefore
\begin{equation*}
  \frac{\partial}{\partial t} \Big( {\bar{E}}+\frac{1}{2}\rho{\bar{\uu}}^{2}\Big)
  +\dive\left(\Big({\bar{E}}+\frac{1}{2}{\bar{\rho}}\,{\bar{\uu}}^{2}\Big) {\bar{\uu}} -{\bar{\bsigma}}\cdot{\bar{\uu}}^{\flat}\right)
  -{\bar{\rho}}\, \mathbf{g}\cdot {\bar{\uu}}=0.
\end{equation*}
Using~\eqref{eq:LinearMomentumCCM}, we furthermore have (omitting the bars, with $\nabla$ still the covariant derivative for the Euclidean metric $q$),
\begin{equation*}
  \dive \big({\bm{\sigma}}\cdot\bm{u}^{\flat} \big)
  =\bm{u} \cdot \dive{\bm{\sigma}} + {\bm{\sigma}}:\mathbf{d} =\frac{1}{2} \rho \left( \frac{\partial \bm{u}^{2}}{\partial t}+ \nabla_{\bm{u}} \bm{u}^{2} \right)+ {\bm{\sigma}}:\mathbf{d} - \rho \, \mathbf{g}\cdot \bm{u} ,
\end{equation*}
since ${\bm{\sigma}}$ is symmetric, and where
\begin{equation*}
  \mathbf{d}:= \frac{1}{2}\left(\nabla \bm{u}^{\flat} +(\nabla \bm{u}^{\flat})^{\star} \right),
\end{equation*}
is the classical strain rate tensor. We end up with the standard expression of \emph{internal energy balance} in Classical Continuum Mechanics \cite{MH1994,LC1985},
\begin{equation*}
  \rho  \left( \frac{\partial e}{\partial t}+ \nabla_{\bm{u}} e \right)= {\bm{\sigma}}:\mathbf{d} ,
\end{equation*}
with no heat transfer, and where $e:=E/\rho$ is the specific internal energy.

\section{Conclusion}

We have revisited Souriau's variational formulation of Relativistic Hyperelasticity. This theory was derived in 1958 with the mindset of Gauge theory: the \emph{perfect matter field} $\Psi$ is somehow similar to the wave function $\psi$ in Quantum Mechanics, but at a macroscopic scale and without the same interpretation. The role of the three-dimensional \emph{body} $\mathcal{B}$, which labels the material points constitutive of the continuous medium under study in the Universe, has been emphasized: it is common to the Relativistic Hyperelasticity theory and to the three-dimensional Classical Continuum Mechanics theory. The body $\mathcal{B}$ is naturally distinguished from a reference configuration $\Omega_{t_{0}}$ in the present Relativistic framework, since $\mathcal{B}$ is embedded into the (non-physical) vector space $V$, whereas $\Omega_{t_{0}}$ is a spacelike submanifold of the Universe $\mM$. In both the Classical and Relativistic frameworks, the body $\mathcal{B}$ is endowed with a volume form, the \emph{mass measure} $\mu$, and a \emph{fixed Riemannian metric} ${\bm{\gamma}}_{0}$. Since this is shared by both theories, we have tried to make a parallel, when possible, and to point out the differences. Our point of view is mainly oriented towards mechanics, rather than astrophysics.

The fundamental observation of Souriau is that the formulation of general covariant constitutive equations for perfect matter involve the metric $g$ only through the \emph{conformation}, defined by
\begin{equation*}
  \mathbf{H}=(T\Psi) g^{-1} (T\Psi)^{\star}.
\end{equation*}
It is a non degenerate contravariant three-dimensional tensor valued function which plays the role of strain, or more precisely of the inverse of the right Cauchy--Green tensor $\mathbf{C}$. The connections between $\mathbf{H}$, the four-dimensional degenerate metric $h=g+\mathbf{U}^{\flat}\otimes \mathbf{U}^{\flat}$ and the four-dimensional degenerate co-metric $h^{\sharp}=g^{-1}+\mathbf{U}\otimes \mathbf{U}$ (considered by Carter and Quintana \cite{CQ1972}, for instance) are given by lemma~\ref{lem:h-H},
\begin{equation*}
  \mathbf{H}=(T\Psi) h^{\sharp} (T\Psi)^{\star},
  \qquad
  h=(T\Psi)^{\star} \mathbf{H}^{-1}T\Psi.
\end{equation*}
Thanks to these formulas, all the definitions of a strain tensor can be expressed using a comparison between the inverse of the conformation $\mathbf{H}^{-1}$ and $\mathbf{H}_{0}^{-1}:={\bm{\gamma}}_{0} \circ \Psi $, where ${\bm{\gamma}}_{0}$ is a reference metric on the body $\mathcal{B}$. Among these definitions, we mention
\begin{equation*}
  \mathfrak E:= \frac{1}{2} \left(\mathbf{H}^{-1}-\mathbf{H}_{0}^{-1}\right)
  \quad \text{and} \quad
  \widehat{\mathfrak E} := -\frac{1}{2}\log \big(\mathbf{H}\, \mathbf{H}_{0}^{-1}\big).
\end{equation*}
The links between the different metrics and strain tensors encountered in the literature, either defined on the World tube $\mathcal{W}$, or on the body $\mathcal{B}$, have been clarified (in \autoref{sec:Conformation-strains}, \autoref{sec:3D-metrics} and \autoref{sec:strains}).

In the framework of Variational General Relativity, the stress-energy tensor $\mathbf{T}$ derives from a \emph{general covariant Lagrangian} (theorem~\ref{thm:relativistic-hyperelasticity} and remark~\ref{rem:T}) and its decompositions allow for the rigorous formulation of \emph{stress tensors},
\begin{itemize}
  \item first, four-dimensional, such as $\mathbf{S}$ or ${\bm{\Sigma}}$ (with a preference for Eckart--Bennoun definition~\eqref{eq:Hyperelasticity-law-Sigma}),
  \item and, then, three-dimensional, such as the \emph{generalized second Piola--Kirchhoff stress tensor} $\bs$, which is always defined by~\eqref{eq:def-s-3D}, and the \emph{generalized Cauchy stress tensor} ${\bm{\sigma}}$, which definition requires the introduction of a spacetime (with a preference for our definition~\eqref{eq:Hyperelasticity-b}).
\end{itemize}
The full Relativistic Hyperelasticity theory is four-dimensional, but its constitutive laws are essentially three-dimensional and very similar to the ones of Classical Continuum Mechanics, a feature which has been used in~\cite{EBT2006} and~\cite{Bro2021}. We have formalized it in \autoref{sec:hyperelasticity-spacetime} and in \autoref{sec:stress-on-the-body}.

By considering the Schwarzschild spacetime (instead of the flat Minkowski spacetime like Souriau did), we have been able to take into account gravity. Following, this time, Künzle~\cite{Kuen1976}, we have recovered the Newton-Cartan formulation of Hyperelasticity in Galilean Relativity, as the limit $c\to \infty$ of our  relativistic formulation in Schwarzschild spacetime.