# Instabilities and geometry of growing tissues

Doron Grossman<EMAIL_ADDRESS>College de France  Jean-
Francois Joanny  College de France

###### Abstract

We derive a complete, coarse grained, continuum model of the 2D vertex model.
It is applicable for different underlying geometries, and allows for
analytical analysis of an otherwise numerical model. Using a geometric
approach and out–of–equilibrium statistical mechanics, we calculate both
mechanical and dynamical instabilities within a tissue, and their dependence
on different variables, including activity, and disorder. Most notably, the
tissue’s response depends on the existence of mechanical residual stresses on
a cellular level. Thus, even freely growing tissues may exhibit a growth
instability depending on food consumption. Using this geometric model we can
readily distinct between elasticity and plasticity in a growing, flowing,
tissue.

Cellular media such as tissues are often described using the so-called vertex
model, in which cells are described as confluent polygons or polyhedra. In
this paper, we consider 2-dimensional tissues, which are a good representation
of squamous epithelia. Each cell is assigned a preferred area and perimeter as
sketched in Fig . 1.a [1, 2]. The energy of the tissue depends on the
difference of the actual area and perimeter from their preferred values.

$\displaystyle E$
$\displaystyle=\sum_{i}\hat{k}_{A}\left(A_{i}-\hat{A}^{0}_{i}\right)^{2}+k_{P}\left(\,{}^{2}\\!P_{i}-\,{}^{2}\\!P^{0}_{i}\right)^{2}.$
(1)

We use a slightly modified model replacing the actual perimeter by
$\,{}^{2}\\!P=\sum_{e}\ell_{e}^{2}$, where $\ell_{e}$ is the length of edge
$e$ of the cell. This model agrees with the standard formulation of the vertex
model to linear order, and we expect it to exhibit a qualitatively similar
behaviour beyond linearity, while it significantly simplifies further
calculations (see appendix A). The area and perimeter moduli are positive
properties of the cells.

In a tissue, cells divide, die, and move past each other exchanging neighbors
via $T_{1}$ transition, as shown in fig. 1). All these transformations change
the edge-network topology, and can relax internal stresses. Recent studies
show subtle solid-solid transitions, due to the appearance of soft-deformation
modes, depending both on the ratio between the reference perimeter and area
$\,{}^{2}\\!P^{0}/A^{0}$, and the magnitude of disorder. [3, 4, 5, 6].The
discrete vertex model typically requires numerical studies [7, 8, 9, 10].
Other theoretical descriptions of tissues include phenomenological continuum
models [11, 12], or otherwise neglecting relaxation and plasticity [13, 5].
One of the most significant difficulties is the distinction between solid-like
(elastic) behavior, and fluid-like behavior [10, 14]. Here we present a
rigorous derivation of a continuum vertex model using an intrinsic, geometric,
approach. The result is an easily generalized, 2-dimensional model, that takes
into account geometry and out - of - equilibrium statistical mechanics, which
exemplifies the difference between elasticity and the plasticity in the vertex
model, and allows for an analytical treatment of the various topological
processes allowed in the tissue.

Given a configuration $\vec{f}\left(x^{\mu}\right)$, describing the position
of cells vertices at coordinates $x^{\mu}$, the (induced) metric is defined by
$g_{\mu\nu}=\partial_{\mu}\vec{f}\cdot\partial_{\nu}\vec{f}$. Assuming cells’
size is the smallest scale in the system, we approximate the distance between
two neighboring vertices by $\ell_{12}^{2}\simeq g_{\mu\nu}\Delta
x_{12}^{\mu}\Delta x_{12}^{\nu}$ ($\Delta
x_{12}^{\mu}=x_{2}^{\mu}-x_{1}^{\nu}$ being the coordinate difference between
vertices). The perimeter is
$\,{}^{2}\\!P_{i}=\sum_{e}\ell_{e}^{2}=g_{\mu\nu}(x_{i})Q^{\mu\nu}_{i}$. We
define the network tensor:

$\displaystyle Q_{i}^{\mu\nu}=\sum_{e}\Delta x_{e}^{\mu}\Delta x_{e}^{\nu},$
(2)

where the sum is taken over all edges $e$ of the $i^{th}$ cell.

The cell area itself is not a simple function of $Q^{\mu\nu}_{i}$, still it
can be written as $A_{i}=s_{i}\sqrt{g(x_{i})}\sqrt{Q_{i}}$ where $g=\det
g_{\mu\nu}$, $Q=\det Q^{\mu\nu}$, and
$s_{i}(Q^{\mu\nu})=\frac{A_{i}}{\sqrt{g(x_{i})Q_{i}}}\neq 0$, is a finite
correction term that depends smoothly on $Q^{\mu\nu}$. This factor can be
factored out, renormalizing the area term in a controlled manner. The energy
may then be written as

$\displaystyle E$
$\displaystyle=\sum_{i}{k}_{A}\left(\sqrt{gQ_{i}}-{A}^{0}_{i}\right)^{2}+k_{P}\left(g_{\mu\nu}Q^{\mu\nu}_{i}-\,{}^{2}\\!P^{0}\right)^{2},$
(3)

where ${k}_{A}$ and ${A}^{0}_{i}$ are the re-normalized values.

In a continuum model [15], a tissue is described by a density function
$\Psi\left(x,A^{0},\,{}^{2}\\!P^{0},Q^{\mu\nu}\right)$, normalized so that
$\int\Psi\,\text{d}[Q,A^{0},\,{}^{2}\\!P^{0}]=\sqrt{g}\rho$, where $\rho$ is
the cell density and
$\,\text{d}[Q,A^{0},\,{}^{2}\\!P^{0}]=\,\text{d}Q^{\mu\nu}\,\text{d}A^{0}\,\text{d}\,{}^{2}\\!P^{0}$
is the integration measure over the dynamic fields. The energy of the
continuum model reads:

$\displaystyle E$
$\displaystyle=\int\int\Psi\left\\{k_{A}\left(\sqrt{gQ}-A^{0}\right)^{2}\right.$
(4)
$\displaystyle\left.+k_{P}\left(g_{\mu\nu}Q^{\mu\nu}-\,{}^{2}\\!P^{0}\right)^{2}\right\\}\,\text{d}[Q,A^{0},\,{}^{2}\\!P^{0}]\,\text{d}^{2}x$

where $\,\text{d}^{2}x=\,\text{d}x^{1}\,\text{d}x^{2}$ is the integration
measure over coordinates. Completing the integration over the dynamic fields,
we get

$\displaystyle E$
$\displaystyle=\int\rho(x)\left\\{k_{A}\left(A\left(x\right)-A^{0}\left(x\right)\right)^{2}\right.$
(5)
$\displaystyle\left.+k_{P}\left(\,{}^{2}\\!P\left(x\right)-\,{}^{2}\\!P^{0}\left(x\right)\right)^{2}+k_{A}\left(\Delta
A^{2}+\Delta{A^{0}}^{2}\right)\right.$
$\displaystyle\left.+k_{P}\left(\Delta\,{}^{2}\\!P^{2}+\Delta{\,{}^{2}\\!P^{0}}^{2}\right)\right\\}\sqrt{g}\,\text{d}x$

where the field $f(x)=\langle f\rangle=\frac{\int\Psi
f\,\text{d}[Q,A^{0},\,{}^{2}\\!P^{0}]}{\int\Psi\,\text{d}[Q,A^{0},\,{}^{2}\\!P^{0}]}$
is the average of the variable $f$, and $\Delta f^{2}=\langle
f^{2}\rangle-\langle f\rangle^{2}$ is the variance, which include spatial
derivative as well statistical, local, variance (see appendix B).

We assume that the elastic stress $\sigma_{el}^{\mu\nu}=-\frac{\delta
E}{\delta g_{\mu\nu}}$ is a fast relaxing variable, and is divergence free (in
the absence of external forces),

$\displaystyle\nabla_{\mu}\sigma^{\mu\nu}_{el}=0.$ (6)

However, the tissue may also grow and relax internal stress. This is achieved
via two type of processes. Solid-like relaxation only involves non-topological
changes of the cellular network, such as growth (change in reference area),
and elongation (change in reference perimeter term). Fluid-like relaxation,
involves topological changes and induces cellular flow and diffusion. We
consider proliferation (division and apoptosis), and $T_{1}$ transitions.

When a tissue behaves as a solid, a coordinate system is always well defined,
and tissue deformations correspond directly to changes in the metric
$g_{\mu\nu}$. These may be caused either by external forces, $F^{\mu}_{ext}$,
so that $\nabla_{\nu}\sigma^{\mu\nu}_{el}=F^{\mu}_{ext}$, or by internal
changes in the reference values, which can be calculated by solving (6), as
$\sigma^{\mu\nu}_{el}$ evolves in time with $A^{0}$ and $\,{}^{2}\\!P^{0}$.

As cells flow, i.e. - change their relative positions, one may need to
redefine the coordinate system of the tissue. In the absence of any external
reference, this may seem as a non-trivial feat. Yet, if topological changes
occur in a small region, while neighboring tissue does not, one may keep the
coordinates on the non-changing surrounding tissue, and deduce that the
coordinates in the region of transition have remained the same as well (see
fig. 1 (b)-(d)). Rather, the quantity that has changed is the network tensor,
$Q^{\mu\nu}(x)$, at that point. This means that eq. (6) is still valid and the
metric is still well defined, even when $\sigma_{el}^{\mu\nu}$ non - trivially
changes via its dependence on the network tensor $Q^{\mu\nu}$. Note, that
within this Lagrangian view, the average cellular flow, as measured by a
Lagrangian (co - moving) viewer is zero. To put simply, this happens since the
mean flow involves large configuration deformations which may be regarded as
change in the metric $g_{\mu\nu}$. Flow is thus described by change in the
energy minimizing metric (which solves eq. (6))

Figure 1:  (a) A cellular tissue before any transition happens, it is easy to
assign coordinates and metric. (b)-(d) A topoligical transition (in this case
a $T_{1}$), outside obscured region the system remains a solid. (d) The final
netwrok.

Tissue fluid-dynamics are described using a Fokker - Plank equation for the
coarse grained density function $\Psi$, over coordinates and variables (a
complete expression, including solid processes, is given in appendixes B, C):

$\displaystyle\partial_{t}\Psi=$
$\displaystyle\sum_{n}\left\\{\nabla_{\mu}\nabla_{\nu}\left[D_{n,x}^{\mu\nu}\Psi\right]-\frac{\partial}{\partial
Q^{\mu\nu}}\left[\partial_{t}Q^{\mu\nu}_{n}\Psi\right]\right.$ (7)
$\displaystyle+\left.\frac{\partial^{2}}{\partial Q^{\mu\nu}\partial
Q^{\lambda\sigma}}\left[D_{i,Q}^{\mu\nu\lambda\sigma}\Psi\right]\right\\}.$

The first term is the usual Fokker - Plank diffusion term with vanishing mean
flow, $D_{n,x}^{\mu\nu}$ is the diffusion coefficient stemming from the nth
type topological process. The second term accounts for the change in network
topology, the third term is noise. The sum is taken over division ($d$),
apoptosis ($a$), and $T_{1}$ transitions.

Eq. (7) must be completed by the requirement that the local cellular density
changes when cells divide or die:
$\partial_{t}\rho=\int\partial_{t}\Psi\,\text{d}[A^{0},\,{}^{2}\\!P^{0},\\{\ell\\}]+\gamma\rho-\frac{\rho}{2}g^{\mu\nu}\partial_{t}g_{\mu\nu}+\zeta_{\rho}$,
where the proliferation rate is
$\gamma=\frac{1}{2}Q^{-1}_{\mu\nu}\left(\partial_{t}Q^{\mu\nu}_{d}+\partial_{t}Q^{\mu\nu}_{a}\right)$.
We consider a white noise with statistics
$\langle\zeta_{\rho}(t,x)\zeta_{\rho}(t^{\prime},x^{\prime})\rangle=D_{\gamma}\rho\delta(t-t^{\prime})\delta(x-x^{\prime})$
directly related to proliferation noise $D_{n,Q}^{\mu\nu\alpha\beta}$.

Using eq. (7) one can write the energy generation rate

$\displaystyle\frac{dE}{dt}=$
$\displaystyle\int\left\\{\mu\gamma+\eta_{\mu\nu}\partial_{t}Q^{\mu\nu}-\sigma_{el}^{\mu\nu}\partial_{t}g_{\mu\nu}+\eta_{D}\right\\}\rho\sqrt{g}\,\text{d}x,$
(8)

where $\mu$ may be interpreted as a cell chemical potential,
$\sigma_{el}^{\mu\nu}$ is the elastic stress, $\eta_{\mu\nu}$ is the
dissipation associated with the change in network, and $\eta_{D}$ is the a
dissipation term due to diffusion (see appendix C for their full expressions).
Note again that the proliferation rate $\gamma$ is directly related to
$\partial_{t}Q^{\mu\nu}$ and is not an independent variable. Using an Onsager
approach, we may now relate the unknown flux ($\partial_{t}Q^{\mu\nu}$) to the
thermodynamic forces $\mu,\sigma_{el},\eta_{\mu\nu}$ which have an explicit
expression using in terms of our variables.

In the following, we make two important simplifications. First, we perform a
mean field approximation considering the local average values of all fields
and we ignore terms involving the spatial derivatives of the fields. Also, for
simplicity we assume that the fields are uniform, limiting ourselves to planar
configurations. Second, we use a minimal relaxation model in which the Onsager
coefficient matrix is diagonal, with the exception of the active
contributions, which are added to all fluxes. A more complex treatment is left
for future work.

The mean energy per cell reads

$\displaystyle E/N_{cell}=$ $\displaystyle
k_{A}\left(\sqrt{gQ}-A^{0}\right)^{2}+k_{P}\left(g_{\mu\nu}Q^{\mu\nu}-\,{}^{2}\\!P^{2}\right)^{2}$
(9) $\displaystyle+E_{\mu\nu\alpha\beta}\left\langle\Delta Q^{\mu\nu}\Delta
Q^{\alpha\beta}\right\rangle,$

where $N_{cell}=\frac{\int\sqrt{g}\,\text{d}S}{A_{cell}}$, $\Delta
E_{\mu\nu\alpha\beta}=k_{A}\sqrt{gQ}\Huge[\Large(\sqrt{gQ}-\frac{1}{4}\Large)\Large(Q^{-1}_{\mu\nu}Q^{-1}_{\alpha\beta}-Q^{-1}_{\mu\alpha}Q^{-1}_{\nu\beta}\Large)+\frac{1}{4}Q^{-1}_{\mu\alpha}Q^{-1}_{\nu\beta}\Huge]+k_{P}g_{\mu\nu}g_{\alpha\beta}$
is the expansion ofthe energy to second order in the fluctuations of
$Q^{\mu\nu}$ with respect to its average value. Here
$N_{cell}\propto\frac{1}{\sqrt{Q}}$ is the total number of cells in the tissue
at a given time.

In our minimal scheme, Onsager relations are obtained for the network tensor,
as the metric tensor is considered a fast variable. The dynamics of the
network tensor are composed of a relaxation term and an active term:

$\displaystyle\partial_{t}Q^{\mu\nu}=-\Gamma^{\mu\nu\alpha\beta}\eta_{\alpha\beta}-\xi
Q^{\mu\nu}.$ (10)

The first term is relaxation, with Onsager coefficient
$\Gamma^{\mu\nu\alpha\beta}$ and the conjugate force to the rate of change of
the network tensor
$\eta_{\alpha\beta}=\frac{1}{2\sqrt{Q}}\left[k_{A}\left(gQ-{A^{0}}^{2}\right)-k_{P}\left(g_{\mu\nu}Q^{\mu\nu}-\,{}^{2}\\!P^{0}\right)^{2}\right]Q^{-1}_{\alpha\beta}+\frac{2}{\sqrt{Q}}k_{P}\left(g_{\mu\nu}Q^{\mu\nu}-\,{}^{2}\\!P^{0}\right)g_{\alpha\beta}$.
The second active term describes a constant proliferation rate depending on
external energy supply, with $\xi>0$ for a growing tissue. A homeostatic
state, is reached whenever $\partial_{t}Q=0$, as this indicates that the
network, on average, does not change. The Onsager tensor,
$\Gamma^{\mu\nu\alpha\beta}$ may depend strongly on the topological process
involved (appendix D). The simplest tensorial form allowed by symmetry is

$\displaystyle\Gamma^{\mu\nu\alpha\beta}=H_{1}\left(\frac{1}{2}\left(Q^{\mu\alpha}Q^{\nu\beta}+Q^{\nu\alpha}Q^{\mu\beta}\right)+HQ^{\alpha\beta}Q^{\mu\nu}\right),$
(11)

where the rate constant $H_{1}$ is taken identical for all the topological
processes (and can be set to $\upsilon=1$). The dimensionless parameter
$H\geq-\frac{1}{2}$ is process dependent. For a pure $T_{1}$ transition, in
which the cell area does not change on average, $H=-\frac{1}{2}$, any
$H>-\frac{1}{2}$ corresponds to some proliferation. While other forms of
$\Gamma^{\mu\nu\alpha\beta}$ are allowed, it is only this form with
$H=-\frac{1}{2}$ that allows for ”pure shear” where the area does not change,
in any geometry. It is therefore very useful to use this same form in other
cases as well. Additionally, it allows for a simple ”one constant” (or
average) limit by considering a single $\Gamma^{\mu\nu\alpha\beta}$ describing
all processes at once.

#### Mechanics-

We first study a free, solid, tissue, for which $Q^{\mu\nu}$ does not change.
Any deformation is then due to a change of the metric $g_{\mu\nu}$. The
elastic stress in the tissue vanishes, so that $\sigma_{el}^{\mu\nu}=0$ is an
energy minimizing solution. Since both $Q^{\mu\nu}$ and $g_{\mu\nu}$ are
positive definite tensors, one can work in a coordinate system in which
$Q^{\mu\nu}=\delta^{\mu\nu}$ . Additionally, without loss of generality, one
can choose the energy scale so that $k_{A}=1$ and the length scale so that
$A^{0}=1$. Solving for $g_{\mu\nu}$, yields, the eigenvalues of $g_{\mu\nu}$,
$g_{1}$ and $g_{2}$:

$\displaystyle\left.\begin{array}[]{cc}g_{1}=\frac{1}{g_{2}}=\frac{\,{}^{2}\\!P^{0}\pm\sqrt{{\,{}^{2}\\!P^{0}}^{2}-4}}{2},&\,{}^{2}\\!P^{0}\geq
2,\\\ \
g_{1}={g_{2}}=\frac{1+2k_{P}\,{}^{2}\\!P^{0}}{1+4k_{P}},&\,{}^{2}\\!P^{0}<2\end{array}\right.$
(14)

These are the expected results known [4, 5, 9] for soft – solids (in the limit
$\,{}^{2}\\!P>2$), where cells assume an elongated shape, along a
spontaneously chosen direction, and for hard – solids when $\,{}^{2}\\!P<2$)
where cells are isotropic, and are under internal residual stress as easily
seen from the fact that $E\neq 0$. The case $\,{}^{2}\\!P=2$ corresponds to a
marginal solid.

#### Fluctuating mechanics-

The effects of fluctuations depend, to some extent, on the nature of the
disorder. We considered two archetypical examples, assuming that the
correlation between different components of $Q^{\mu\nu}$ vanish. The first one
is isotropic fluctuations, $\langle\Delta Q^{\mu\nu}\Delta
Q^{\alpha\beta}\rangle=\frac{\phi}{2}\left(g^{\mu\alpha}g^{\nu\beta}+g^{\nu\alpha}g^{\mu\beta}\right)$,
$\phi$ being the fluctuation amplitude. This is a very intuitive model,
suggesting that disorder is related to tissue shape. The second example
assumes that whatever mechanisms govern the network $Q^{\mu\nu}$ gives rise to
fluctuations, and thus the fluctuation should be should be proportional to
$\frac{\phi}{2}\left(O^{\mu\alpha}O^{\nu\beta}+O^{\nu\alpha}O^{\mu\beta}\right)$,
where $\left(O^{2}\right)^{\mu\nu}=Q^{\mu\nu}$. We term this second example
”multiplicative”. In both cases, the resulting effect is a shift of the
critical transition between an isotropic and symmetry-broken tissues to a
value $\,{}^{2}\\!P^{0}_{crit}>2$.

$\displaystyle\,{}^{2}\\!P^{0}_{crit}\simeq
2+\phi\times\left\\{\begin{array}[]{cc}\frac{1+6k_{P}}{2k_{P}}&\text{isotropic}\\\
\frac{1+16k_{P}}{8k_{P}}&\text{multiplicative}\\\ \end{array}\right.$ (17)

These results are in accordance with other similar accounts of changing the
shape parameters transition between ordered and disordered tissues [16, 6],
suggesting that fluctuations play a significant role in tissue mechanics and
dynamics.

#### External force-

Within the mean field approach, one can calculate the non-linear responses to
an external force, easily recovering the results of Ref. [5] (see appendix E),
and calculate further non-linear response coefficients such as Poisson’s
ratio.

#### Relaxation Dynamics-

The model also allows for an analytical calculation of tissue dynamics. We
first consider a tissue strongly adhered to a solid substrate, in the absence
of growth $\xi=0$. As the actual size of the tissue is assumed given, we
choose a constant metric, $g_{\mu\nu}=\delta_{\mu\nu}$, and let $Q^{\mu\nu}$
evolve. In this case, internal elastic stresses may develop, but are balanced
by forces from the substrate. Working with the eigenvalues of $Q^{\mu\nu}$,
$q_{1}$ and $q_{2}$, and using eq. (11)

$\displaystyle\dot{q}_{i}=$
$\displaystyle-\frac{\left(1+2H\right)q_{i}}{2\sqrt{q_{1}q_{2}}}\left[q_{1}q_{2}-1-k_{P}\left(q_{1}+q_{2}-\,{}^{2}\\!P^{0}\right)^{2}\right]$
(18)
$\displaystyle-\frac{2}{\sqrt{q_{1}q_{2}}}k_{p}\left(q_{1}+q_{2}-P^{0}\right)\left[q_{i}+H\left(q_{1}+q_{2}\right)\right]q_{i},$

where $i\in\left(1,2\right)$. Stable solutions depend on $\,{}^{2}\\!P^{0}$
(and therefore on the existence of residual stresses)

$\displaystyle\left.\begin{array}[]{cc}q_{1}=1/q_{2}=\frac{\,{}^{2}\\!P^{0}\pm\sqrt{{\,{}^{2}\\!P^{0}}^{2}-4}}{2},&\,{}^{2}\\!P^{0}\geq
2\\\
q_{1}=q_{2}=\sqrt{\frac{1+k_{p}{\,{}^{2}\\!P^{0}}^{2}}{1+4k_{P}}},&\,{}^{2}\\!P^{0}<2\end{array}\right.$
(21)

Interestingly, the solution for $\,{}^{2}\\!P^{0}>2$ corresponds to the same
solution (in different coordinates) as the pure mechanical problem. The
solution for $\,{}^{2}\\!P^{0}<2$ differs. The case of a tissue with no
proliferation is singular because the only topological transitions are the
pure $T_{1}$ transitions, which conserve area. It is treated in appendix F

#### Growing tissue-

A tissue grows if the active growth rate $\xi>0$. The number of cells is
$N_{cell}\propto\frac{1}{\sqrt{Q}}$ so that
$\partial_{t}N_{cells}=-\frac{N_{cells}}{2}Q^{-1}_{\mu\nu}\partial_{t}Q^{\mu\nu}$.
When considering a free tissue, one has to find both the metric $g_{\mu\nu}$
and the network tensor $Q^{\mu\nu}$, by solving equations (6) and (10)
simultaneously. For a free tissue we get $\partial_{t}N_{cells}=\gamma
N_{cells}$. The growth rate $\gamma$ depends on the existence of residual
stresses in the tissue:

$\displaystyle\gamma=\left\\{\begin{array}[]{cc}\xi,&\,{}^{2}\\!P^{0}\geq 2\\\
\xi-\frac{1+2H}{2\sqrt{Q}}\mathcal{E},&\,{}^{2}\\!P^{0}<2\end{array}\right.$
(24)

where
$\mathcal{E}=\left(\frac{1+2k_{P}\,{}^{2}\\!P^{0}}{1+4k_{P}}-1\right)^{2}+k_{P}\left(2\frac{1+2k_{P}\,{}^{2}\\!P^{0}}{1+4k_{P}}-\,{}^{2}\\!P^{0}\right)^{2}$
is the average residual energy per cell. The active growth rate, $\xi$ must be
large enough to allow continuing tissue growth, otherwise it dies under
mechanical regulation of growth. In this case, solving eq. (6), results in an
exponential change of the metric with time: $g_{\mu\nu}(t)=e^{\gamma
t}g_{\mu\nu}(0)$.

#### Confined growth-

Finally, when a tissue proliferates against an external hard boundary, the
pressure inside is expected eventually to be high enough that the death rate
cancels the active division rate at some finite pressure. Thus we search for a
homeostatic state in which $\partial_{t}Q^{\mu\nu}=0$. As the tissue is
bounded by an external wall, a coordinate system in which
$g_{\mu\nu}=\delta_{\mu\nu}$ may be chosen. As above, the existence of both
internal and residual stresses inside the tissue is important and drives
cells’ death.

In the hard solid state, $\,{}^{2}\\!P^{0}<2$, the network tensor in the
homeostatic state, $Q_{h}^{\mu\nu}$, exhibits a transition occurring at
$\xi^{*}=\frac{(1+2H)(4-{\,{}^{2}\\!P^{0}}^{2})}{4\,{}^{2}\\!P^{0}}$

$\displaystyle\left.\begin{array}[]{cc}q_{1}=q_{2}=\frac{\sqrt{\xi^{2}+(1+2H)^{2}(1+4k_{P})(1+k_{P}{\,{}^{2}\\!P^{0}}^{2})}-\xi}{(1+2H)(1+4k_{P})},&\xi<\xi^{*}\\\
q_{1,2}=\frac{\,{}^{2}\\!P^{0}}{2}\pm\frac{1}{2}\sqrt{{\,{}^{2}\\!P^{0}}^{2}-4+\frac{8\xi\left(\sqrt{\xi^{2}+\left(1+2H\right)^{2}}-\xi\right)}{(1+2H)^{2}}}&\xi>\xi^{*}\end{array}\right.$
(27)

When $\,{}^{2}\\!P^{0}\geq 2$ (no residual stress), however, there is no
transition, and $Q_{h}^{\mu\nu}$ is given by the second line of eq. (27), for
all $\xi$. In both cases, the stability of the symmetry-broken solution
($q_{2}=\,{}^{2}\\!P^{0}-q_{1}$) is dynamically driven, since for large $\xi$
values, it is energetically less favorable than the isotropic solution. In
fig. 2, we plot the homeostatic pressure, defined as
$p_{h}=-g_{\mu\nu}\sigma_{el}^{\mu\nu}$ when $Q^{\mu\nu}=Q^{\mu\nu}_{h}$.

Figure 2:  Homeostatic pressure of a tissue, as a function of the active
division rate, $\xi$, for $\,{}^{2}\\!P^{0}=1.5<2$ (blue) and
$\,{}^{2}\\!P^{0}=3>2$(red). When $\xi=0$, the pressure is negative for
$\,{}^{2}\\!P^{0}<2$ and vanishes for any $\,{}^{2}\\!P^{0}\geq 2$. Pressure
for either $\,{}^{2}\\!P^{0}>2$ or $\,{}^{2}\\!P^{0}<2$ whiel $\xi>\xi^{*}$
does not depend on $\,{}^{2}\\!P^{0}$. Here $H=0$, $k_{P}=1$, $\xi^{*}\simeq
0.2917$. $p_{h}(\xi>\xi^{*})$ is independent of $k_{P}$.

Concluding, a wide and intricate range of complex behaviors that depend
strongly on the existence of residual stresses at the cellular level is seen
in this study. These exemplify the importance of microscopics on the
macroscopic behavior of a tissue. A special focus should be given regarding
the unique predictions of the homeostatic state.

Relaxing mean field, or adding more complete Onsager relations are bound to
shed new insights to the physics of living matter. When uniformity is relaxed,
the requirement that spatial derivatives vanish can almost always be achieved
in many cases, recovering the usual defect pattern on sphere, for example.
Nevertheless, when considering evolving, non-flat geometries, one must include
bending term, which are absent in this work.

#### Acknowledgments

This paper is dedicated to the memory of a former teacher, Erez Barnoy, who
recently passed away.

## References

  * Honda [1983] H. Honda, Geometrical models for cells in tissues, International review of cytology 81, 191 (1983).
  * Fletcher  _et al._ [2014] A. G. Fletcher, M. Osterfield, R. E. Baker, and S. Y. Shvartsman, Vertex models of epithelial morphogenesis, Biophysical journal 106, 2291 (2014).
  * Farhadifar  _et al._ [2007] R. Farhadifar, J.-C. Röper, B. Aigouy, S. Eaton, and F. Jülicher, The influence of cell mechanics, cell-cell interactions, and proliferation on epithelial packing, Current Biology 17, 2095 (2007).
  * Staple  _et al._ [2010] D. B. Staple, R. Farhadifar, J.-C. Röper, B. Aigouy, S. Eaton, and F. Jülicher, Mechanics and remodelling of cell packings in epithelia, The European Physical Journal E 33, 117 (2010).
  * Moshe  _et al._ [2018] M. Moshe, M. J. Bowick, and M. C. Marchetti, Geometric frustration and solid-solid transitions in model 2d tissue, Physical review letters 120, 268105 (2018).
  * Sahu  _et al._ [2019] P. Sahu, J. Kang, G. Erdemci-Tandogan, and M. L. Manning, Nonlinear analysis of the fluid-solid transition in a model for ordered biological tissues, arXiv preprint arXiv:1905.12714 (2019).
  * Chiou  _et al._ [2012] K. K. Chiou, L. Hufnagel, and B. I. Shraiman, Mechanical stress inference for two dimensional cell arrays, PLoS computational biology 8, e1002512 (2012).
  * Merkel and Manning [2017] M. Merkel and M. L. Manning, Using cell deformation and motion to predict forces and collective behavior in morphogenesis, in  _Seminars in cell & developmental biology_, Vol. 67 (Elsevier, 2017) pp. 161–169.
  * Popović  _et al._ [2021] M. Popović, V. Druelle, N. A. Dye, F. Jülicher, and M. Wyart, Inferring the flow properties of epithelial tissues from their geometry, New Journal of Physics 23, 033004 (2021).
  * Guirao  _et al._ [2015] B. Guirao, S. U. Rigaud, F. Bosveld, A. Bailles, J. Lopez-Gay, S. Ishihara, K. Sugimura, F. Graner, and Y. Bellaïche, Unified quantitative characterization of epithelial tissue development, Elife 4, e08519 (2015).
  * Ranft  _et al._ [2010] J. Ranft, M. Basan, J. Elgeti, J.-F. Joanny, J. Prost, and F. Jülicher, Fluidization of tissues by cell division and apoptosis, Proceedings of the National Academy of Sciences 107, 20863 (2010).
  * Hannezo  _et al._ [2014] E. Hannezo, J. Prost, and J.-F. Joanny, Theory of epithelial sheet morphology in three dimensions, Proceedings of the National Academy of Sciences 111, 27 (2014).
  * Murisic  _et al._ [2015] N. Murisic, V. Hakim, I. G. Kevrekidis, S. Y. Shvartsman, and B. Audoly, From discrete to continuum models of three-dimensional deformations in epithelial sheets, Biophysical journal 109, 154 (2015).
  * Merkel  _et al._ [2017] M. Merkel, R. Etournay, M. Popović, G. Salbreux, S. Eaton, and F. Jülicher, Triangles bridge the scales: Quantifying cellular contributions to tissue deformation, Physical Review E 95, 032401 (2017).
  * Gorban [2006] A. N. Gorban, Basic types of coarse-graining, in  _Model Reduction and Coarse-Graining Approaches for Multiscale Phenomena_ (Springer, 2006) pp. 117–176.
  * Bi  _et al._ [2015] D. Bi, J. Lopez, J. M. Schwarz, and M. L. Manning, A density-independent rigidity transition in biological tissues, Nature Physics 11, 1074 (2015).

##  Appendix A Appendix A - Energy Functional

Typically, the energy of the vertex model is written:

$\displaystyle E$
$\displaystyle=\sum_{i}\hat{k}_{A}\left(A_{i}-\hat{A}^{0}_{i}\right)^{2}+\tilde{k}_{P}\left(P_{i}-P^{0}_{i}\right)^{2},$
(28)

Focusing on the perimeter term:

$\displaystyle E_{P}$
$\displaystyle=\sum_{i}\tilde{k}_{P}\left(P_{i}-P^{0}_{i}\right)=\sum_{i}\tilde{K}_{P}\left(\frac{P-P^{0}}{P^{0}}\right)^{2}$
(29)

where we defined $\tilde{K}_{P}=\tilde{k}_{P}{P^{0}}^{2}$ is the elastic
modulus, and $\frac{P-P^{0}}{P^{0}}$ is the generalized strain along the
perimeter. Since the strain is uniform along the perimeter, one can brake it
onto smaller segments

$\displaystyle\left(\frac{P-P^{0}}{P^{0}}\right)^{2}=\sum_{e}\left(\frac{\ell_{e}-\ell_{e}^{0}}{\ell_{e}^{0}}\right)^{2}$
(30)

where $\ell_{e}$ is the actual length of an edge, and $\ell_{e}^{0}$ is a
”reference” length, satisfying $\sum_{e}\ell_{e}^{0}=P^{0}$. Thus we may write

$\displaystyle E_{P}$ $\displaystyle=\sum_{i}\sum_{e\in
i}\tilde{K}_{P}\left(\frac{\ell_{e}-\ell_{e}^{0}}{\ell_{e}^{0}}\right)^{2},$
(31)

where the sum is taken over all edges associated with the $i^{th}$ cell. For
small strains such that $\ell_{e}\sim\ell_{e}^{0}$ one may write approximate

$\displaystyle\frac{\ell_{e}-\ell_{e}^{0}}{\ell_{e}^{0}}\simeq\frac{\ell_{e}^{2}-{\ell_{e}^{0}}^{2}}{2{\ell_{e}^{0}}^{2}}.$
(32)

Thus, finally

$\displaystyle E_{P}$ $\displaystyle=\sum_{i}\sum_{e\in
i}\tilde{K}_{P}\left(\frac{\ell_{e}-\ell_{e}^{0}}{\ell_{e}^{0}}\right)^{2}\simeq\sum_{i}\sum_{e\in
i}\tilde{K}_{P}\left(\frac{\ell_{e}^{2}-{\ell_{e}^{0}}^{2}}{4{\ell_{e}^{0}}^{2}}\right)^{2}$
(33)
$\displaystyle\sum_{i}\frac{\tilde{K}_{P}}{4{\,{}^{2}\\!P^{0}}^{2}}\left(\frac{\,{}^{2}\\!P_{i}-{\,{}^{2}\\!P_{i}^{0}}^{2}}{{\,{}^{2}\\!P_{i}^{0}}^{2}}\right)^{2}$

where used the same uniform strain condition as before. This time under the
constraint $\,{}^{2}\\!P^{0}=\sum_{e}{\ell_{e}^{0}}^{2}$. Redefining
$k_{P}=\frac{\tilde{K}_{P}}{4{\,{}^{2}\\!P^{0}}^{2}}=\frac{\tilde{k}_{P}{{P^{0}}^{2}}}{4{\,{}^{2}\\!P^{0}}^{2}}$,
results in eq. (1). Since we assume $\ell_{0}$ are given at this point, this
should not pose any problem. However, eventually we relax this condition and
only limit the sum. While limiting $\,{}^{2}\\!P^{0}$ is not the same as
limiting $P^{0}$, as both are distances after all.

###  A.1 Appendix B - Derivatives in Coarse-Grained models

Coarse-graining is the averaging of quantities over some scale $\ell$. This is
done by taking a kernel function $\phi_{\ell}(x)=\phi_{\ell}(-x)$ so that
$\int\phi_{\ell}(x)dx=1$ and
$\phi_{\ell}(x)\xrightarrow{x\rightarrow\infty}0$. Coarse-graining a function
$f(x)$ is given by the convulsion with the kernel-

$\displaystyle f_{\ell}(x)=\int\phi_{\ell}(y)f(x-y)dy.$ (34)

Assuming $\phi_{\ell}$ decays fast enough, relative to the scale of change of
$f(x)$, one can write to second order in $\ell$-

$\displaystyle f_{\ell}(x)$
$\displaystyle=\int\phi_{\ell}(y)f(x-y)dy\simeq\int\phi_{\ell}(y)\left[f(x)-y\partial_{x}f(x)+\frac{1}{2}y^{2}\left(\partial_{x}f(x)\right)^{2}\right]dy=f(x)+C\ell^{2}\left(\partial_{x}f(x)\right)^{2},$
(35)

where $C$ is some constant, assumed small.

The coarse-graining of the square of function gives a similar result-

$\displaystyle(f^{2})_{\ell}(x)$
$\displaystyle=\int\phi_{\ell}(y)f^{2}(x-y)dy\simeq\int\phi_{\ell}(y)\left(f(x)-y\partial_{x}f(x)\right)^{2}dy=f^{2}(x)+C_{2}\ell^{2}\left(\partial_{x}f(x)\right)^{2},$
(36)

where $C_{2}$ is again some constant.

When coarse-graining the vertex model, we essentially start from a very
singular cell distribution - the density, for example, is given by

$\rho(x)=\sum\delta(x-x_{i}),$

$x_{i}$ being the position of the $i^{th}$ cell. Thus we assumes some scale
over which a system realization is already coarse-grained. While in principle
$\ell$ is an arbitrary scale, when it is large enough, the same expressions
above results with a natural correlation scale of the system (if exists) .
Therefore

$\displaystyle\langle
A^{2}\rangle=\langle(A-A(x)+A(x))^{2}\rangle=\langle(\Delta
A+A(x))^{2}\rangle=\langle\Delta
A^{2}\rangle+\langle(A(x))^{2}\rangle=\langle\Delta
A^{2}\rangle+(A(x))^{2}+C_{3}\ell^{2}\left(\partial_{x}A(x)\right)^{2},$ (37)

where $C_{3}$ is some constant, and $\Delta$ .

###  A.2 Appendix C - Full Coarse-Grained model

The full Fokker - Planck equation is given by:

$\displaystyle\partial_{t}\Psi=$ $\displaystyle\frac{\partial^{2}}{\partial
x^{\mu}\partial x^{\nu}}\left[D_{x}^{\mu\nu}\Psi\right]$ (38)
$\displaystyle-\frac{\partial}{\partial
A^{0}}\left[\tilde{K}\Psi\right]+\frac{\partial^{2}}{(\partial
A^{0})^{2}}\left[D_{\tilde{K}}\Psi\right]$
$\displaystyle-\frac{\partial}{\partial\,{}^{2}\\!P^{0}}\left[\tilde{\varepsilon}\Psi\right]+\frac{\partial^{2}}{(\partial\,{}^{2}\\!P^{0})^{2}}\left[D_{\tilde{\varepsilon}}\Psi\right]$
$\displaystyle+\sum_{i}\left\\{-\frac{\partial}{\partial
Q^{\mu\nu}}\left[\gamma_{i}\,\delta_{i}Q^{\mu\nu}\Psi\right]+\frac{\partial^{2}}{\partial
Q^{\mu\nu}\partial
Q^{\lambda\sigma}}\left[D_{i,Q}^{\mu\nu\lambda\sigma}\Psi\right]\right\\}.$

$D_{x}^{\mu\nu}$ is the total diffusion coefficient stemming from the
different topological processes possible. The next lines describe the effect
of noisy growth and topological transitions on the density function. Sum is
taken over division ($d$), apoptosis ($a$), and $T_{1}$ transitions.
$\tilde{K}$ is the average growth rate of a cell, $\tilde{\varepsilon}$ is the
average elongation, and $\partial_{t}Q^{\mu\nu}$ is the average change in the
network. $D_{\tilde{K}/\tilde{\varepsilon}/Q}$’s relate to noise.

Beginning from eq. (38) , we can compute the following dynamic equations for
the variables $\alpha=\int A^{0}\Psi(x,\chi)\,\text{d}\chi$ (where
$\chi=\\{A^{0},\,{}^{2}\\!P^{0},Q\\}$),
$\,{}^{2}\\!\tilde{\varpi}=\frac{\int\,{}^{2}\\!P_{0}\Psi\,\text{d}\chi}{g_{\mu\nu}Q^{\mu\nu}\int\Psi\,\text{d}\chi}$

$\displaystyle\frac{\partial\alpha}{\partial t}$
$\displaystyle=\left(\gamma_{a}-\gamma_{d}\right)\alpha+\nabla_{\mu}\nabla_{\nu}\left[D_{x}^{\mu\nu}\alpha\right]+\tilde{K}\rho-\frac{\alpha}{2}g^{\mu\nu}\partial_{t}g_{\mu\nu}$
(39) $\displaystyle\frac{\partial\,{}^{2}\\!\tilde{\varpi}}{\,\text{d}t}$
$\displaystyle=\frac{\partial}{\partial
t}\left(\frac{\,{}^{2}\\!P^{0}}{\,{}^{2}\\!P}\right)=\frac{\partial_{t}{\,{}^{2}\\!P}^{0}}{\,{}^{2}\\!P}-\frac{\,{}^{2}\\!P^{0}}{\,{}^{2}\\!P}\frac{\partial_{t}g_{\mu\nu}Q^{\mu\nu}+g_{\mu\nu}\partial_{t}Q^{\mu\nu}}{g_{\mu\nu}Q^{\mu\nu}}$
(40)
$\displaystyle=\frac{1}{\,{}^{2}\\!P}\left(\tilde{\varepsilon}+\frac{\nabla_{\mu}\nabla_{\nu}\left(D^{\mu\nu}_{x}\rho\,{}^{2}\\!P^{0}\right)-\,{}^{2}\\!P^{0}\nabla_{\mu}\nabla_{\nu}\left(D^{\mu\nu}_{x}\rho\right)}{\rho}\right)-\,{}^{2}\\!\tilde{\varpi}\frac{\partial_{t}g_{\mu\nu}Q^{\mu\nu}+g_{\mu\nu}\partial_{t}Q^{\mu\nu}}{g_{\mu\nu}Q^{\mu\nu}}$

The metric $g$ is given by the instantaneous solution for the equation

$\nabla_{\mu}\sigma^{\mu\nu}_{el}=f_{ext}^{\nu}$

where $\sigma^{\mu\nu}_{el}=\frac{\delta E}{\delta g_{\mu\nu}}$ is the stress
tensor, and $f_{ext}^{\nu}$ is an external force acting on the tissue. Thus in
each time step we solve for $g$ and allow $E$ to change via the change in
$\Psi$.

Confluence of cells means $\rho=\frac{1}{\langle A\rangle}$. Using this
relation and the time dependence of the second moments-

$\displaystyle\partial_{t}\langle A_{0}^{2}\rangle$ $\displaystyle=2\langle
A_{0}\rangle\tilde{K}+2D_{\tilde{K}}+\frac{1}{\rho}\left[\nabla_{\mu}\nabla_{\nu}\left(D^{\mu\nu}_{x}\langle
A_{0}^{2}\rangle\rho\right)-\langle
A_{0}^{2}\rangle\nabla_{\mu}\nabla_{\nu}\left(D^{\mu\nu}_{x}\rho\right)\right]$
(42) $\displaystyle\partial_{t}\langle A^{2}\rangle$ $\displaystyle=\langle
A\rangle^{2}\left\\{g^{\mu\nu}\partial_{t}g_{\mu\nu}-2\gamma-\frac{2}{\rho}\nabla_{\mu}\nabla_{\nu}\left[D^{\mu\nu}_{x}\rho\right]\right\\}+\partial_{t}\left(\frac{\langle\Delta\rho^{2}\rangle}{\rho^{4}}\right)$
(43) $\displaystyle\partial_{t}\langle(\,{}^{2}\\!P^{0})^{2}\rangle$
$\displaystyle=2\langle\,{}^{2}\\!P^{0}\rangle\tilde{\varepsilon}+2D_{\tilde{\varepsilon}}+\frac{1}{\rho}\left[\nabla_{\mu}\nabla_{\nu}\left(D^{\mu\nu}_{x}\langle(\,{}^{2}\\!P^{0})^{2}\rangle\rho\right)-\langle(\,{}^{2}\\!P^{0})^{2}\rangle\nabla_{\mu}\nabla_{\nu}\left(D^{\mu\nu}_{x}\rho\right)\right]$
(44) $\displaystyle\partial_{t}\langle\,{}^{2}\\!P^{2}\rangle$
$\displaystyle=2\left(g_{\mu\nu}\partial_{t}g_{\alpha\beta}+2g_{\mu\nu}g_{\alpha\beta}\right)\langle
Q^{\mu\nu}\partial_{t}Q^{\alpha\beta}\rangle+2g_{\mu\nu}g_{\alpha\beta}D_{Q}^{\mu\nu\alpha\beta},$
(45)

we can write the energy production rate. It is given fully by

$\displaystyle\frac{dE}{dt}=\int\left\\{\left(\frac{\partial\mathcal{E}}{\partial
t}\right)_{g}+\left(\frac{\partial\mathcal{E}}{\partial
t}\right)_{T_{1}}+\left(\frac{\partial\mathcal{E}}{\partial
t}\right)_{\gamma}+\left(\frac{\partial\mathcal{E}}{\partial
t}\right)_{growth}+\left(\frac{\partial\mathcal{E}}{\partial
t}\right)_{disorder}+\left(\frac{\partial\mathcal{E}}{\partial
t}\right)_{diff}\right\\}\rho\sqrt{g}\,\text{d}x$ (46)

where

$\displaystyle\left(\frac{\partial\mathcal{E}}{\partial t}\right)_{g}$
$\displaystyle=$
$\displaystyle\,\left[K_{A}\left(1-\alpha\right)g^{\mu\nu}+2K_{P}\left(1-\,{}^{2}\\!\tilde{\varpi}\right)\frac{Q^{\mu\nu}}{\,{}^{2}\\!P}\right]\partial_{t}g_{\mu\nu}$
(47) $\displaystyle\left(\frac{\partial\mathcal{E}}{\partial
t}\right)_{T_{1}}$ $\displaystyle=$
$\displaystyle\,2K_{P}\left(1-\,{}^{2}\\!\tilde{\varpi}\right)\frac{g_{\mu\nu}}{\,{}^{2}\\!P}~{}\left(\partial_{t}Q^{\mu\nu}\right)_{T_{1}}$
(48) $\displaystyle\left(\frac{\partial\mathcal{E}}{\partial
t}\right)_{\gamma_{d/a}}$ $\displaystyle=$
$\displaystyle\,\left[-K_{A}\left(1-\alpha^{2}\right)+K_{P}\left(1-\,{}^{2}\\!\tilde{\varpi}\right)^{2}\right]\gamma_{d/a}+2K_{P}\left(1-\,{}^{2}\\!\tilde{\varpi}\right)\frac{g_{\mu\nu}}{\,{}^{2}\\!P}\left(\partial_{t}Q^{\mu\nu}\right)_{\gamma_{d/a}}$
(49) $\displaystyle\left(\frac{\partial\mathcal{E}}{\partial
t}\right)_{growth}$ $\displaystyle=$
$\displaystyle\,-2K_{A}\left(1-\alpha\right)\rho\tilde{K}-2K_{P}\left(1-\,{}^{2}\\!\tilde{\varpi}\right)\frac{\varepsilon}{\,{}^{2}\\!P}+2K_{A}D_{\tilde{K}}\rho^{2}+2K_{P}\frac{D_{\tilde{\varepsilon}}}{\,{}^{2}\\!P^{2}}$
(50) $\displaystyle\left(\frac{\partial\mathcal{E}}{\partial
t}\right)_{disorder}$ $\displaystyle=$
$\displaystyle\,-K_{A}\rho^{2}\partial_{t}\left(\frac{\langle\Delta\rho^{2}\rangle}{\rho^{4}}\right)+2K_{P}\frac{\langle\Delta
Q^{\mu\nu}\Delta
Q^{\alpha\beta}\rangle}{\,{}^{2}\\!P^{2}}g_{\mu\nu}\partial_{t}g_{\alpha\beta}+2K_{p}\frac{g_{\mu\nu}g_{\alpha\beta}}{\,{}^{2}\\!P^{2}}\langle\Delta
Q^{\mu\nu}\partial_{t}Q^{\alpha\beta}\rangle$ (51)
$\displaystyle\,+\left[K_{A}\rho^{2}\left(\langle\Delta\rho^{2}\rangle+\langle\Delta
A_{0}^{2}\rangle\right)+K_{P}\frac{1}{\,{}^{2}\\!P^{2}}\left(\langle\Delta\,{}^{2}\\!P^{2}\rangle+\langle\Delta\,{}^{2}\\!P_{0}^{2}\rangle\right)\right]\left(\gamma_{d}-\gamma_{a}\right)$
$\displaystyle\left(\frac{\partial\mathcal{E}}{\partial t}\right)_{diff}$
$\displaystyle=$
$\displaystyle\,\left(-K_{A}\left(1-\alpha^{2}\right)+K_{P}\left(1-\,{}^{2}\\!\tilde{\varpi}\right)^{2}\right)\frac{1}{\rho}\nabla_{\mu}\nabla_{\nu}\left(D^{\mu\nu}\rho\right)$
(52) $\displaystyle\,+K_{A}\left[\rho\mathcal{D}\left(\langle
A_{0}^{2}\rangle\right)-2\mathcal{D}\left(\langle
A_{0}\rangle\right)\right]+K_{P}\left[\frac{1}{\rho\,{}^{2}\\!P^{2}}\mathcal{D}\left(\langle
P_{0}^{2}\rangle\right)-2\frac{1}{\rho\,{}^{2}\\!P}\mathcal{D}\left(\langle
P_{0}\rangle\right)\right]$

Here $K_{A}=k_{A}gQ$, $K_{P}=k_{P}\left(g_{\mu\nu}Q^{\mu\nu}\right)^{2}$, and
$\gamma_{i}=-\frac{1}{2}Q^{-1}_{\mu\nu}\partial_{t}Q^{\mu\nu}_{i}$.
Additionally defined the ”covariant” diffusion functional, of a scalar $f$ is

$\displaystyle\mathcal{D}\left(f\right)=\nabla_{\mu}\nabla_{\nu}\left(D^{\mu\nu}\rho
f\right)-f\nabla_{\mu}\nabla_{\nu}\left(D^{\mu\nu}\rho\right),$ (53)

which describes how a value $f$, carried by the cell density, diffuses
regardless of the effect of diffusion on the density $\rho$. Indeed, when $f$
is position independent, this terms is zero.

Equation (46) can be written in a simple form

$\displaystyle\frac{dE}{dt}=\int\left\\{\gamma\mu-\sigma_{el}^{\mu\nu}\partial_{t}g_{\mu\nu}-\sigma_{A}\tilde{K}-\sigma_{P}\tilde{\varepsilon}+\eta_{\beta\gamma}\partial_{t}Q^{\mu\nu}+\eta_{D}\right\\}\rho\sqrt{g}\,\text{d}x$
(54)

where $\mu$ can be interpreted as a ”chemical potential”,
$\sigma_{el}^{\mu\nu}$ is the elastic stress, $\sigma_{A}$ is an active
”pressure”, $\sigma_{P}$ is an active ”line tension”, $\eta_{\beta\gamma}$ is
the flow - stress, and $\eta_{D}$ is the energy dissipation due to diffusion,
and noise. They are given by (omitting ”disorder” terms for simplicity)

$\displaystyle\mu$
$\displaystyle=-K_{A}\left(1-\alpha^{2}\right)+K_{P}\left(1-\,{}^{2}\\!\tilde{\varpi}\right)^{2}+-2K_{P}\left(1-\,{}^{2}\\!\tilde{\varpi}\right)\frac{g_{\mu\nu}}{\,{}^{2}\\!P}\frac{\left(\delta
Q^{\mu\nu}\right)}{Q^{-1}_{\alpha\beta}\partial_{t}Q^{\alpha\beta}}$ (55)
$\displaystyle\sigma_{el}^{\mu\nu}$
$\displaystyle=-K_{A}\left(1-\alpha\right)g^{\mu\nu}-2K_{P}\left(1-\,{}^{2}\\!\tilde{\varpi}\right)\frac{Q^{\mu\nu}}{\,{}^{2}\\!P}$
(56) $\displaystyle\sigma_{A}$
$\displaystyle=-2K_{A}\rho\left(1-\alpha\right)$ (57)
$\displaystyle\sigma_{P}$
$\displaystyle=-2K_{P}\frac{1}{\,{}^{2}\\!P}\left(1-\,{}^{2}\\!\tilde{\varpi}\right)$
(58) $\displaystyle\eta_{\mu\nu}$
$\displaystyle=2K_{P}\left(1-\,{}^{2}\\!\tilde{\varpi}\right)\frac{g_{\mu\nu}}{\,{}^{2}\\!P}$
(59) $\displaystyle\eta_{D}$
$\displaystyle=\left(-K_{A}\left(1-\alpha^{2}\right)+K_{P}\left(1-\,{}^{2}\\!\tilde{\varpi}\right)^{2}\right)\frac{1}{\rho}\nabla_{\mu}\nabla_{\nu}\left(D^{\mu\nu}\rho\right)$
(60) $\displaystyle\,+K_{A}\left[\rho\mathcal{D}\left(\langle
A_{0}^{2}\rangle\right)-2\mathcal{D}\left(\langle
A_{0}\rangle\right)\right]+K_{P}\left[\frac{1}{\rho\,{}^{2}\\!P^{2}}\mathcal{D}\left(\langle
P_{0}^{2}\rangle\right)-2\frac{1}{\rho\,{}^{2}\\!P}\mathcal{D}\left(\langle
P_{0}\rangle\right)\right]$

Using an Onsager approach, we may now relate the unknown fluxes
($\partial_{t}Q^{\mu\nu},\tilde{K},\tilde{\varepsilon}$) to the potentials
$\mu,\sigma_{el},\eta_{\mu\nu}$ which have an explicit expression using our
variables. Such models in general could be very complex.

###  A.3 Appendix D -Choice of $\Gamma^{\mu\nu\alpha\beta}$

At this point we are missing explicit relations between the rates and and the
stress (or other expressions). One way to do so be to apply a minimal Onsager
appraoch so that the rate of change is linear in its conjugate, with the least
couplings. Thus, a naive guess would be

$\displaystyle\left(\partial_{t}Q^{\mu\nu}\right)_{i}=\Gamma^{\mu\nu\alpha\beta}_{i}\frac{\delta
E}{\delta Q^{\alpha\beta}}$ (61)

where the second term appear since $\rho$ is (non-trivially) related to
$Q^{\mu\nu}$. $\Gamma_{i}^{\mu\nu\alpha\beta}$ is a symmetric Onsager tensor
relating the geometry to changes in network.

It is immediately clear that $\Gamma^{\mu\nu\alpha\beta}$ is symmetric to
exchange of $\mu\leftrightarrow\nu$ and $\alpha\leftrightarrow\beta$. From the
contribution to the energy generation rate, it is also clear that it must be
symmetric to $(\mu,\nu)\leftrightarrow(\alpha,\beta)$.

The most general expression satisfying the these symmetries is

$\displaystyle\Gamma^{\mu\nu\alpha\beta}$
$\displaystyle=\frac{H_{1}}{2}\left(Q^{\mu\alpha}Q^{\nu\beta}+Q^{\nu\alpha}Q^{\mu\beta}\right)+H_{2}Q^{\alpha\beta}Q^{\mu\nu}$
(62)
$\displaystyle+\frac{H_{3}}{2}\left((Q^{-1})^{\mu\alpha}(Q^{-1})^{\nu\beta}+(Q^{-1})^{\nu\alpha}(Q^{-1})^{\mu\beta}\right)+H_{4}(Q^{-1})^{\alpha\beta}(Q^{-1})^{\mu\nu}$
$\displaystyle+\frac{H_{4}}{4}\left(Q^{\mu\alpha}(Q^{-1})^{\nu\beta}+(Q^{-1})^{\mu\alpha}Q^{\nu\beta}+Q^{\nu\alpha}(Q^{-1})^{\mu\beta}+(Q^{-1})^{\nu\alpha}Q^{\mu\beta}\right)+\frac{H_{5}}{2}\left(Q^{\alpha\beta}(Q^{-1})^{\mu\nu}+(Q^{-1})^{\alpha\beta}Q^{\mu\nu}\right)$
$\displaystyle+\frac{G_{1}}{2}\left(g^{\mu\alpha}g^{\nu\beta}+g^{\nu\alpha}g^{\mu\beta}\right)+G_{2}g^{\alpha\beta}g^{\mu\nu}$
$\displaystyle+\frac{G_{3}}{4}\left(g^{\mu\alpha}Q^{\nu\beta}+Q^{\mu\alpha}g^{\nu\beta}+g^{\nu\alpha}Q^{\mu\beta}+Q^{\nu\alpha}g^{\mu\beta}\right)+\frac{G_{4}}{2}\left(g^{\alpha\beta}Q^{\mu\nu}+Q^{\alpha\beta}g^{\mu\nu}\right)$
$\displaystyle+\frac{G_{5}}{4}\left(g^{\mu\alpha}(Q^{-1})^{\nu\beta}+(Q^{-1})^{\mu\alpha}g^{\nu\beta}+g^{\nu\alpha}(Q^{-1})^{\mu\beta}+(Q^{-1})^{\nu\alpha}g^{\mu\beta}\right)+\frac{G_{6}}{2}\left(g^{\alpha\beta}(Q^{-1})^{\mu\nu}+(Q^{-1})^{\alpha\beta}g^{\mu\nu}\right)$

However, terms including the $Q^{-1}$, have negligible contributions for large
$Q$’s and are thus unlikely. In any case, just as with terms involving
$g^{\mu\nu}$, these expressions cannot contribute to a $T_{1}$ transition as
they they contribute either isotropically ($g^{\mu\nu}$) or inversely
($Q^{-1}$, making shorter dimensions even shorter ). We thus adopt a minimal
scheme-

$\displaystyle\Gamma^{\mu\nu\alpha\beta}=H_{1}\left(\frac{1}{2}\left(Q^{\mu\alpha}Q^{\nu\beta}+Q^{\nu\alpha}Q^{\mu\beta}\right)+HQ^{\alpha\beta}Q^{\mu\nu}\right).$
(63)

where $H_{2}=H_{1}H$ This suggest that $\partial_{t}Q^{\mu\nu}\propto
H_{1}(Q^{\mu}_{\alpha}Q^{\alpha\nu}+HQ^{\mu\nu}\,{}^{2}\\!P)$, thus, a
suitable choice of the coefficient $H$ can account for different transitions.
Without loss of generality, one can set $H_{1}=1$. The choice

$\displaystyle H=-1/2$ (64)

results, at the locally flat, Q - diagonal frame, where $Q=Diag[q_{1},q_{2}]$,
with
$\partial_{t}Q\propto\frac{1}{q_{1}+q_{2}}Diag[q_{1}(q_{1}-q_{2}),-q_{2}(q_{1}-q_{2})]$
and $\partial_{t}E\propto(q_{1}-q_{2})^{2}$ is expected for a $T_{1}$
transition.

Since $\rho$ and $Q$ are related, our knowledge of the dynamics of $\rho$ can
be used determine $H_{2}$ by the flow and fluctuation terms to those derived
from $\partial_{t}Q^{\mu\nu}$.

For $H=0$, at the locally flat, $Q$ \- diagonal frame $\left(\delta
Q^{\mu\nu}\right)_{\gamma}\propto Diag[q_{1}^{2},q_{2}^{2}]/(q_{1}+q_{2})$.
Note that the choice of $\left(\delta Q^{\mu\nu}\right)_{\gamma}$ results in a
negative contribution to the energy whenever $\gamma>0$ and
$\,{}^{2}\\!\tilde{\varpi}<1$ (perimeter too long), or $\gamma<0$ and
$\,{}^{2}\\!\tilde{\varpi}>1$ (too short). While for the cases of $\gamma>0$
and $\,{}^{2}\\!\tilde{\varpi}>1$, or $\gamma<0$ and
$\,{}^{2}\\!\tilde{\varpi}<1$, this contribution is negative but small. The
overall contribution must be negative of course. This is in qualitative
agreement with a realistic case (since division only shortens, and appoptosis
lengthens), but quantitatively it differs. Consider the case of division when
the cells are too short, realistically we will only change the short dimension
of the cell so as to not loose too much length, while in this modelling
shortening happens in both dimensions.

###  A.4 Appendix E - External Force

We now turn to solve the case of finite sudden strain. And we will work within
the mean field approximation (where our dynamical fields, $g_{\mu\nu}$ and
$Q^{\mu\nu}$ do not vary through space). We thus consider a square tissue,
with some initial $Q^{\mu\nu}$, located such that one edge is given by the
coordinate $y=0$, the other by $y=Y$, and similarly there’s an edge at $x=0$
and $x=X$. At time zero we set the tissue’s configuration $\vec{f}(x,y)$ so
that $\left|\vec{f}(x=X)-\vec{f}(x=0)\right|=L$. We can incorporate this
constraint to an effective energy using a Lagrange multiplier

$\displaystyle E_{ext}$
$\displaystyle=\int\lambda\left(\left|\vec{f}(x=X)-\vec{f}(x=0)\right|-L\right)\,\text{d}y=\int\,\text{d}y\lambda\left(\left|\int\limits_{0}^{X}\partial_{x}\vec{f}\,\text{d}x\right|-L\right)\,\text{d}y$
(65)
$\displaystyle=\int\lambda\left(\sqrt{\int\int\vec{f}^{\prime}(x)\cdot\vec{f}^{\prime}(x^{\prime})\,\text{d}x^{\prime}\,\text{d}x}-L\right)\,\text{d}y.$

Within the mean field approach
$\vec{f}^{\prime}(x)=\partial_{x}\vec{f}](x)=const.$ So that
$\vec{f}^{\prime}(x)\cdot\vec{f}^{\prime}(x^{\prime})=g_{11}$ we can write

$\displaystyle E_{ext}$
$\displaystyle=\int\lambda\left(\sqrt{X^{2}g_{11}}-L\right)\,\text{d}y=\int\lambda\left(X\sqrt{g_{11}}-X\sqrt{G}\right)\,\text{d}y=\int{\lambda}X\left(\sqrt{g_{11}}-\sqrt{G}\right)\,\text{d}y$
(66)
$\displaystyle=\int{\lambda}\left(\sqrt{g_{11}}-\sqrt{G}\right)\,\text{d}x\,\text{d}y$

where $G$ is the length scale of the stretched tissue.

Taking the variation with respect to $g$ an $\lambda$ we derive the following
equations:

$\displaystyle\sqrt{g_{11}}=\sqrt{G}$ (67)
$\displaystyle\sigma^{\mu\nu}=\sigma_{el}^{\mu\nu}+\frac{\lambda}{2\sqrt{g_{11}}}\left(\begin{array}[]{cc}1&0\\\
0&0\end{array}\right)$ (70)

where $\sigma_{el}^{\mu\nu}$ is the elastic stress given above. In principle,
in order to find the metric $g_{\mu\nu}$ minimizing the elastic energy we need
to find the stress so that on the boundaries $x=0$ and $x=X$, the stress
balances the force. In the mean field approximation the conditions at the
boundary impose $\sigma^{\mu\nu}=0$.

We start by assuming $Q=diag(q_{1},q_{2})$ and $g=diag(g_{11},g_{22})$, are
aligned along the principal stretching direction. Finding $g_{22}$ and
$\lambda$ by solving eq. (67) for given $q_{1},q_{2},G$, we can explicitly
solve the problem. By defining strain with respect to the free (non stretched)
tissue- $\Delta=g_{11}-g_{11}^{free}$. We find that compatible systems
($P^{0}>2$) exhibit a bi-stability corresponding to exchange of
$g_{11}^{free}$ and $g_{22}^{free}$, i.e - for $\Delta^{*}$ =
$g_{22}^{free}-g_{11}^{free}$. In figure 3(a) we plot the elastic energy as a
function of $\Delta$, for three cases - compatible case $P^{0}>2$ (green),
marginal case $P^{0}=2$ (orange), and incompatible case $P^{0}<2$. The fact
that $E(\Delta=0,P^{0}<2)\neq 0$ is a clear sign of residual stresses.

|
---|---
|

(a)

(b)

(c)

(d)

Figure 3: Tissue energy (a), $g_{22}$ (b), force (c), and Poisson’s ratio
($\nu$, d), for the compatible case (green line), marginal case (red), and
incompatible case(blue), as indicated in the legend, all with $k_{p}=1$.
Additionally, different values of $k_{p}$ are plotted (dashed orange, and
dotted black lines) for $P^{0}=2.2$ (compatible). As a function of
$\Delta=g_{11}-g_{11}^{free}$. Note that for the solid green, dashed red, and
dotted graphs (all with $\,{}^{2}\\!P_{0}>2$), there are two energetically
favored states - with $\Delta=0$, and $\Delta=\Delta^{*}$, which in this case
$\delta*\simeq 0.9165$

Poisson’s ration is typically defined as the derivative of the transverse
strain vs the longitudinal strain $\nu=-\frac{d\epsilon_{y}}{d\epsilon_{x}}$,
near mechanical equilibrium. In our notation this translates into

$\displaystyle\nu$ $\displaystyle=-\frac{d\log g_{22}}{d\log\Delta}$ (71)

$\Delta=g_{11}-g_{11}^{free}$ the tissue elongation. Calculating this value
for arbitrarily large differences of $g_{11}$, yields figure 3(d). The most
interesting feature of it is that far from $\Delta=0$, Poisson’s ratio may be
arbitrarily large on some finite region of $\Delta$. Not surprisingly, for
$\Delta\rightarrow\infty$, $\nu\rightarrow 1$.

###  A.5 Appendix F - Relaxation of non-active tissue via $T_{1}$ transition

When only $T_{1}$ transitions are allowed, $H=-\frac{1}{2}$, we look for the
stationary solutions of

$\displaystyle\dot{q}_{i}=$
$\displaystyle-\frac{2}{\sqrt{q_{1}q_{2}}}k_{p}\left(q_{1}+q_{2}-P^{0}\right)\left[q_{i}-\frac{1}{2}\left(q_{1}+q_{2}\right)\right]q_{i}.$
(72)

Those are characterized by the initial area $A(0)=\sqrt{q_{1}(0)q_{2}(0)}$,
rather than $\,{}^{2}\\!P^{0}$.

$\displaystyle\left.\begin{array}[]{cc}q_{1}=q_{2}=A(0),&A(0)>\frac{P^{0}}{2}\\\
q_{1}=\frac{{A(0)}^{2}}{q_{2}}=\frac{1}{2}\left(P^{0}\pm\sqrt{{P^{0}}^{2}-4{A(0)}^{2}}\right),&A(0)\leq\frac{P^{0}}{2}\end{array}\right.$
(75)


# The matter Lagrangian of a non-perfect fluid

Sergio Mendoza<EMAIL_ADDRESS>Sarahí Silva<EMAIL_ADDRESS>Universidad Nacional Autónoma de México, Instituto de Astronomía, Ciudad
Universitaria AP 70-264, Ciudad de México 04510, México

(28th August 2024)

###### Abstract

We show that the matter Lagrangian of a non-perfect fluid takes the negative
or positive value of the total energy density of the fluid, which is composed
of the the rest plus internal energy densities. The ambiguity of the sign
depends on the definition of the matter action and the chosen signature.

General relativity and gravitation; Fundamental problems and general
formalism; Canonical formalism, Lagrangians and variational principles

##  I Introduction

Since the mid 1950’s it has become important to compute the value of the
matter Lagrangian of an ideal fluid [see e.g. 1, and references therein]. Over
the years, different approaches have shown that this value can be either (up
to a sign) the pressure or its full energy density, this last one composed of
the addition between the rest energy density plus the thermodynamical internal
energy density of the fluid. Mendoza and Silva [1] showed that the only
possible value that the matter Lagrangian of an ideal fluid can only be the
complete energy density of the fluid (up to a sign depending on definitions of
the matter action and the signature of the metric).

Using calculus of variations, in this article we show that the matter
Lagrangian $\mathcal{L}_{\text{matt}}$ of a non-perfect fluid depends on the
total energy density of the fluid $e$ in the form
$\mathcal{L}_{\text{matt}}=\mp e$, where the sign of the matter Lagrangian
depends on the definition of the matter action and on the signature of the
metric. In other words, the matter Lagrangian of a non-perfect fluid consists
of the rest mass energy density plus the purely thermodynamical internal
energy density of the fluid. In Section II we present the matter action and
the definition of the energy-momentum through the variations of this action.
In Section III we introduce some basic thermodynamic relations and the energy
momentum tensor of a non perfect fluid in the Eckart frame. In Section IV we
use the mass continuity equation and the energy equation in the Eckart frame
in order to find the value of the matter Lagrangian of a non-perfect fluid.
Finally in Section V we discuss our results and show that the same result for
the matter Lagrangian case is also valid in the Landau frame of the fluid.

##  II Action

The matter action $S_{\text{matt}}$ is defined as [see e.g. 2]:

$S_{\text{matt}}=\pm\frac{1}{c}\int{\mathcal{L}_{\text{matt}}\,\sqrt{-g}\,\mathrm{d}^{4}x},$
(1)

where $\mathcal{L}_{\text{matt}}$ represents the matter Lagrangian, $c$ is the
speed of light an $g$ is the non-positive determinant of the metric tensor
$g_{\alpha\beta}$. The choice of a $+$ or $-$ sign in equation (1) is a
definition convention. In this work we use a signature ($+,-,-,-$) for the
metric. Greek space-time indices vary from $0$ to $4$ and spatial Latin
indices vary from $1$ to $3$ and we use the Einstein’s summation convention
over repeated indices.

The metric, Hilbert or Belifante-Rosenfeld energy-momentum tensor
$T_{\alpha\beta}$ is defined through the variations $\delta S$ of the matter
action (1). Assuming that the matter Lagrangian is a function of the metric
tensor $g^{\alpha\beta}$ only and not of its derivatives $\partial
g^{\alpha\beta}/\partial x^{\lambda}$, the metric energy-momentum tensor is
defined as [see e.g. 2, 3, 4, 5]:

$T_{\alpha\beta}=\pm\frac{2}{\sqrt{-g}}\frac{\delta\left(\sqrt{-g}\,\mathcal{L}_{\text{matt}}\right)}{\delta
g^{\alpha\beta}},$ (2)

or, equivalently [see e.g. 1]:

$T_{\alpha\beta}=\pm 2\frac{\delta\mathcal{L}_{\text{matt}}}{\delta
g^{\alpha\beta}}\mp g_{\alpha\beta}\mathcal{L}_{\text{matt}},$ (3)

and in contravariant representation:

$T^{\alpha\beta}=\mp 2\frac{\delta\mathcal{L}_{\text{matt}}}{\delta
g_{\alpha\beta}}\mp g^{\alpha\beta}\mathcal{L}_{\text{matt}}.$ (4)

##  III Non-Perfect Fluids

The first law of thermodynamics for relativistic fluids can be written as [see
e.g. 6]:

$\mathrm{d}\left(\frac{e}{\rho}\right)=T\mathrm{d}\left(\frac{\sigma}{\rho}\right)-p\mathrm{d}\left(\frac{1}{\rho}\right),$
(5)

where $\sigma=\rho s$ is the entropy density of the fluid with $s$ the entropy
per unit mass and $\rho$ the fluid or gas (baryonic) mass density. The total
energy density $e=\rho c^{2}+\xi$ contains the rest energy density $\rho
c^{2}$ and a pure internal energy density term $\xi$. The pressure and
temperature are represented by $p$ and $T$, respectively.

In the Eckart frame [7] for relativistic fluids, the energy-momentum tensor of
a non-perfect fluid is written as [see e.g. 8]:

$T^{\alpha\beta}=T^{\alpha\beta}_{\text{pf}}+T^{\alpha\beta}_{\text{npf}},$
(6)

where $T^{\alpha\beta}_{\text{pf}}$ and $T^{\alpha\beta}_{\text{npf}}$ refer
to the perfect and non-perfect fluid contribution, respectively. Equation (6)
can be written explicitly as [see e.g. 8]:

$T^{\alpha\beta}=eu^{\alpha}u^{\beta}+(p+\Pi)h^{\alpha\beta}+q^{\alpha}u^{\beta}+q^{\beta}u^{\alpha}+\pi^{\alpha\beta},$
(7)

where $\Pi,q^{\alpha}$ and $\pi^{\alpha\beta}$ are, the viscous bulk pressure,
the heat flux and the anisotropic stress tensor respectively. In equation (7)
the projection tensor is represented by
$h^{\alpha\beta}=u^{\alpha}u^{\beta}-g^{\alpha\beta}$. The covariant
divergence of the complete energy-momentum tensor (6) is null, i.e.:

$\nabla_{\alpha}T^{\alpha\beta}=0.$ (8)

##  IV Matter Lagrangian for a Non-Perfect Fluid

In the Eckart frame the mass density current $\rho u^{\alpha}$ is parallel to
the four velocity and in the continuity equation takes the same form as for a
perfect fluid, i.e. [see e.g. 8]:

$\nabla_{\alpha}\left(\rho u^{\alpha}\right)=0,$ (9)

For this reason the variation of the mass density $\rho$ satisfies [see e.g.
1]:

$\frac{\delta\rho}{\delta
g^{\alpha\beta}}=\frac{1}{2}\rho\left(g_{\alpha\beta}-u_{\alpha}u_{\beta}\right),$
(10)

for a non-perfect fluid in the Eckart frame.

In order to obtain the variation of the entropy per unit mass $s$ in terms of
the variation of the metric tensor $g^{\alpha\beta}$ let us multiply equation
(8) by $u_{\beta}$ so that the following energy equation is obtained:

$\displaystyle
u_{\beta}\nabla_{\alpha}T^{\alpha\beta}=\nabla_{\alpha}((e+p)u^{\alpha})-u^{\alpha}\nabla_{\alpha}p+\Pi\nabla_{\alpha}u^{\alpha}$
$\displaystyle+\nabla_{\alpha}q^{\alpha}+u^{\alpha}u_{\beta}\nabla_{\alpha}q^{\beta}+u_{\beta}\nabla_{\alpha}\pi^{\alpha\beta}.$
(11)

Since the relativistic enthalpy $\omega$ is defined as:

$\omega:=e+p=\rho(\mu+Ts),$ (12)

where $\mu:=\omega/\rho-Ts$ is the relativistic chemical potential, it follows
from the first law of thermodynamics (5) that:

$\mathrm{d}\mu=\frac{1}{\rho}\mathrm{d}p-s\mathrm{d}T.$ (13)

Using this, we can rewrite the energy equation (11) as:

$\begin{split}u^{\alpha}\nabla_{\alpha}s=&-\frac{\Pi}{T\rho}\nabla_{\alpha}u^{\alpha}-\frac{1}{T\rho}\nabla_{\alpha}q^{\alpha}-\frac{1}{T\rho}u_{\beta}u^{\alpha}\nabla_{\alpha}q^{\beta}\\\
&-\frac{1}{T\rho}u_{\beta}\nabla_{\alpha}\pi^{\alpha\beta}\end{split}$ (14)

Let us now assume that the specific entropy $s$, the four velocity
$u^{\alpha}$ and thermodynamic fluxes $q^{\alpha}$ and $\pi^{\alpha\beta}$
depend only on the metric tensor $g^{\alpha\beta}$ and not on its derivatives
$\partial g^{\alpha\beta}/\partial x^{\lambda}$, so that:

$\frac{\partial}{\partial x^{\mu}}=\frac{\partial g^{\alpha\beta}}{\partial
x^{\mu}}\frac{\partial}{\partial g^{\alpha\beta}}.$ (15)

Using this operator (15) and the following relations:

$u_{\mu}\frac{\partial u^{\mu}}{\partial
g^{\alpha\beta}}=\frac{1}{2}u_{\alpha}u_{\beta},$ (16)

$\frac{\partial\sqrt{-g}}{\partial
g^{\alpha\beta}}=-\frac{1}{2}\sqrt{-g}g_{\alpha\beta},$ (17)

$\frac{\partial g^{\mu\nu}}{\partial
g^{\alpha\beta}}=\frac{1}{2}\left(\delta^{\mu}_{\alpha}\delta^{\nu}_{\beta}+\delta^{\mu}_{\beta}\delta^{\nu}_{\alpha}\right),$
(18)

$\frac{\partial g_{\mu\nu}}{\partial
g^{\alpha\beta}}=-\frac{1}{2}\left(g_{\mu\alpha}g_{\nu\beta}+g_{\mu\beta}g_{\nu\alpha}\right),$
(19)

$u_{\mu}\frac{\partial q^{\mu}}{\partial
g^{\alpha\beta}}=\frac{1}{4}\left(q_{\alpha}u_{\beta}+q_{\beta}u_{\alpha}\right),$
(20)

in equation (14) it follows that:

$\begin{split}\frac{\delta s}{\delta g^{\alpha\beta}}=\frac{\partial
s}{\partial g^{\alpha\beta}}&=\frac{1}{2}\frac{\Pi}{\rho
T}\left(g_{\alpha\beta}-u_{\alpha}u_{\beta}\right)\\\
&-\frac{1}{2}\frac{1}{\rho
T}\left(q_{\alpha}u_{\beta}+u_{\alpha}q_{\beta}+\pi_{\alpha\beta}\right),\end{split}$
(21)

for a non-perfect fluid.

Let us now choose the specific entropy $s$ and the mass density $\rho$ as the
dependent variables for any given thermodynamical quantity, so that:

$\frac{\delta}{\delta g^{\alpha\beta}}=\frac{\delta\rho}{\delta
g^{\alpha\beta}}\evaluated{\frac{\partial}{\partial\rho}}_{s}+\frac{\delta
s}{\delta g^{\alpha\beta}}\evaluated{\frac{\partial}{\partial s}}_{\rho}.$
(22)

Substitution of the operator (22) into equation (3) and using relations (10)
and (21), it follows that:

$\begin{split}&T_{\alpha\beta}=\pm\left\\{\rho\evaluated{\frac{\partial\mathcal{L}_{\text{matt}}}{\partial\rho}}_{s}+\frac{\Pi}{\rho
T}\evaluated{\frac{\partial\mathcal{L}_{\text{matt}}}{\partial
s}}_{\rho}\right\\}\left(g_{\alpha\beta}-u_{\alpha}u_{\beta}\right)\\\
&\mp\frac{1}{\rho
T}\evaluated{\frac{\partial\mathcal{L}_{\text{matt}}}{\partial
s}}_{\rho}\left(q_{\alpha}u_{\beta}+u_{\alpha}q_{\beta}+\pi_{\alpha\beta}\right)\mp
g_{\alpha\beta}\mathcal{L}_{\text{matt}},\end{split}$ (23)

for a non-perfect fluid. Direct comparison of equations (7) and (23) imply
that:

$\begin{split}&eu_{\alpha}u_{\beta}+(p+\Pi)h_{\alpha\beta}+q_{\alpha}u_{\beta}+q_{\beta}u_{\alpha}+\pi_{\alpha\beta}=\\\
&\pm\left\\{\rho\evaluated{\frac{\partial\mathcal{L}_{\text{matt}}}{\partial\rho}}_{s}+\frac{\Pi}{\rho
T}\evaluated{\frac{\partial\mathcal{L}_{\text{matt}}}{\partial
s}}_{\rho}\right\\}\left(g_{\alpha\beta}-u_{\alpha}u_{\beta}\right)\\\
&\mp\frac{1}{\rho
T}\evaluated{\frac{\partial\mathcal{L}_{\text{matt}}}{\partial
s}}_{\rho}\left(q_{\alpha}u_{\beta}+u_{\alpha}q_{\beta}+\pi_{\alpha\beta}\right)\mp
g_{\alpha\beta}\mathcal{L}_{\text{matt}}.\end{split}$ (24)

The previous equation (24) can be solved for $\mathcal{L}_{\text{matt}}$ in
two different ways:

  1. (i)

By multiplication of equation (24) by $u^{\alpha}u^{\beta}$ to obtain:

$\mathcal{L}_{\text{matt}}=\mp e.$ (25)

  2. (ii)

In order to obtain a system of equations for $\mathcal{L}_{\text{matt}}$, let
us first multiply equation (24) by the metric tensor $g^{\alpha\beta}$ to
obtain:

$e-3\left(p+\Pi\right)=\pm
3\rho\evaluated{\frac{\partial\mathcal{L}_{\text{matt}}}{\partial\rho}}_{s}\pm
3\frac{\Pi}{\rho
T}\evaluated{\frac{\partial\mathcal{L}_{\text{matt}}}{\partial s}}_{\rho}\mp
4\mathcal{L}_{\text{matt}}.$ (26)

and multiply equation (24) by the projection tensor
$h^{\alpha}{}_{\lambda}:=u^{\alpha}u_{\lambda}-\delta^{\alpha}{}_{\lambda}$ to
obtain:

$\begin{split}&-(p+\Pi)(u_{\lambda}u_{\beta}-g_{\lambda\beta})-(q_{\lambda}u_{\beta}+\pi_{\lambda\beta})=\\\
&\mp\left(\mathcal{L}_{\text{matt}}-\rho\evaluated{\frac{\partial\mathcal{L}_{\text{matt}}}{\partial\rho}}_{s}-\frac{\Pi}{\rho
T}\evaluated{\frac{\partial\mathcal{L}_{\text{matt}}}{\partial
s}}_{\rho}\right)(u_{\lambda}u_{\beta}-g_{\lambda\beta})\\\ &\pm\frac{1}{\rho
T}\evaluated{\frac{\partial\mathcal{L}_{\text{matt}}}{\partial
s}}_{\rho}\left(q_{\lambda}u_{\beta}+\pi_{\lambda\beta}\right),\end{split}$
(27)

Equating terms with $(u_{\lambda}u_{\beta}-g_{\lambda\beta)}$ and the ones
with $(q_{\lambda}u_{\beta}+\pi_{\lambda\beta})$ in equation (27), the
following two relations are obtained:

$\displaystyle
p+\Pi=\pm\left(\mathcal{L}_{\text{matt}}-\rho\evaluated{\frac{\partial\mathcal{L}_{\text{matt}}}{\partial\rho}}_{s}-\frac{\Pi}{\rho
T}\evaluated{\frac{\partial\mathcal{L}_{\text{matt}}}{\partial
s}}_{\rho}\right),$ (28) and $\displaystyle\frac{1}{\rho
T}\evaluated{\frac{\partial\mathcal{L}_{\text{matt}}}{\partial s}}_{\rho}=\mp
1.$ (29)

The solution to the system of equations (26), (28) and (29) recovers the same
value of $\mathcal{L}_{\text{matt}}=\mp e$ given in relation (25).

With the value of $\mathcal{L}_{\text{matt}}=\mp e$ the set of equations
(28)-(29) can be written as:

$\evaluated{\frac{\partial e}{\partial\rho}}_{s}=\frac{e+p}{\rho},$ (30)

and

$\evaluated{\frac{\partial e}{\partial s}}_{\rho}=\rho T,$ (31)

which are consistent with the first law of thermodynamics (5).

##  V Final remarks

In this work we have shown that the value of the matter Lagrangian for a non-
perfect fluid is

$\mathcal{L}_{\text{matt}}=\mp e,$ (32)

where the sing depends on the definition of the matter action and on the
signature of the metric111For the choice of a ($-,+,+,+$) signature, equation
(10) turns into [see e.g. 1]:  $\frac{\delta\rho}{\delta
g^{\alpha\beta}}=\frac{1}{2}\rho\left(g_{\alpha\beta}+u_{\alpha}u_{\beta}\right),$
and equation (21) is now:  $\frac{\delta s}{\delta
g^{\alpha\beta}}=\frac{1}{2}\frac{\Pi}{\rho
T}\left(g_{\alpha\beta}+u_{\alpha}u_{\beta}\right)+\frac{1}{2}\frac{1}{\rho
T}\left(q_{\alpha}u_{\beta}+u_{\alpha}q_{\beta}+\pi_{\alpha\beta}\right).$
(33) With all these results it follows that:  $\mathcal{L}_{\text{matt}}=\pm
e.$ (34) 222In the Landau frame the energy momentum tensor $T_{\alpha\beta}$
and the particle flux density $\rho_{\alpha}$ of a non-perfect fluid are
respectively given by [see e.g. 6]:
$T_{\alpha\beta}=(e+p)u_{\alpha}u_{\beta}-pg_{\alpha\beta}+\tau_{\alpha\beta},$
(35) $\rho_{\alpha}=\rho u_{\alpha}+\nu_{\alpha},$ (36) where
$\tau_{\alpha\beta}$ is the viscosity tensor and $\nu_{\alpha}$ is the thermal
conduction of the fluid. Analogously to the development shown in section IV,
the variation of the mass density $\rho$ and the entropy per unit mass $s$ can
be written as:  $\frac{\delta\rho}{\delta
g^{\alpha\beta}}=\frac{1}{2}\rho(g_{\alpha\beta}-u_{\alpha}u_{\beta})-\frac{1}{4}(\nu_{\alpha}u_{\beta}+u_{\alpha}\nu_{\beta}),$
(37) and  $\frac{\delta s}{\delta
g^{\alpha\beta}}=\frac{1}{4\rho}\left(s+\frac{\mu}{T}\right)(\nu_{\alpha}u_{\beta}+u_{\alpha}\nu_{\beta})-\frac{1}{2\rho
T}\tau_{\alpha\beta},$ (38) respectively. With all this, the same result of
equation (32) is also obtained in the Landau frame.. This value is consistent
with the first law of thermodynamics (5).

Note that the same result for the matter Lagrangian is obtained for an ideal
fluid [see e.g. 1]. The difference resides in the value of the total energy
density for each fluid, which is the addition between the rest mass energy
density plus the pure internal energy density.

## Acknowledgements

This work was supported by a PAPIIT DGAPA-UNAM grant IN110522. SM and SS
acknowledge support from CONAHCyT (26344,1228643).

## References

  * Mendoza and Silva [2021] S. Mendoza and S. Silva, International Journal of Geometric Methods in Modern Physics 18, 2150059 (2021).
  * Landau [2013] L. Landau,  _The Classical Theory of Fields_ , Course of Theoretical Physics (Elsevier Science, 2013).
  * Harko and Lobo [2018] T. Harko and F. Lobo,  _Extensions of f(R) Gravity_ , Cambridge Monographs on Mathem (Cambridge University Press, 2018).
  * Franklin [2017] J. Franklin,  _Classical Field Theory_ (Cambridge University Press, 2017).
  * Năstase [2019] H. Năstase,  _Classical Field Theory_ (Cambridge University Press, 2019).
  * Landau and Lifshitz [2013] L. Landau and E. Lifshitz,  _Fluid Mechanics_ , Course of Theoretical Physics No. v. 6 (Elsevier Science, 2013).
  * Eckart [1948] C. Eckart, Physical Review 73, 373 (1948).
  * Rezzolla and Zanotti [2013] L. Rezzolla and O. Zanotti,  _Relativistic hydrodynamics_ (Oxford University Press, USA, 2013).