A

a) b) Figure 1. a) Exemplary results of Neural Volumetric Blendshapes (brown)
compared to linear blendshapes (blue). Among others, volume preservation,
wrinkles, and avoided self-intersections result in more realistic and
immersive facial animations. b) An example fit of the novel layered head
model. The model encapsulates the skin, the skull, and the muscles with wraps
for which we present a data-driven fitting algorithm. The space between the
wraps can be canonically tetrahedralized.

# Neural Volumetric Blendshapes: Computationally Efficient Physics-Based
Facial Blendshapes

Nicolas Wagner  ,  Ulrich Schwanecke  and  Mario Botsch

###### Abstract.

Computationally weak systems and demanding graphical applications are still
mostly dependent on linear blendshapes for facial animations. The accompanying
artifacts such as self-intersections, loss of volume, or missing soft tissue
elasticity can be avoided by using physics-based animation models. However,
these are cumbersome to implement and require immense computational effort. We
propose neural volumetric blendshapes, an approach that combines the
advantages of physics-based simulations with realtime runtimes even on
consumer-grade CPUs. To this end, we present a neural network that efficiently
approximates the involved volumetric simulations and generalizes across human
identities as well as facial expressions. Our approach can be used on top of
any linear blendshape system and, hence, can be deployed straightforwardly.
Furthermore, it only requires a single neutral face mesh as input in the
minimal setting. Along with the design of the network, we introduce a pipeline
for the challenging creation of anatomically and physically plausible training
data. Part of the pipeline is a novel hybrid regressor that densely positions
a skull within a skin surface while avoiding intersections. The fidelity of
all parts of the data generation pipeline as well as the accuracy and
efficiency of the network are evaluated in this work. Upon publication, the
trained models and associated code will be released.

††journal: TOG

##  1\. Introduction

At present, research in the field of head avatars and facial animation is
mainly concerned with obtaining photorealistic results through neural networks
(Cao et al., 2022; Athar et al., 2022) which can be operated on
computationally rich systems. What usually falls short, however, is the
inclusion of less capable hardware setups that mostly rely on the use of
linear blendshapes (Lewis et al., 2014). Although linear blendshapes have been
intensively researched and improved over the past decades, there are still
known shortcomings like physically implausible distortions, loss of volume,
anatomically illegal expressions, missing volumetric elasticity, or self-
intersections. Physics-based simulations have been proposed that overcome most
artifacts of linear blendshapes (Ichim et al., 2017, 2016; Cong, 2016; Choi et
al., 2022), nonetheless, they are usually laborious to handle and
computationally expensive, too. Hybrid models that try to combine the best of
both worlds are either not sophisticated enough in the quality of the
simulated physical properties (Barrielle et al., 2016) or are still too
inefficient to be used on slower devices (Ichim et al., 2016).

An auspicious approach in the latter category is physics-based volumetric
blendshapes (Ichim et al., 2016). These can be animated with physical
plausibility, anatomical constraints can be taken into account, self-
intersections can be prevented, and the control is identical to linear
blendshapes. Although the level of detail is slightly sacrificed in comparison
to other physics simulations (Cong, 2016; Ichim et al., 2017), volumetric
blendshapes can still only be used at low frame rates. We improve on this
approach with _neural volumetric blendshapes_ that approximate the involved
physics-based simulation with a neural network. At this, the efficient network
design proposed in this work allows for inference on consumer-grade hardware.
Further, neural volumetric blendshapes include an adapted interpolation scheme
which is truly volume-preserving. A biologically necessary constraints which
previously has been ignored. The network is paired with a training data
generation pipeline which produces training instances that reflect the
previously mentioned advantages of physics-based simulations and facilitates a
generalization of the network across different human identities. Thus, in
contrast to other recent work that tries to approximate physics-based facial
simulations by neural networks (Srinivasan et al., 2021; Yang et al., 2022;
Choi et al., 2022), we avoid time-consuming individualizations for
starightforward deployment. Also, our animation system does not require
multiple or even sequences of optical scans. Instead, the network trained in
this work induces an anatomical plausible deformation transfer such that our
system is instantly applicable to any neutral head surface or on top of
arbitrary linear blendshape rigs.

The aforementioned data generation pipeline is unique in that, to the best of
our knowledge, there are no other comprehensive datasets that relate a broad
range of head shapes in diverse facial expressions to the underlying
anatomical and physical characteristics. To curate such a dataset for the
first time, we bring together multiple data sources such as CT data to reflect
the anatomy of heads, 3D reconstructions of images in the wild to collect
diverse head shapes, or facial expressions in the form of recorded blendshape
weights from dyadic conversational situations. The result is a dataset of
heads with neutral and non-neutral expressions represented by a novel
standardized model, encompassing not only the skin surface but also the
biomechanical anatomy.

The _key novelties and contributions_ we present in this paper can be
summarized as follows:

  * •

A novel _layered head model_ (LHM) representing the skin surface as well as
the entire biomechanical anatomy.

  * •

A data-driven procedure for _fitting the LHM_ to neutral skin surfaces.

  * •

An _inverse physics-based simulation_ which fits the LHM to skin surfaces of
non-neutral facial expressions.

  * •

A novel neural network design that approximates physics-based simulations to
efficiently implement neural volumetric blendshapes.

  * •

A pipeline for _creating training data_ of the neural network that among other
parts includes the LHM and the associated fitting procedures.

##  2\. Related Work

###  2.1. Personalized Anatomic Models

Algorithms that create personalized anatomical models can essentially be
distinguished according to two paradigms: _heuristic-based_ and _data-driven_.
Considering heuristic-based approaches, Anatomy Transfer (Ali-Hamadi et al.,
2013) applies a space warp to a template anatomical structure to fit a target
skin surface. The skull and other bones are only deformed by an affine
transformation. A similar idea is proposed by Gilles et al. (Gilles et al.,
2010). While they also implement a statistical validation of bone shapes, the
statistics are collected from artificially deformed bones. In (Ichim et al.,
2016; Kadleček et al., 2016), an inverse physics simulation was used to
reconstruct anatomical structures from multiple 3D expression scans. Saito et
al. (Saito et al., 2015) simulate the growth of soft tissue, muscles, and
bones. A musculoskeletal biomechanical model is fitted from sparse
measurements in (Schleicher et al., 2021) but not qualitatively evaluated.

There are only a few data-driven approaches because combined data sets of
surface scans and MRI, DXA, or CT images are hard to obtain for various
reasons (e.g. data privacy or unnecessary radiation exposure). The recent work
of OSSO (Keller et al., 2022) predicts full body skeletons from 2000 DXA
images that do not carry precise 3D information. Further, bones are positioned
within a body by predicting only three anchor points per bone group and not
avoiding intersections between skin and skull. A model that prevents skin-
skull intersections and also considers muscles is based on fitting
encapsulating wraps instead of the anatomy itself (Komaritzan et al., 2021).
However, no accurate algorithm based on medical imaging but a BMI regressor
(Maalin et al., 2021) is used to position the wraps. A much more accurate,
pure face model, was developed by Achenbach et al. (Achenbach et al., 2018).
Here, CT scans are combined with optical scans by a multilinear model (MLM)
which can map from skulls to faces and vice versa. As before, no self-
intersections are prevented and only bones are fitted. Building on the data
from (Achenbach et al., 2018) and following the idea of a layered body model
(Komaritzan et al., 2021), we create a statistical layered head model
including musculature that avoids self-intersections.

###  2.2. Physics-Based Facial Animation

A variety of techniques for animating faces have been developed in the past
(Ichim et al., 2015; Bradley et al., 2010; Zhang et al., 2008; Parke, 1991).
Data-driven models (Lewis et al., 2014; Ichim et al., 2016; Lewis et al.,
2005), which have recently been significantly improved by deep learning (Zheng
et al., 2022; Garbin et al., 2022; Cao et al., 2022; Feng et al., 2021; Song
et al., 2020; Athar et al., 2022), are certainly dominant. Due to their
simplicity and speed, especially linear blendshapes (Lewis et al., 2014) are
still used most often in demanding applications and whenever no
computationally rich hardware is available. Physics-based models have been
developed for a long time (Sifakis et al., 2005) and avoid artifacts like
implausible contortions and self-intersections, but due to their complexity
and computational effort, they are rarely used. Hybrid approaches add surface-
based physics to linear blendshapes for more detailed facial expressions.
(Barrielle et al., 2016; Bickel et al., 2008). However, by construction they
can not model volumetric effects.

The pioneering work of Sifakis et al. (Sifakis et al., 2005) is the first
fully physics-based facial animation. The simulation is conducted on a
personalized tetrahedron mesh which can only be of a limited resolution due to
a necessary dense optimization problem. With Phace (Ichim et al., 2017), this
problem was overcome by an improved physics simulation. The art-directed
muscle model (Cong and Fedkiw, 2019; Bao et al., 2019; Cong, 2016)
additionally represents muscles as B-splines and allows control of expressions
via trajectories of spline control points. A solely inverse model for
determining the physical properties of faces was proposed in (Kadleček and
Kavan, 2019). Neural soft-tissue dynamics (Santesteban et al., 2020; Casas and
Otaduy, 2018) extend the SMPL (Loper et al., 2015) body model with secondary
motion. Recently, (Srinivasan et al., 2021; Yang et al., 2022; Choi et al.,
2022) adapted neural soft-body dynamics to learn the physical properties of a
particular person. However, these approaches must be retrained for new objects
and are slow in inference. With volumetric blendshapes (Ichim et al., 2016), a
hybrid approach has been presented that combines the structure of linear
blendshapes with physical and anatomical plausibility. We extend this work to
neural volumetric blendshapes to make physical plausibility realtime capable
while maintaining the control structure of standard linear blendshapes.

##  3\. Method

The basis of our neural volumetric blendshapes approach (Section 3.2) is a
novel layered head representation (Section 2). Starting from there, we design
a physics-based facial animation model (Section 3.3) and distill it into a
defining dataset (Section 3.4). With this dataset, a neural network that makes
the animation model realtime capable is trained (Section 3.5).

###  3.1. Layered Head Model

Figure 2. All components of the layered head model template $\mathcal{T}$.
Skin $S_{\mathcal{T}}$, skin wrap $\hat{S}_{\mathcal{T}}$, muscles
$M_{\mathcal{T}}$, muscles wrap $\hat{M}_{\mathcal{T}}$, skull
$B_{\mathcal{T}}$, and the skull wrap $\hat{B}_{\mathcal{T}}$.

We represent a head $\mathcal{H}=\rho_{\mathcal{H}}(\mathcal{T})$ with neutral
expression through a component-wise transformation (see Section 3.4 for
details) of a layered head model (LHM) template

(1)
$\mathcal{T}=(S_{\mathcal{T}},B_{\mathcal{T}},M_{\mathcal{T}},\hat{S}_{\mathcal{T}},\hat{B}_{\mathcal{T}},\hat{M}_{\mathcal{T}}),$

that consists of six triangle meshes. $S_{\mathcal{T}}$ describes the skin
surface including the eyes, the mouth cavity, and the tongue,
$B_{\mathcal{T}}$ the surface of all skull bones and teeth, $M_{\mathcal{T}}$
the surface of all muscles and the cartilages of the ears and nose.
$\hat{S}_{\mathcal{T}}$ is the skin layer, i.e. a closed wrap enveloping
$S_{\mathcal{T}}$, $\hat{B}_{\mathcal{T}}$ the skull layer that envelopes
$B_{\mathcal{T}}$, and $\hat{M}_{\mathcal{T}}$ the muscle layer that envelopes
$M_{\mathcal{T}}$. Other anatomical structures are omitted for simplicity. The
template structures $S_{\mathcal{T}},B_{\mathcal{T}}$ and $T_{\mathcal{T}}$
were designed by an experienced digital artist. The skin, skull, and muscle
layers $\hat{S}_{\mathcal{T}},\hat{B}_{\mathcal{T}}$, and
$\hat{M}_{\mathcal{T}}$ were generated by shrink-wrapping a sphere as close as
possible to the corresponding surfaces without intersections and follow the
same triangulation. The complete template is shown in Figure 2.

The skull and muscle layers are further massaged such that the quads of all
prisms that can be spanned between corresponding faces of the skin, muscle,
and skull layer are as rectangular as possible while preserving the original
shapes. To that end, we determine for the skull layer

(2) $\displaystyle\arg\min_{X}$ $\displaystyle
w_{r}{E}_{r}(X,\hat{B}_{\mathcal{T}})+w_{d}E_{d}(X,\hat{B}_{\mathcal{T}}),$

within the efficient projective dynamics (Bouaziz et al., 2014) optimization
framework, initializing with $X=\hat{B}_{\mathcal{T}}$. Here, $E_{d}$ is the
two-sided Hausdorff distance to the non-massaged shape and

(3) $\displaystyle E_{r}(X,\hat{B}_{\mathcal{T}})=$
$\displaystyle\sum_{(x_{i}^{0},x_{i}^{1})\in X}(90^{\circ}-\angle
x_{i}^{0}x_{i}^{1}\hat{b}_{i}^{1})^{2}+(90^{\circ}-\angle
x_{i}^{1}x_{i}^{0}\hat{b}_{i}^{0})^{2}$
$\displaystyle\quad\quad\quad+(90^{\circ}-\angle\hat{b}_{i}^{0}\hat{b}_{i}^{1}x_{i}^{1})^{2}+(90^{\circ}-\angle\hat{b}_{i}^{1}\hat{b}_{i}^{0}x_{i}^{0})^{2},$

induces the rectangular prism shapes. After the optimization we set
$\hat{B}_{\mathcal{T}}=X$. The same optimization is run for the muscle layer.

Mesh | $S_{\mathcal{T}}$ | $B_{\mathcal{T}}$ | $M_{\mathcal{T}}$ | $\hat{S}_{\mathcal{T}}$
---|---|---|---|---
#Vertices | 21875 | 14572 | 16388 | 29039
#Faces / #Tets | 42738 | 28856 | 32370 | 58056
|  |  |  |
Mesh | $\hat{B}_{\mathcal{T}}$ | $\hat{M}_{\mathcal{T}}$ | $\mathbb{S}_{\mathcal{T}}$ | $\mathbb{M}_{\mathcal{T}}$
#Vertices | 29039 | 29039 | 49852
#Faces / #Tets | 58056 | 58056 | 123429 | 73681
Table 1. An overview of the sizes of the individual template LHM $\mathcal{T}$
components.

The wrapping layers of the LHM allow for at least two significant advantages.
On the one hand, they provide a simplified representation of the skin surface,
the musculature, and the skull, which we exploit in Section 3.4 for
determining $\rho_{\mathcal{H}}$. On the other hand, they can be used for
topologically and semantically consistent tetrahedralization of the head
volume by splitting the prisms between the layers into tetrahedrons.

More precisely, the construction of the LHM also defines a soft tissue tet
mesh $\mathbb{S}_{\mathcal{T}}$ (i.e. between the skin and the muscle layer)
and a muscle tissue tet mesh $\mathbb{M}_{\mathcal{T}}$ (i.e. between the
muscle and the skull layer). The massage of the muscle and skull layers
ensures nicely shaped tets with minimal shearing. Further,
$\mathbb{S}_{\mathcal{T}}$ can be fine-tuned by removing all vertices and
connected tets outside of $S_{\mathcal{T}}$, adding the vertices of
$S_{\mathcal{T}}$, and connecting them via Delaunay tetrahedralization. The
sizes of all template components are given in Table 1.

###  3.2. Neural Volumetric Blendshapes

Building on the LHM representation, we can now derive neural volumetric
blendshapes. For this, the classical concept of linear blendshapes is reviewed
first. Thereupon, we derive volumetric blendshapes as well as physics-based
facial animation and elaborate on the accompanying problems. Finally, we
introduce neural volumetric as a solution to the previously identified
shortcomings.

For a specific head
$\mathcal{H}=(S_{\mathcal{H}},B_{\mathcal{H}},M_{\mathcal{H}},\hat{S}_{\mathcal{H}},\hat{B}_{\mathcal{H}},\hat{M}_{\mathcal{H}})$,
a linear blendshape rig of $n$ blendshapes is a set of expression surfaces

(4) $\\{{S}_{\mathcal{H}}^{i}\\}_{i=1}^{n}$

which often includes the neutral surface $S_{\mathcal{H}}$. A novel expression
is formed through linear interpolation as

(5) $S_{\mathcal{H}}^{exp}=\sum_{i}w_{exp}^{i}S_{\mathcal{H}}^{i},$

where the $w_{exp}^{i}$ weight the share of each blendshape in the expression.

The corresponding volumetric blendshape rig can be defined as

(6)
$\\{V^{i}_{\mathcal{H}}\\}_{i=1}^{n}=\\{(\nabla\mathbb{S}_{\mathcal{H}}^{i},\nabla\mathbb{M}_{\mathcal{H}}^{i},B_{\mathcal{H}}^{i})\\}_{i=1}^{n}.$

At this, $\nabla\mathbb{S}_{\mathcal{H}}^{i}$,
$\nabla\mathbb{M}_{\mathcal{H}}^{i}$ describe soft and muscle tissue
deformations as vectors of $3\times 3$ deformation gradients.
$B_{\mathcal{H}}^{i}$ describes the skull after rigid motion of jaw and
cranium. An anatomical plausible inverse physics model $\phi^{\dagger}$ (see
Section 3.3 for details) relates the linear to the volumetric rig as

(7)
$V^{i}_{\mathcal{H}}=\phi^{\dagger}(S_{\mathcal{H}}^{i},\mathbb{S}_{\mathcal{H}},\mathbb{M}_{\mathcal{H}},B_{\mathcal{H}}).$

Expressed in words, $\phi^{\dagger}$ determines the volumetric deformations of
a head that cause the surface deformation of a facial expression. Vice versa,
a forward physics model $\phi$ acts as an left inverse to $\phi^{\dagger}$ and
maps volumetric deformations back to surface deformations as

(8)
$\bar{S}_{\mathcal{H}}^{i}=\phi(\mathbb{S}_{\mathcal{H}},\mathbb{M}_{\mathcal{H}},V^{i}_{\mathcal{H}}).$

In reality, ${S}_{\mathcal{H}}^{i}$ is often times not an anatomically legal
facial expression with respect to $\phi^{\dagger}$. As a consequence,
$\phi^{\dagger}$ is build such that $\bar{S}_{\mathcal{H}}^{i}$ is an
anatomical reachable expression which, in an euclidean sense, is close but not
necessarily equal to ${S}_{\mathcal{H}}^{i}$.

The major advantage of volumetric blendshapes rigs is that novel expressions
calculated as

(9)
$S_{\mathcal{H}}^{exp}=\phi(\mathbb{S}_{\mathcal{H}},\mathbb{M}_{\mathcal{H}},\bigoplus_{i}w_{exp}^{i}\otimes
V^{i}_{\mathcal{H}})$

can be guaranteed to be anatomically plausible through a suitable choice of
non-linear interpolation functions $\bigoplus$ and $\otimes$. Adaptions of
volumetric blendshape rigs and Equation 9 are the foundation of many state-of-
the-art physics-based facial animation systems (Ichim et al., 2017; Cong and
Fedkiw, 2019). The common approach of these systems to non-linear
interpolation is to separately interpolate the stretch and the rotation
components of the deformation gradients to ensure the biologically necessary
volume-preservation. However, considering the stretching components, the more
complex volume-preserving interpolation (Arsigny et al., 2007; Jung et al.,
2015) is mostly ignored and only a linear interpolation is carried out.
Considering the rotation components, it is has not been discussed, yet, how
these should be extrapolated for facial animations if blend weights do not
form a convex combination.

In our approach to improve on volumetric blendshapes, we therefore propose a
hybrid approach for creating novel expressions as

(10)
$S_{\mathcal{H}}^{exp}=\phi(\phi^{\dagger}(\sum_{i}w_{exp}^{i}\phi(V_{\mathcal{H}}^{i}),\mathbb{S}_{\mathcal{H}},\mathbb{M}_{\mathcal{H}},B_{\mathcal{H}})).$

By doing so, the straightforward and fast interpolation of linear blendshapes
can be used while maintaining the advantages of anatomical and physical
plausibility. What remains is the issue that high-quality physics-based
animations mostly have in common, i.e. slow runtimes. As a countermeasure, we
propose neural volumetric blendshapes. Additionally to the redesigned
expression evaluation of Equation 10, neural volumetric blendshapes use an
efficient neural network $f$ that approximates $\phi(\phi^{\dagger}(\cdot))$.

In the following, we first depict the construction of $\phi$ and
$\phi^{\dagger}$, subsequently structure a pipeline for creating a dataset
that reflects the physics-models, and finally design and learn $f$.

a) b) Figure 3. a) The training scheme of the skin to skull layer distances
regressor which is used in the layered head model fitting. b) Procedural
overview of the layered head model fitting algorithm. Orange frames indicate
inputs, blue frames outputs.

###  3.3. Physics-Based Models

We realize the anatomically plausible inverse physics $\phi^{\dagger}$ and the
left inverse $\phi$ as projective dynamics energies $E_{\phi^{\dagger}}$ and
$E_{\phi}$, respectively. Starting with $\phi^{\dagger}$, separate terms for
soft tissue, muscle tissue, the skin, the skull, and auxiliary components are
applied. Considering the soft tissue $\mathbb{S}_{\mathcal{H}}$, we closely
follow the model of (Ichim et al., 2017) and impose

(11)
$E_{{\mathbb{S}_{\mathcal{H}}}}=w_{vol}\sum_{t\in\mathbb{S}_{\mathcal{H}}}E_{vol}(t)+w_{soft}\sum_{t\in\mathbb{S}_{\mathcal{H}}}E_{str}(t),$

which for a tet $t$ changes in volume

(12) $\displaystyle E_{vol}(t)=(\det(\text{F}(t))-1)^{2}$

and strain

(13) $\displaystyle E_{str}(t)=\min_{R\in
SO(3)}\lvert\lvert\text{F}(t)-R\rvert\rvert^{2}_{F}$

are penalized. $\text{F}(t)$ calculates the deformation gradient of a tet $t$.
The matrix $R\in SO(3)$ denotes the optimal rotation given the Frobenius norm.

To reflect the biological structure of the skin, we additionally impose a
dedicated strain constraint

(14) $E_{{S_{\mathcal{H}}}}=\sum_{t\in{S}_{\mathcal{H}}}E_{str}(t)$

on the triangles of the skin.

For the muscle tets $\mathbb{M}_{\mathcal{H}}$ we follow the arguments of
(Kadleček and Kavan, 2019) that capturing fiber directions for tetrahedralized
muscles is in general too restrictive. Hence, only a volume-preservation term

(15)
$E_{{\mathbb{M}_{\mathcal{H}}}}=w_{vol}\sum_{t\in\mathbb{M}_{\mathcal{H}}}E_{vol}(t)$

is applied. The weight $w_{vol}$ is set sufficiently high such that the
volume-preservation is almost a hard constraint.

The skull is not tetrahedralized as it is assumed to be non-deformable even
though it is rigidly movable. The non-deformability of the skull is
represented by

(16) $E_{B_{\mathcal{H}}}\\!=w_{skull}\\!\\!\sum_{t\in
B_{\mathcal{H}}}\\!E_{str}(t)+w_{c}\\!\\!\sum_{x\in
B_{\mathcal{H}}}\\!E_{c}(x,B_{\mathcal{H}}),$

i.e. a strain $E_{str}$ on the triangles and mean curvature regularization

(17) $E_{c}(x,B_{\mathcal{H}})=A\lvert\lvert\Delta x-R\Delta
b\rvert\rvert^{2}$

on the vertices of the skull. The matrix $R\in SO(3)$ denotes the optimal
rotation keeping the vertex Laplacian $\Delta x$ as close as possible to its
initial value $\Delta b$. The vertex Laplacian is discretized using the
cotangent weights and the Voronoi areas $A$ (Botsch et al., 2010). Like
$w_{vol}$, the weights $w_{skull},w_{c}$ are set such that the non-
deformability is almost a hard constraint, too.

To connect the muscle tets as well as the eyes to the skull, connecting tets
are introduced. For the muscle tets, each skull vertex connects to the closest
three vertices in $\mathbb{M}_{\mathcal{H}}$ to form a connecting tet. For the
eyes, connecting tets are formed by connecting each eye vertex to the three
closest vertices in $B_{\mathcal{H}}$. On these connecting tets, the energy
$E_{con}$ with the same constraints as in Equation 11 is imposed as almost
hard constraints. By this design, the jaw and the cranium are moved
independent from each other by muscle activation but the eyes remain rigid and
move only with the cranium.

Finally, the energy

(18) $E_{inv}=\sum_{x\in S_{\mathcal{H}}}E_{tar}(x,{S}^{exp}_{\mathcal{H}})$

of soft Dirichlet constraints

(19) $\displaystyle E_{tar}(x,{S}^{exp}_{\mathcal{H}})=\lvert\lvert
x-{s}\rvert\rvert^{2},$

is added, attracting each vertex $x$ of the skin surface $S_{\mathcal{H}}$ of
the simulated head to a unique vertex $s$ from the target expression
${S}^{exp}_{\mathcal{H}}$. As previously mentioned, the expression
${S}^{exp}_{\mathcal{H}}$ might be anatomically illegal. As a countermeasure,
we can impose a maximum strain by balancing $w_{soft}$ with $w_{tar}$. Thus,
together with the almost hard constraints and by the construction of
projective dynamics, $\phi^{\dagger}$ always results in a plausible expression
$\bar{S}^{exp}_{\mathcal{H}}$ close to the target one. The weighted sum of the
aforementioned energies gives the total energy

(20) $\displaystyle E_{\phi^{\dagger}}=$
$\displaystyle\>\>w_{{\mathbb{S}_{\mathcal{H}}}}E_{{\mathbb{S}_{\mathcal{H}}}}+w_{{\mathbb{M}_{\mathcal{H}}}}E_{{\mathbb{M}_{\mathcal{H}}}}+w_{{B_{\mathcal{H}}}}E_{{B_{\mathcal{H}}}}$
$\displaystyle+w_{{{S}_{\mathcal{H}}}}E_{{{S}_{\mathcal{H}}}}+w_{con}E_{con}+w_{inv}E_{inv}$

of the backward model $\phi^{\dagger}$.

The forward model $\phi$ is considerably simpler in structure and is realized
as the energy

(21) $\displaystyle E_{\phi}=$
$\displaystyle{}\sum_{t\in\mathbb{S}_{\mathcal{H}}}E_{dg}(t,\nabla\mathbb{S}_{\mathcal{H}}^{exp})+\sum_{t\in\mathbb{M}_{\mathcal{H}}}E_{dg}(t,\nabla\mathbb{M}_{\mathcal{H}}^{exp})$
$\displaystyle+w_{tar}\sum_{x\in
B_{\mathcal{H}}}E_{tar}(x,{B}^{exp}_{\mathcal{H}})$

that aims to achieve given deformation gradients
$(\nabla\mathbb{S}_{\mathcal{H}}^{exp},\nabla\mathbb{M}_{\mathcal{H}}^{exp})$
through

(22) $\displaystyle E_{dg}(t)=\min_{R\in
SO(3)}\lvert\lvert\text{F}(t)-R\rvert\rvert^{2}_{F}.$

while matching the vertices of the targeted skull ${B}^{exp}_{\mathcal{H}}$.

In both models we resolve self-intersections between colliding lips or teeth
in a subsequent projective dynamics update. In the second update, colliding
vertices are resolved as in (Komaritzan and Botsch, 2018). The distinctive
feature here is that no collision-gaps can occur after dissolving self-
intersections.

###  3.4. Generation of Training Data

To approximate $\phi(\phi^{\dagger}(\cdot))$ by a neural network $f$, a
defining training dataset is required in first place. By the construction of
Equation 10, this training dataset must consist of instances that relate
diverse facial expressions created through linear blendshapes to the
corresponding $\phi(\phi^{\dagger}(\cdot))$-generated surfaces. This dataset
must also cover a variety of distinct head shapes to train $f$ to be as
generally applicable as $\phi(\phi^{\dagger}(\cdot))$.

In the following, we describe a pipeline for creating such training instances
which can be roughly divided into two high-level steps. First, in order to
evaluate $\phi(\phi^{\dagger}(\cdot))$ on a reasonable domain, it is necessary
to model a head $\mathcal{H}$ from an extensive head-model and determine the
corresponding neutral soft and muscle tissue tet meshes
$\mathbb{S}_{\mathcal{H}}$ and $\mathbb{M}_{\mathcal{H}}$. Second,
$S_{\mathcal{H}}$ has to be deformed to an expression
${S}_{\mathcal{H}}^{exp}$ and mapped to the anatomically and physically
plausible $\bar{S}_{\mathcal{H}}^{exp}$. A more detailed algorithmic
description is given in Algorithm 1.

Sampling Head Shapes We start off the first part of the pipeline by randomly
drawing a neutral skin surface $S_{\mathcal{H}}$ from DECA (Feng et al.,
2021), one of the most comprehensive high-resolution face models currently
available. More specifically, we randomly draw an image from the Flickr-Faces-
HQ (Karras et al., 2019) dataset and let DECA determine the corresponding
neutral head shape. Further, a precomputed mapping is applied to adapt the
DECA topology to our template.

Figure 4. The MLP architecture that realizes $f$ and approximates
$\phi(\phi^{\dagger}(\cdot))$. At this, $f$ is designed as an adapter such
that the optional part left of the dotted line can be exchanged. Algorithm 1
Data Generation
---
1\. a) |  Draw a random image $I_{\mathcal{H}}$ of a head $\mathcal{H}$ from the
FlickerHQ dataset.
1. b) |  Calculate skin surface $S_{\mathcal{H}}=\mathit{DECA}(I)$ with neutral expression parameters.
1. c) |  Find LHM transformation $\rho_{\mathcal{H}}$ from $S_{\mathcal{H}}$, build up tet meshes $\mathbb{S}_{\mathcal{H}}$ and $\mathbb{M}_{\mathcal{H}}$.
2\. a) |  Create ARKit blendshapes $\\{S_{\mathcal{H}}^{i}\\}_{i=1}^{52}$ from $S_{\mathcal{H}}$ with deformation transfer.
2. b) |  Sample weights $w_{exp}=\\{w^{i}_{exp}\\}_{i=1}^{52}$ from dyadic recordings and calculate $S_{\mathcal{H}}^{exp}=\sum_{i=1}^{52}w^{i}_{exp}S_{\mathcal{H}}^{i}$.
1. c) |  Get $\bar{S}_{\mathcal{H}}^{exp}$ from $\phi(\phi^{\dagger}(S_{\mathcal{H}}^{exp},\mathbb{S}_{\mathcal{H}},\mathbb{M}_{\mathcal{H}},B_{\mathcal{H}}))$.

Fitting the LHM Next, the template LHM $\mathcal{T}$ is aligned with the skin
surface $S_{\mathcal{H}}$ by finding $\rho_{\mathcal{H}}$ that maps each of
the template components individually. For this, we rely on a hybrid approach
that is largely data-driven but also based on heuristics that ensure anatomic
plausibility and avoids intersections.

As the first of the remaining five template meshes, we fit the skin layer by
setting

(23)
$\hat{S}_{\mathcal{H}}=\text{rbf}_{S_{\mathcal{H}},S_{\mathcal{T}}}(\hat{S}_{\mathcal{T}}).$

The RBF function is a space warp based on triharmonic radial basis functions
(Botsch and Kobbelt, 2005) that is calculated from the template skin surface
to the target and then applied to the template skin layer. By construction,
the skin layer will be warped semantically consistent and stick close to the
targeted skin surface.

Next, we fit the skull layer $\hat{B}_{\mathcal{H}}$ by invoking a linear
regressor $D$ that predicts the distances from the vertices of
$\hat{S}_{\mathcal{H}}$ to the corresponding vertices of
$\hat{B}_{\mathcal{H}}$ and subsequently minimizing with projective dynamics

(24) $\displaystyle\arg\min_{X}w_{r}{E}_{r}(X,\hat{S}_{\mathcal{T}})$
$\displaystyle+w_{d_{2}}E_{d_{2}}(\hat{S}_{\mathcal{H}},X,D(\hat{S}_{\mathcal{H}}))$
$\displaystyle+w_{c}E_{c}(X,\hat{B}_{\mathcal{T}}),$

where $E_{d_{2}}$ is implemented as edge strains to ensure that the predicted
distances are adhered to. Apart from $E_{d_{2}}$, the same regularizing terms
as in Equation LABEL:eq::rect and 17 are used. The optimization is initialized
with $X=\hat{S}_{\mathcal{H}}-D(S_{\mathcal{H}})\cdot
n(\hat{S}_{\mathcal{H}})$ where $n(\hat{S}_{\mathcal{H}})$ are area-weighted
vertex normals. $D$ is trained on the dataset of (Gietzen et al., 2019)
(SKULLS) that relates MRI skull measurements to skin surface scans. To ease
the learning task, we learn the regressor between PCAs of the skin layers and
the skin-to-skull-layer distances. Predicted lengths are set to a minimum
value if they fall below a threshold, thus, avoiding skin-skull intersections
and numerical issues in downstream physics-based simulations. In Figure 3 a),
the linear regressor training is visualized.

The muscle layer $\hat{\mathcal{M}}_{\mathcal{H}}$ is fitted by positioning
its vertices at the same absolute distances between the corresponding skin and
skull layer vertices as in the template, and only passing on a small relative
amount $w_{rel}$ of the distance changes compared to the template. This
approach assumes that the muscle mass in the facial area is only moderately
affected by body weight and skull size.

The skull mesh is placed by setting

(25)
$B_{\mathcal{H}}=\text{rbf}_{\hat{B}_{\mathcal{H}},\hat{B}_{\mathcal{T}}}(B_{\mathcal{T}}).$

The properties of the RBF space warp ensure that the skull mesh remains within
the skull layer if the layer is of sufficient resolution. The muscle mesh
could be placed in a similar fashion but is not needed in our pipeline any
further.

Finally, the soft and muscle tissue tet meshes $\mathbb{S}_{\mathcal{H}}$ and
$\mathbb{M}_{\mathcal{H}}$ can be constructed as described in Section 2. On
average, the complete fitting pipeline takes only about 500ms on an Intel i5
12600K processor. Figure 3 b) visualizes the overall fitting process.

Sampling Expressions In the second part of the pipeline, the actual training
instance consisting of $({S}_{\mathcal{H}}^{exp},\bar{S}_{\mathcal{H}}^{exp})$
is created. Beginning with deformation transfer (Botsch et al., 2006) to
transfer ARKit111https://developer.apple.com/ surface-based blendshapes to
$S_{\mathcal{H}}$. We then create the expression $S_{\mathcal{H}}^{exp}$ by
linearly blending the blendshapes with weights that are obtained from 8 around
10 minutes long dyadic conversations recorded with a custom iOS app.

Training Instance Finally,
$\bar{S}_{\mathcal{H}}^{exp}=\phi(\phi^{\dagger}(S_{\mathcal{H}}^{exp},\mathbb{S}_{\mathcal{H}},\mathbb{M}_{\mathcal{H}},B_{\mathcal{H}}))$
can be computed. Computing one training instance takes approximately 40
seconds on a AMD threadripper pro 3995wx.

###  3.5. Neural Architecture

After creating training data, we continue by building a suitable neural
network architecture for $f$ that is able to closely approximate
$\phi(\phi^{\dagger}(\cdot))$. To maintain the efficiency of linear
blendshapes and to make $f$ fast, $f$ must be of minimal complexity, both in
structure and in the number of parameters.

We start the construction of $f$ by designing input and output as modest as
feasible. Regarding the input, it is apparent that $f$ must be fed the
targeted linearly generated expression surface $S_{\mathcal{H}}^{exp}$.
However, $\phi^{\dagger}$ uses the partially tetrahedralized representation of
the LHM, too, which we want to refrain from to avoid simulation dependencies
for fitting the LHM. Therefore, we input the neutral head surface
$S_{\mathcal{H}}$ and expect $f$ to learn the fitting of the LHM implicitly.
As a side effect, the processing of the high resolution soft and muscle tissue
tet meshes $\mathbb{S}_{\mathcal{H}},\mathbb{M}_{\mathcal{H}}$ can be avoided,
as well. Both inputted surfaces, the neutral and the targeted one, are
tokenized with the help of respective encoders. The separate processing of
$S_{\mathcal{H}}$ only has to be performed once and does not have to be
repeated during an animation. To further simplify the encoding of
$S_{\mathcal{H}}^{exp}$, we input the vertex-wise deltas between the
expression and the neutral surface instead of the expression surface itself.
These deltas are mostly Gaussian distributed and can therefore be learned more
easily by the neural network (Li et al., 2020). In the same spirit, the tokens
are processed by a decoder that outputs only the vertex-wise deformations from
the targeted expression $S_{\mathcal{H}}^{exp}$ to the anatomically and
physically plausible expression $\bar{S}_{\mathcal{H}}^{exp}$.

Both the encoder of the input and the decoder of the output are implemented as
a fully connected networks. The complete architecture is depicted in Figure 4.
Please note that we evaluated many other network architectures such as set
transformers (Lee et al., 2019), convolutional networks on geometry images,
graph neural networks (Scarselli et al., 2008), or implicit architectures
(Mildenhall et al., 2021), but all have exhibited substantially slower
inference speeds while reaching a similar accuracy.

Our design of $f$ is not only fast, i.e. runs only 8 ms per frame on a
consumer grade Intel i5 12600K but it is also straightforward to deploy.
First, by its design as an adapter, it is independent of the linear blendshape
rig that is used to generate the targeted expressions. Second, a single
neutral surface of a head is sufficient on which deformation transfer can be
applied. Building on this, an anatomically plausible animation can be
performed with $f$ as exemplary shown in Figure 4. Finally, no more volumetric
information must be processed in complex simulation frameworks. Due to this
and because only simple fully connected layers are used, the mesh is portable
and easy to use on many different (computationally weak) devices.

##  4\. Experiments

Before demonstrating the accuracy and efficiency (Section 4.3) of neural
volumetric blendshapes, we first evaluate the fitting precision of the LHM
(Section 4.1) as well as the quality of the proposed physical models (Section
4.2).

###  4.1. LHM Fitting

Figure 5. The per-vertex mean L2-error of our LHM fitting algorithm in a
leave-one-out validation on SKULLS (Gietzen et al., 2019). Larger errors only
appear in the back part of the face and have only a minor impact on facial
animations.

The fitting of the LHM is mainly composed of the data-driven positioning of
the skull layer and the subsequent heuristic fitting of the muscle layer. We
evaluate the crucial fitting of the skull layer with the SKULLS (Gietzen et
al., 2019) dataset. Since this dataset consists of 43 instances only, a leave-
one-out validation is performed in which the vertex-wise L2 errors are
measured. The results are compared to the multilinear model as originally used
for SKULLS in (Achenbach et al., 2018).

Figure 6. Exemplary fits of the LHM components skull wrap, muscle wrap, and
skull.

Both models cannot achieve a medical-grade positioning with errors between
approximately 2 and 4 mm. The MLM achieves a higher precision with a mean
error of 1.98mm than our approach that dispositions the skull by 3.83mm on
average. However, the MLM cannot prevent collisions that might crash physics-
based simulations. Also, our fitting algorithm produces large errors only in
regions that are of less importance for facial simulations as can be seen in
Figure 5. The errors are predominately distributed in the back area of the
skull, since here the rectangular constraints of our fitting procedure can
presumably no longer be aligned well to the skin layer. The following section
demonstrates in downstream physics-based simulations that the prediction
quality in the frontal face region is sufficiently adequate for detailed
facial animations. Figure 6 displays fitting examples.

###  4.2. Physics-Based Simulation

Figure 7. Exemplary results of the proposed physical models. In the top row, muscle deformation gradients have been extracted with $\phi^{\dagger}$ and reapplied using $\phi$ to the corresponding neutral head. In the bottom row, the deformation gradients are reapplied to a different identity. $w_{c}$ | $w_{str}^{skull}$ | $w_{B}$ | $w_{str}^{tet}$ | $w_{vol}$ | $w_{con}$ | $w_{inv}$ |
---|---|---|---|---|---|---|---
$10^{3}$ | $10^{3}$ | $10^{3}$ | $10^{2}$ | $10^{3}$ | $10^{3}$ | $10$ |
$w_{{\mathbb{S}_{\mathcal{H}}}}$ | $w_{{\mathbb{M}_{\mathcal{H}}}}$ | $w_{{B_{\mathcal{H}}}}$ | $w_{c}$ | $w_{r}$ | $w_{d}$ | $w_{d_{2}}$ |
$1$ | $1$ | $10$ | $0.1$ | $1.0$ | $10.0$ | $10.0$ |
Table 2. The weights used to implement $\phi,\phi^{\dagger}$, and the LHM
fitting. Figure 8. Exemplary test results of our neural volumetric
blendshapes. Besides resolving self-collisions as, the physical and anatomical
constraints anchored in our approach lead to more realistic and immersive
expressions while being computationally efficient.

To investigate the quality of $\phi^{\dagger}$ and $\phi$ we perform an
experiment on nearly 300 high-resolution optical facial expression scans from
the proprietary 3DScanstore222https://www.3dscanstore.com dataset. First, we
fit the LHM to each of the scans using $\phi^{\dagger}$ and store only the
resulting muscle deformation gradients and bone movements. Then, in turn, we
apply these to the neutral scans of the corresponding individuals through
$\phi$ and simulate the soft tissue as $\phi^{\dagger}$ does. Only if our
physic models and the previous LHM fitting operate with reasonable precision,
do we expect a low L2-loss between the original expression and the simulated
expression as well as visually appealing results.

Across all expression scans, we observe an average L2 reconstruction error of
only 1.6 mm. The reproduction quality of the physical models is also reflected
in the visual results of Figure 7. Further, the accuracy of the LHM fitting
investigated in the previous section is again underlined. Additionally, we
inspected expression retargeting by applying the extracted muscle and bone
transformations to other identities. Again, visual appealing examples can be
found Figure 7. The weights used to realize $\phi,\phi^{\dagger}$, and the LHM
fitting in this and the following experiments are stated in Table 2.

###  4.3. Neural Volumetric Blendshapes

Figure 9. Exemplary results of our collision resolution. Collisions are not
only avoided but also simulated in a physically plausible manner.

To train and evaluate $f$, we assemble a dataset of 50000 training and test
instances by using the pipeline from Section 3.4. At this, we sample 10000
different head shapes and compute 5 different expressions per head. For
training, the Adam optimizer performs 200k update steps with a learning rate
of 0.0001. The learning rate is linearly decreased to 0.00005 over the course
of training and a batch size of 128 is used. In total, the training
specifications result in an approximate runtime of 8 hours on an NVIDIA A6000.
The comparatively short training time can straightforwardly be explained by
the less noisy training data as usually encountered for image-based deep
learning models and the efficient network design. The identities are used 90%
for training and 10% for testing.

First, we evaluate the efficiency of $f$. The initial runtime of
$\phi(\phi^{\dagger}(\cdot))$ with 6 global projective dynamics updates as in
(Ichim et al., 2016) is about 950ms on a workstation AMD threadripper pro
3995wx with 128 cores (implemented with ShapeOp (Deuss et al., 2015)). Our
implementation of $f$ only takes 8ms on a consumer-grade Intel i5 12600K
(implemented with PyTorch333https://pytorch.org). Moreover, a forward pass on
an NVIDIA RTX3090 is calculated in less than one millisecond. Thus, the neural
volumetric blendshapes are suitable for realtime applications on weaker
hardware setups and also offer advantages when many facial animations are to
be executed in parallel on a GPU.

The time advantage alone is not useful if $f$ is not an adequate
approximation. We evaluate the approximation quality by measuring the mean L2
error on the test dataset. More precisely, we randomly generate 5 training-
test splits of the generated dataset and report the average of the test losses
after training on each training set. We achieve an average test error of only
0.13mm, meaning that we can replicate $\phi(\phi^{\dagger}(\cdot))$ very
accurately and $f$ generalizes well across head shapes and facial expressions.
Another open question is the temporal consistency of $f$ which is inherent in
physics-based simulations. To this end, we invite the reader to watch the
attached video in the supplementary material. Finally, visual results of
neural volumetric blendshapes as shown in Figure 9 exhibit the same
improvements towards more realistic facial expressions as much slower physics-
based simulations before. These include the resolution of self-intersections,
the simulation of wrinkles, the adherence to anatomical boundaries, and the
realism-enhancing constraining of volume preservation.

##  5\. Conclusion

In this work, we presented neural volumetric blendshapes, an efficient
realization of physics-based facial animations even on consumer-grade
hardware. Our approach yields more realistic results than the commonly used
linear blendshapes, since artifacts such as volume-loss, self-intersections,
or anatomical flaws are avoided and effects like volumetric elasticity and
wrinkles are added. It is also convenient to use, since existing surface-based
blendshape rigs can be extended into anatomically and physically plausible
animations with no integration overhead while keeping the actuation as before.

We aim to improve neural volumetric blendshapes in at least two directions. On
the one hand, with an even more accurate anatomical model that represents e.g.
trachea and esophagus more precisely. On the other hand, recent results
(Romero et al., 2022) show that contact deformations can also be learned.
Since people touch their faces dozens of times (Spille et al., 2021) a day,
adding contact-handling for more realistic gestures may improve the immersion
significantly.

## References

  * (1)
  * Achenbach et al. (2018) Jascha Achenbach, Robert Brylka, Thomas Gietzen, Katja zum Hebel, Elmar Schömer, Ralf Schulze, Mario Botsch, and Ulrich Schwanecke. 2018\.  A multilinear model for bidirectional craniofacial reconstruction. In _Proceedings of the Eurographics Workshop on Visual Computing for Biology and Medicine_. 67–76.
  * Ali-Hamadi et al. (2013) Dicko Ali-Hamadi, Tiantian Liu, Benjamin Gilles, Ladislav Kavan, François Faure, Olivier Palombi, and Marie-Paule Cani. 2013.  Anatomy transfer.  _ACM transactions on graphics (TOG)_ 32, 6 (2013), 1–8.
  * Arsigny et al. (2007) Vincent Arsigny, Pierre Fillard, Xavier Pennec, and Nicholas Ayache. 2007.  Geometric means in a novel vector space structure on symmetric positive-definite matrices.  _SIAM journal on matrix analysis and applications_ 29, 1 (2007), 328–347.
  * Athar et al. (2022) ShahRukh Athar, Zexiang Xu, Kalyan Sunkavalli, Eli Shechtman, and Zhixin Shu. 2022.  RigNeRF: Fully Controllable Neural 3D Portraits. In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_. 20364–20373.
  * Bao et al. (2019) Michael Bao, Matthew Cong, Stéphane Grabli, and Ronald Fedkiw. 2019.  High-quality face capture using anatomical muscles. In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_. 10802–10811.
  * Barrielle et al. (2016) Vincent Barrielle, Nicolas Stoiber, and Cédric Cagniart. 2016\.  Blendforces: A dynamic framework for facial animation. In _Computer Graphics Forum_ , Vol. 35. Wiley Online Library, 341–352.
  * Bickel et al. (2008) Bernd Bickel, Manuel Lang, Mario Botsch, Miguel A Otaduy, and Markus H Gross. 2008.  Pose-Space Animation and Transfer of Facial Details.. In _Symposium on Computer Animation_. 57–66.
  * Botsch and Kobbelt (2005) Mario Botsch and Leif Kobbelt. 2005.  Real-time shape editing using radial basis functions. In _Computer graphics forum_ , Vol. 24. Blackwell Publishing, Inc Oxford, UK and Boston, USA, 611–621.
  * Botsch et al. (2010) Mario Botsch, Leif Kobbelt, Mark Pauly, Pierre Alliez, and Bruno Lévy. 2010.  _Polygon mesh processing_.  CRC press.
  * Botsch et al. (2006) Mario Botsch, Robert Sumner, Mark Pauly, and Markus Gross. 2006\.  Deformation transfer for detail-preserving surface editing. In _Vision, Modeling & Visualization_. Citeseer, 357–364.
  * Bouaziz et al. (2014) Sofien Bouaziz, Sebastian Martin, Tiantian Liu, Ladislav Kavan, and Mark Pauly. 2014.  Projective dynamics: Fusing constraint projections for fast simulation.  _ACM transactions on graphics (TOG)_ 33, 4 (2014), 1–11.
  * Bradley et al. (2010) Derek Bradley, Wolfgang Heidrich, Tiberiu Popa, and Alla Sheffer. 2010.  High resolution passive facial performance capture.  In _ACM SIGGRAPH 2010 papers_. 1–10.
  * Cao et al. (2022) Chen Cao, Tomas Simon, Jin Kyu Kim, Gabe Schwartz, Michael Zollhoefer, Shun-Suke Saito, Stephen Lombardi, Shih-En Wei, Danielle Belko, Shoou-I Yu, et al. 2022\.  Authentic volumetric avatars from a phone scan.  _ACM Transactions on Graphics (TOG)_ 41, 4 (2022), 1–19.
  * Casas and Otaduy (2018) Dan Casas and Miguel A Otaduy. 2018.  Learning nonlinear soft-tissue dynamics for interactive avatars.  _Proceedings of the ACM on Computer Graphics and Interactive Techniques_ 1, 1 (2018), 1–15.
  * Choi et al. (2022) Byungkuk Choi, Haekwang Eom, Benjamin Mouscadet, Stephen Cullingford, Kurt Ma, Stefanie Gassel, Suzi Kim, Andrew Moffat, Millicent Maier, Marco Revelant, et al. 2022\.  Animatomy: an Animator-centric, Anatomically Inspired System for 3D Facial Modeling, Animation and Transfer. In _SIGGRAPH Asia 2022 Conference Papers_. 1–9.
  * Cong and Fedkiw (2019) Matthew Cong and Ronald Fedkiw. 2019.  Muscle-based facial retargeting with anatomical constraints.  In _ACM SIGGRAPH 2019 Talks_. 1–2.
  * Cong (2016) Matthew Deying Cong. 2016\.  _Art-directed muscle simulation for high-end facial animation_.  Stanford University.
  * Deuss et al. (2015) Mario Deuss, Anders Holden Deleuran, Sofien Bouaziz, Bailin Deng, Daniel Piker, and Mark Pauly. 2015\.  ShapeOp—a robust and extensible geometric modelling paradigm.  In _Modelling Behaviour_. Springer, 505–515.
  * Feng et al. (2021) Yao Feng, Haiwen Feng, Michael J Black, and Timo Bolkart. 2021\.  Learning an animatable detailed 3D face model from in-the-wild images.  _ACM Transactions on Graphics (ToG)_ 40, 4 (2021), 1–13.
  * Garbin et al. (2022) Stephan J Garbin, Marek Kowalski, Virginia Estellers, Stanislaw Szymanowicz, Shideh Rezaeifar, Jingjing Shen, Matthew Johnson, and Julien Valentin. 2022.  VolTeMorph: Realtime, Controllable and Generalisable Animation of Volumetric Representations.  _arXiv preprint arXiv:2208.00949_ (2022).
  * Gietzen et al. (2019) Thomas Gietzen, Robert Brylka, Jascha Achenbach, Katja Zum Hebel, Elmar Schömer, Mario Botsch, Ulrich Schwanecke, and Ralf Schulze. 2019.  A method for automatic forensic facial reconstruction based on dense statistics of soft tissue thickness.  _PloS one_ 14, 1 (2019), e0210257.
  * Gilles et al. (2010) Benjamin Gilles, Lionel Reveret, and Dinesh K Pai. 2010\.  Creating and animating subject-specific anatomical models. In _Computer Graphics Forum_ , Vol. 29. Wiley Online Library, 2340–2351.
  * Ichim et al. (2015) Alexandru Eugen Ichim, Sofien Bouaziz, and Mark Pauly. 2015\.  Dynamic 3D avatar creation from hand-held video input.  _ACM Transactions on Graphics (ToG)_ 34, 4 (2015), 1–14.
  * Ichim et al. (2017) Alexandru-Eugen Ichim, Petr Kadleček, Ladislav Kavan, and Mark Pauly. 2017.  Phace: Physics-based face modeling and animation.  _ACM Transactions on Graphics (TOG)_ 36, 4 (2017), 1–14.
  * Ichim et al. (2016) Alexandru Eugen Ichim, Ladislav Kavan, Merlin Nimier-David, and Mark Pauly. 2016.  Building and animating user-specific volumetric face rigs.. In _Symposium on Computer Animation_. 107–117.
  * Jung et al. (2015) Sungkyu Jung, Armin Schwartzman, and David Groisser. 2015\.  Scaling-rotation distance and interpolation of symmetric positive-definite matrices.  _SIAM J. Matrix Anal. Appl._ 36, 3 (2015), 1180–1201.
  * Kadleček et al. (2016) Petr Kadleček, Alexandru-Eugen Ichim, Tiantian Liu, Jaroslav Křivánek, and Ladislav Kavan. 2016.  Reconstructing personalized anatomical models for physics-based body animation.  _ACM Transactions on Graphics (TOG)_ 35, 6 (2016), 1–13.
  * Kadleček and Kavan (2019) Petr Kadleček and Ladislav Kavan. 2019.  Building accurate physics-based face models from data.  _Proceedings of the ACM on Computer Graphics and Interactive Techniques_ 2, 2 (2019), 1–16.
  * Karras et al. (2019) Tero Karras, Samuli Laine, and Timo Aila. 2019.  A style-based generator architecture for generative adversarial networks. In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_. 4401–4410.
  * Keller et al. (2022) Marilyn Keller, Silvia Zuffi, Michael J Black, and Sergi Pujades. 2022.  OSSO: Obtaining Skeletal Shape from Outside. In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_. 20492–20501.
  * Komaritzan and Botsch (2018) Martin Komaritzan and Mario Botsch. 2018.  Projective skinning.  _Proceedings of the ACM on Computer Graphics and Interactive Techniques_ 1, 1 (2018), 1–19.
  * Komaritzan et al. (2021) Martin Komaritzan, Stephan Wenninger, and Mario Botsch. 2021\.  Inside Humans: Creating a Simple Layered Anatomical Model from Human Surface Scans.  _Frontiers in Virtual Reality_ 2 (2021), 694244.
  * Lee et al. (2019) Juho Lee, Yoonho Lee, Jungtaek Kim, Adam Kosiorek, Seungjin Choi, and Yee Whye Teh. 2019\.  Set transformer: A framework for attention-based permutation-invariant neural networks. In _International conference on machine learning_. PMLR, 3744–3753.
  * Lewis et al. (2014) John P Lewis, Ken Anjyo, Taehyun Rhee, Mengjie Zhang, Frederic H Pighin, and Zhigang Deng. 2014\.  Practice and theory of blendshape facial models.  _Eurographics (State of the Art Reports)_ 1, 8 (2014), 2\.
  * Lewis et al. (2005) John P Lewis, Jonathan Mooser, Zhigang Deng, and Ulrich Neumann. 2005.  Reducing blendshape interference by selected motion attenuation. In _Proceedings of the 2005 symposium on Interactive 3D graphics and games_. 25–29.
  * Li et al. (2020) Jiaman Li, Zhengfei Kuang, Yajie Zhao, Mingming He, Karl Bladin, and Hao Li. 2020\.  Dynamic facial asset and rig generation from a single scan.  _ACM Trans. Graph._ 39, 6 (2020), 215–1.
  * Loper et al. (2015) Matthew Loper, Naureen Mahmood, Javier Romero, Gerard Pons-Moll, and Michael J. Black. 2015\.  SMPL: A Skinned Multi-Person Linear Model.  _ACM Trans. Graphics (Proc. SIGGRAPH Asia)_ 34, 6 (Oct. 2015), 248:1–248:16.
  * Maalin et al. (2021) Nadia Maalin, Sophie Mohamed, Robin SS Kramer, Piers L Cornelissen, Daniel Martin, and Martin J Tovée. 2021.  Beyond BMI for self-estimates of body size and shape: A new method for developing stimuli correctly calibrated for body composition.  _Behavior Research Methods_ 53, 3 (2021), 1308–1321.
  * Mildenhall et al. (2021) Ben Mildenhall, Pratul P Srinivasan, Matthew Tancik, Jonathan T Barron, Ravi Ramamoorthi, and Ren Ng. 2021\.  Nerf: Representing scenes as neural radiance fields for view synthesis.  _Commun. ACM_ 65, 1 (2021), 99–106.
  * Parke (1991) Frederic I Parke. 1991\.  Control parameterization for facial animation. In _Computer Animation’91_. Springer, 3–14.
  * Romero et al. (2022) Cristian Romero, Dan Casas, Maurizio M Chiaramonte, and Miguel A Otaduy. 2022.  Contact-centric deformation learning.  _ACM Transactions on Graphics (TOG)_ 41, 4 (2022), 1–11.
  * Saito et al. (2015) Shunsuke Saito, Zi-Ye Zhou, and Ladislav Kavan. 2015\.  Computational bodybuilding: Anatomically-based modeling of human bodies.  _ACM Transactions on Graphics (TOG)_ 34, 4 (2015), 1–12.
  * Santesteban et al. (2020) Igor Santesteban, Elena Garces, Miguel A Otaduy, and Dan Casas. 2020.  SoftSMPL: Data-driven Modeling of Nonlinear Soft-tissue Dynamics for Parametric Humans. In _Computer Graphics Forum_ , Vol. 39. Wiley Online Library, 65–75.
  * Scarselli et al. (2008) Franco Scarselli, Marco Gori, Ah Chung Tsoi, Markus Hagenbuchner, and Gabriele Monfardini. 2008\.  The graph neural network model.  _IEEE transactions on neural networks_ 20, 1 (2008), 61–80.
  * Schleicher et al. (2021) Robert Schleicher, Marlies Nitschke, Jana Martschinke, Marc Stamminger, Björn M Eskofier, Jochen Klucken, and Anne D Koelewijn. 2021\.  BASH: Biomechanical Animated Skinned Human for Visualization of Kinematics and Muscle Activity.. In _VISIGRAPP (1: GRAPP)_. 25–36.
  * Sifakis et al. (2005) Eftychios Sifakis, Igor Neverov, and Ronald Fedkiw. 2005\.  Automatic determination of facial muscle activations from sparse motion capture marker data.  In _ACM SIGGRAPH 2005 Papers_. 417–425.
  * Song et al. (2020) Steven L Song, Weiqi Shi, and Michael Reed. 2020.  Accurate face rig approximation with deep differential subspace reconstruction.  _ACM Transactions on Graphics (TOG)_ 39, 4 (2020), 34–1.
  * Spille et al. (2021) Jente L Spille, Martin Grunwald, Sven Martin, and Stephanie M Mueller. 2021.  Stop touching your face! A systematic review of triggers, characteristics, regulatory functions and neuro-physiology of facial self touch.  _Neuroscience & Biobehavioral Reviews_ 128 (2021), 102–116.
  * Srinivasan et al. (2021) Sangeetha Grama Srinivasan, Qisi Wang, Junior Rojas, Gergely Klár, Ladislav Kavan, and Eftychios Sifakis. 2021.  Learning active quasistatic physics-based models from data.  _ACM Transactions on Graphics (TOG)_ 40, 4 (2021), 1–14.
  * Yang et al. (2022) Lingchen Yang, Byungsoo Kim, Gaspard Zoss, Baran Gözcü, Markus Gross, and Barbara Solenthaler. 2022\.  Implicit neural representation for physics-driven actuated soft bodies.  _ACM Transactions on Graphics (TOG)_ 41, 4 (2022), 1–10.
  * Zhang et al. (2008) Li Zhang, Noah Snavely, Brian Curless, and Steven M Seitz. 2008\.  Spacetime faces: High-resolution capture for~ modeling and animation.  In _Data-Driven 3D Facial Animation_. Springer, 248–276.
  * Zheng et al. (2022) Yufeng Zheng, Victoria Fernández Abrevaya, Marcel C Bühler, Xu Chen, Michael J Black, and Otmar Hilliges. 2022.  Im avatar: Implicit morphable head avatars from videos. In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_. 13545–13555.


∎

e1e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
11institutetext: Departamento de Física, Universidad de Antofagasta, Aptdo
02800, Chile.

# Wave zone in the Hořava-Lifshitz theory at the kinetic-conformal point in
the low energy regime

J. Mestra-Páeze1,addr1  J.M. Peñae2,addr1  A. Restucciae3,addr1

(Received: date / Accepted: date)

###### Abstract

We show that in the Hořava-Lifshitz theory at the kinetic-conformal point, in
the low energy regime, a wave zone for asymptotically flat fields can be
consistently defined. In it, the physical degrees of freedom, the transverse
traceless tensorial modes, satisfy a linear wave equation. The Newtonian
contributions, among which there are terms which manifestly break the
relativistic invariance, are non-trivial but do not obstruct the free
propagation (radiation) of the physical degrees of freedom. For an appropriate
value of the couplings of the theory, the wave equation becomes the
relativistic one in agreement with the propagation of the gravitational
radiation in the wave zone of General Relativity. Previously to the wave zone
analysis, and in general grounds, we obtain the physical Hamiltonian of the
Hořava-Lifshitz theory at the kinetic-conformal point in the constrained
submanifold. We determine the canonical physical degrees of freedom in a
particular coordinate system. They are well defined fuctions of the
transverse-traceless modes of the metric and coincide with them in the wave
zone and also at linearized level.

###### Keywords:

Gravitational-Waves Hořava-Lifshitz Wave-Zone.

††journal: Eur. Phys. J. C

##  1 Introduction

The detection of gravitational waves has opened a new era in the study of
physics Miller2019 ; Nitz2021 . Multi-messenger astronomy will be decisive in
the study of astrophysical and cosmological events and can lead to the
discovery of new phenomena in extreme situations beyond the reach of
experimental tests that we now carry out.

Hořava-Lifshitz gravity is a recent proposal for a candidate to ultraviolet
completion of General Relativity (GR) horava2009 ; Blas2010 . This theory
models the gravity as 4-dimensional differentiable manifold with a foliation-
structure of co-dimension one. The foliation-leaves are 3-dimensional
Riemannian submanifolds. In addition, the time and space scale in different
ways, consequently the relativistic symmetry is manifestly broken. The
anisotropic scaling allows to include interaction terms with high spatial
derivatives in the potential, without breaking the symmetry of the action
under diffeomorphisms that preserve the foliation while keeping the second-
order time derivatives of the kinetic term. The theory contains several
coupling constants. There is only one in the kinetic term of the Hořava-
Lifshitz action, it is dimensionless and plays a relevant role in the theory.
When its value is $\lambda=1/3$, the so-called kinetic-conformal point, the
theory propagates, at the linearized level, the same degrees of freedom of
linearized General Relativity and with an appropriate choice of coupling
parameters, it is consistent with low energy experiments EmirGumrukcuoglu2018
; BellorinRestucciaSotomayor_2013 . At linearized level the non-propagating
components of the metric and the lapse function become zero in both theories
as a consequence of the constraints.

The original Hořava-Lifshitz gravity horava2009  suffers from a strong
coupling problem associated to the scalar mode of the theory. An improved
formulation was obtained by the inclusion of a new interaction term proposed
by Blas, Pujolàs and Sibiryakov (BPS), quadratic in derivatives and compatible
with the foliation preserving diffeomorphisms Blas2010 ; Blass2010a  and the
corresponding contributions to the potential. This extended formulation is
free from the strong coupling problem of the original Hořava-Lifshitz version.
It is a weakly coupled anisotropic description of gravity.

The Hořava-Lifshitz gravity for $\lambda=1/3$, with the inclusion of the BPS
interaction term, does not propagate, at any energy scale, the scalar mode
which occurs for $\lambda\neq 1/3$. Therefore, it does not have any strong
coupling problem. The inclusion of the BPS interaction term is essential to
have a formulation of the second class constraints in terms of strictly
elliptic partial differential equations, which allow a consistent elimination
of non-physical modes. Besides, the field equations for the $\lambda=1/3$
model, evaluated at $\alpha=0$ and $\beta=1$ coincide exactly with the General
Relativity field equations in the gauge $\pi=0$ (an admissible gauge condition
outside the black hole horizon in GR). This gauge was used in the ADM analysis
of the dynamics of GR ADM2008 . Moreover, for the models with $\lambda\neq
1/3$ the only dependence on $\lambda$ in the Hamiltonian is through the term
proportional to $\pi^{2}$ given by

$\frac{N}{\sqrt{g}}\frac{\lambda}{(3\lambda–1)}\pi^{2}.$ (1)

On the other hand, in the Hamiltonian of the $\lambda=1/3$ model there is no
$\pi^{2}$ term, since there is a second class constraint $\pi=0$ which arises
directly as a primary constraint from the canonical formulation of the theory.
The term $\pi^{2}$ could in principle be generated from the renormalization
group flow, since it is an admissible term in the Hamiltonian under the
foliation preserving diffeomorphisms. However, the second class constraint
which has to be imposed at all times prevents its appearance. Consequently,
the BPS extension of the $\lambda=1/3$ Hořava-Lifshitz gravity is a viable
model, describing anisotropic gravity.

There are then three non projectable viable models of anisotropic gravity: the
healthy extension of Hořava-Lifshitz Blas2010  for $\lambda\neq 1/3$, which
describes the propagation of transverse-traceless tensorial modes together
with a scalar one, the $U(1)$ symmetric models Horava2010 ; WangAnzohongWu2011
; daSilva2011 ; Kluson2011 ; Linetal2011  for the projectable models and
ZhuEtal2011  for the non projectable models, and finally the $\lambda=1/3$
model of Hořava-Lifshitz gravity, the conformal kinetic model, which with the
inclusion of the BPS interaction term, propagates only the transverse-
traceless tensorial modes. We will consider in this paper this latest model.
The theory ends up being power counting renormalizable and unitary
Bellorin_Restuccia_2016B ; charmousis2009 ; visser2009 ; papazoglou2010 ;
orlando2009 ; shu2009 ; benedetti2014 ; contillo2013 ; d2014asymptotic ;
d2015covariant ; barvinsky2016 ; wang2017 ; shin2017 ; pospelov2012 .

Beyond the linearized formulation, GR has a well defined wave zone, on which
the physical degrees of freedom propagate freely on a nontrivial Newtonian
background. This is a nontrivial property of GR not necessarily valid for
nonlinear theories. Hence, we may wonder if such a property is also valid for
Hořava-Lifshitz gravity, and if it is the case what are the effects of the
anisotropy, that is, the breaking of the Lorentz symmetry, on the wave zone.

Arnowitt, Deser and Misner (ADM) proved in the 60’s that General Relativity
has a well defined wave zone Arnowitt1961 . In this space-time region the
metric components $g^{TT}_{ij}$ at order $\mathcal{O}(1/r)$ satisfy the same
wave equation as in the linearized theory around a Minkowski space-time. In
the wave zone there also exists a Newtonian background at order
$\mathcal{O}(1/r)$ that does not prevent the $g^{TT}_{ij}$ modes, the physical
degrees of freedom, to propagate as free radiation Arnowitt1961 ; ADM2008 .

We prove in this work that in the $(\lambda=1/3)$-Hořava-Lifshitz gravity at
low energies, although the relativistic symmetry is broken, there is a wave
zone where the physical degrees of freedom, the transverse traceless tensorial
modes, propagate freely satisfying a wave equation, without any interaction
with the nontrivial Newtonian background as in GR. This means that the low
energy results obtained for the wave propagation in GR, in the wave zone, are
also valid for the Hořava-Lifshitz gravity theory.

We emphasize that although in the wave zone the dynamical equation satisfied
by the physical degrees of freedom is the same as the one obtained from a
linearized analysis, in the wave zone there are nontrivial Newtonian
components of the metric (background). This Newtonian background, a nonlinear
effect of the analysis, determine the energy and momentum of the gravitational
wave. This situation also happens in GR, but with different expressions for
the energy and momentum Blas2011 .

We show, using the local symmetries, that the Hořava-Lifshitz theory, at the
kinetic-conformal point, can be canonically reduced to a Hamiltonian
formulation in terms of the physical degrees of freedom, in a particular
coordinate frame. The Hamiltonian density explicitly depends on the Lorentz
violating term characterized by the coupling parameter $\alpha$.

In section 2, we review the Hořava-Lifshitz theory at the kinetic-conformal
point. In section 3, we show that there is a Hamiltonian formulation in terms
of the physical degrees of freedom. We obtain the Hamiltonian density in a
particular coordinate system. In section 4, we obtain the formulation of the
theory in the wave zone. Finally, in section 5, we give the conclusions of our
paper.

##  2 Foliation, geometry and Hořava-Lifshitz gravity

Let $M$ be a 4-dimensional differentiable manifold. $M$ has a codimension-one
foliation structure $(M,\mathcal{F})$ if the maximal atlas
$\mathcal{F}\equiv(U_{\mathcal{A}},\varphi_{\mathcal{A}})$ i.e
$M=\cup\,U_{\mathcal{A}}$, where $U_{\mathcal{A}}$ is a family of open subsets
of $M$ and $\varphi_{\mathcal{A}}:U_{\mathcal{A}}\rightarrow
D_{\mathcal{A}}\subset\mathbb{R}^{1}\times\mathbb{R}^{3}$ are diffeomorphism
such that if $U_{i}\cap U_{j}\neq\emptyset$ the transition of charts is
defined by

$\displaystyle\varphi_{i}\circ\varphi^{-1}_{j}:\varphi_{j}(U_{i}\cap
U_{j})\rightarrow\varphi_{i}(U_{i}\cap U_{j}),$ (2)
$\displaystyle(t,x)\rightarrow(\tilde{t}(t),\tilde{x}(t,x)).$ (3)

The couple $(M,\mathcal{F})$ and its equivalents under the diffeomorphisms,
that preserve the foliation structure, $\mathcal{F}_{Diff}$, provide the
geometrical structure of the Hořava-Lifshitz theory where space and time scale
anisotropically $t\rightarrow b^{z}t$ and $x\rightarrow bx$. We remark that
$M$ is the disjoint union of 3-dimensional Riemannian manifolds
$(\Sigma_{t},g_{ij})$, $g_{ij}(t,x)=\frac{\partial\tilde{x}^{l}}{\partial
x^{i}}\frac{\partial\tilde{x}^{m}}{\partial
x^{j}}\tilde{g}_{lm}(\tilde{t},\tilde{x})$, where the following geometric
objects compatibles with the foliation structure are introduced: a proper time
defined through the introduction of the lapse $N$ and a shift of the spatial
coordinates defined through $N^{i}$ in order to have a contravariant
transformation law under $\mathcal{F}_{Diff}$.

$\displaystyle\tilde{N}(\tilde{t},\tilde{x})d\tilde{t}=N(t,x)dt,$ (4)
$\displaystyle
d\tilde{x}^{i}+\tilde{N}^{i}(\tilde{t},\tilde{x})d\tilde{t}=\frac{\partial\tilde{x}^{i}}{\partial
x^{j}}[dx^{j}+N^{j}(t,x)dt]\,,$ (5)

we emphasize that there is not a space-time metric on $M$. The metric on the
leaves of the foliation $g_{ij}$ and the fields $N$ and $N^{i}$ are used to
describe the evolution of the gravitational field. They scale anisotropically
as $g_{ij}\rightarrow b^{0}g_{ij}$, $N\rightarrow b^{0}N$ and
$N^{i}\rightarrow b^{z-1}N^{i}$.

Taking into account the anisotropic scaling and foliation structure the
proposal incorporates terms with high spatial derivatives in the potential
without breaking the symmetry under $\mathcal{F}_{Diff}$. The Hamiltonian of
the Hořava-Lifshitz gravity theory at the kinetic-conformal point is given by
RestucciaTello2020 ; BellorinRestucciaTello2018 ; Restuccia-Tello_2020 ;
BellorinRestucciaSotomayor_2013

$\displaystyle
H=\int_{\Sigma_{t}}d^{3}x\bigg{\\{}N\sqrt{g}\left[\frac{\pi^{ij}\pi_{ij}}{g}-\mathcal{V}\left(g_{ij},N\right)\right]$
$\displaystyle-N_{j}H^{j}-\sigma P_{N}-\mu\pi\bigg{\\}}+\beta E_{ADM}\,,$ (6)

here $\pi^{ij}$ is the canonical conjugate of $g_{ij}$, $N_{i}\equiv
g_{ij}N^{j}$, $\sigma$ and $\mu$ are Lagrange multipliers. The surface
integral

$\displaystyle E_{ADM}$ $\displaystyle\equiv$
$\displaystyle\oint_{\partial\Sigma_{t}}\left(\partial_{j}g_{ij}-\partial_{i}g_{jj}\right)dS_{i}\,,$
(7)

is added in order to ensure the Fréchet differentiability of the Hamiltonian,
see Regge-Teitelboim1974  where this idea was introduced and it was shown that
$E_{ADM}$ is the ADM energy in GR. The potential, up to terms quadratic on the
Riemann tensor and the vector field $a_{i}$, is
$\mathcal{V}=\mathcal{V}^{(1)}+\mathcal{V}^{(2)}+\mathcal{V}^{(3)}$, with
Bellorin_Restuccia_2016B

$\displaystyle\mathcal{V}^{(1)}$ $\displaystyle=$ $\displaystyle\beta R+\alpha
a_{i}a^{i},$ (8) $\displaystyle\mathcal{V}^{(2)}$ $\displaystyle=$
$\displaystyle\alpha_{1}R\nabla_{i}a^{i}+\alpha_{2}\nabla_{i}a_{j}\nabla^{i}a^{j}+\beta_{1}R_{ij}R^{ij}+\beta_{2}R^{2}$
(9) $\displaystyle\mathcal{V}^{(3)}$ $\displaystyle=$
$\displaystyle\alpha_{3}\nabla^{2}R\nabla_{i}a^{i}+\alpha_{4}\nabla^{2}a_{i}\nabla^{2}a^{i}+\beta_{3}\nabla_{i}R_{jk}\nabla^{i}R^{jk}$
(10) $\displaystyle+\beta_{4}\nabla_{i}R\nabla^{i}R\,,$

where $a_{i}\equiv\frac{1}{N}\partial_{i}N$, $\nabla_{i}$ represents the
affine connection, the covariant derivative constructed with the Riemannian
metric on the leaves, $\alpha$’s and $\beta$’s are coupling constants. The
potential also contains terms of the same order in the spatial derivatives as
the ones explicitly shown, but of cubic order or greater on the Riemann tensor
and the vector field $a_{i}$. Since they do not contribute to the dominant
order in the wave zone they have not been presented in (9) and (10). The
primary constraints of the theory are

$\displaystyle H^{j}=2\nabla_{i}\pi^{ij}=0,\quad
P_{N}=0,\quad\pi=g_{ij}\pi^{ij}=0\,.$ (11)

If we consider only the low energy potential, $\mathcal{V}^{(1)}$, the time
preservation of primary constrains imply the following secondary constraints:

$\displaystyle H_{P}$ $\displaystyle\equiv$
$\displaystyle\frac{3}{2}\frac{1}{\sqrt{g}}\pi^{ij}\pi_{ij}+\frac{1}{2}\sqrt{g}\beta
R$ (12)
$\displaystyle+\sqrt{g}\left(\frac{\alpha}{2}-2\beta\right)a_{i}a^{i}-2\beta\sqrt{g}\nabla^{i}a_{i}=0,$
$\displaystyle H_{N}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{\sqrt{g}}\left(\pi^{ij}\pi_{ij}-\beta
gR\right)+\alpha\sqrt{g}a_{i}a^{i}+$ (13) $\displaystyle
2\alpha\sqrt{g}\nabla_{i}a^{i}=0\,,$

which together with the last two constraints in (11) are second class
constraints, while $H^{j}=0$ is a first class constraint.

The evolution equations are,

$\displaystyle\partial_{t}{g}_{ij}$ $\displaystyle=$
$\displaystyle\frac{2N}{\sqrt{g}}\pi_{ij}+2\nabla_{(i}N_{j)}-\mu g_{ij}\,,$
(14) $\displaystyle\partial_{t}{\pi}^{ij}$ $\displaystyle=$
$\displaystyle\frac{N}{2}\frac{g^{ij}}{\sqrt{g}}\pi^{kl}\pi_{kl}-\frac{2N}{\sqrt{g}}\pi^{ik}\pi^{j}{}_{k}+N\sqrt{g}\beta\left(\frac{R}{2}g^{ij}-R^{ij}\right)$
(15) $\displaystyle-\alpha
N\sqrt{g}\left(a^{i}a^{j}-\frac{1}{2}g^{ij}a_{k}a^{k}\right)-\nabla_{k}\left[2\pi^{k(i}N^{j)}-\pi^{ij}N^{k}\right]$
$\displaystyle+\beta\sqrt{g}\left[\nabla^{(i}\nabla^{j)}N-g^{ij}\nabla^{2}N\right]+\mu\pi^{ij}\,.\qquad$

##  3 Energy in Hořava-Lifshitz gravity

If we calculate the Hamiltonian (5) in the constrained submanifold it reduces
to a surface term. We are going to show in this section that the physical
Hamiltonian, that is, the Hamiltonian expressed solely in terms of the
physical degrees of freedom is the mentioned surface term. To do so, we
evaluate the Lagrangian on the constrained submanifold and express the kinetic
terms solely in terms of the physical degrees of freedom. From it we obtain
the physical Hamiltonian which ends up being the surface term $E$. This is a
way to show that the gravitational energy is given by the surface term $E$. An
interesting point in our argument will be that no integration by parts needs
to be performed.

In order to identify which field components propagate and which are static
ones, we will use an orthogonal linear decomposition in transverse and
longitudinal parts, the ADM decomposition Arnowitt1961 ; ADM2008 . A symmetric
tensor that vanishes at infinity can be expressed as

$f_{ij}=f_{ij}^{TT}+f_{ij}^{T}+2\partial_{(j}f_{i)}\,.$ (16)

The transverse part
$f_{ij}^{T}\equiv\frac{1}{2}[\delta_{ij}f^{T}-\frac{1}{{\Delta}}\partial_{i}\partial_{j}f^{T}]$
is divergence-free. The $f_{ij}^{TT}$-part is divergence-free and trace-free.
The remaining term $2\partial_{(j}f_{i)}$ is its longitudinal part.
$f^{T}=\delta^{ij}f_{ij}^{T}$ is the trace of the transverse part of $f_{ij}$.
$\frac{1}{\Delta}$ is the inverse of the flat space Laplacian, defined on the
space of functions which vanish at infinity.

It is known that the Hamiltonian of the Hořava-Lifshitz action can be
rewritten as a linear combination of constraints Bellorin_Restuccia_2016B
provided two surface terms are included:

$E=-\beta\oint_{\partial\Sigma_{t}}dS_{i}g^{T}_{,i}-2\alpha\oint_{\partial\Sigma_{t}}dS_{i}N_{,i}\,\,,$
(17)

and that the physical degrees of freedom of the Hořava-Lifshitz gravity at the
linearized level are the TT tensorial modes. Notice that the expression of
$g^{T}$ and $N$ in terms of the physical degrees of freedom can be obtained
from the two second class constraints (13) and (12), respectively.

The Lagrangian evaluated on the submanifold of constraints reduces to the
kinetic terms plus surface terms

$L=\int dtd^{3}x\,\pi^{ij}\,\partial_{t}g_{ij}-E\,.$ (18)

We assume the following flat asymptotic behaviour, as in GR:

$\displaystyle g_{ij}-\delta_{ij}=\mathcal{O}(1/r)\,,\quad\partial
g_{ij}=\mathcal{O}(1/r^{2})\,,$ (19)
$\displaystyle\pi^{ij}=\mathcal{O}(1/r^{2})\,,\quad\partial\pi^{ij}=\mathcal{O}(1/r^{3})\,,$
(20) $\displaystyle N-1=\mathcal{O}(1/r)\,,\quad\partial
N=\mathcal{O}(1/r^{2})\,,$ (21) $\displaystyle
N_{i}=\mathcal{O}(1/r)\,,\quad\partial N_{i}=\mathcal{O}(1/r^{2})\,.$ (22)

Furthermore, we can fix, using the spacelike diffeomorphisms
$\mathcal{F}_{Diff}$, the coordinate condition ADM2008

$g_{i}=x^{i}+\left(\frac{1}{4\Delta}\right)g^{T}_{,i}\,,$ (23)

or in differential form

$g_{ij,jkk}-\frac{1}{4}g_{kj,kji}-\frac{1}{4}g_{jj,kki}=0\,.$ (24)

We will use an equivalent decomposition to (16) but reorganized in a different
way. We consider the following orthogonal decomposition of $h_{ij}\equiv
g_{ij}-\delta_{ij}$,

$h_{ij}=h_{ij}^{T\tau}+h_{ij}^{\tau}+\frac{1}{3}\delta_{ij}h,$ (25)

where $\delta_{ij}h_{ij}^{T\tau}=\delta_{ij}h_{ij}^{\tau}=0$, ’$\tau$’ means
traceless with respect to $\delta_{ij}$. It follows that $h=h_{ii}$. By
definition, $\partial_{i}h_{ij}^{T\tau}=0$, ’T’ means transverse with respect
to $\partial_{i}\,$, and $h^{\tau}_{ij}$ is defined in terms of a vector field
$W_{j}$

$h_{ij}^{\tau}=\partial_{i}W_{j}+\partial_{j}W_{i}-\frac{2}{3}\delta_{ij}\partial_{k}W_{k}.$
(26)

This decomposition exists and it is unique. It is analogous to the York
decomposition, however it is not covariant. It can also be rewritten as

$\left(1+\frac{1}{3}h\right)h_{ij}=h_{ij}^{T\tau}+h_{ij}^{\tau}+\frac{1}{3}g_{ij}h.$
(27)

Comparing the decomposition (25) with the T+L ADM decomposition (16) we obtain

$h_{ij}^{T\tau}=h_{ij}^{TT}\,,\quad
W_{i}=h_{i}-\frac{1}{4}\frac{\partial_{i}h^{T}}{\Delta}\,,\quad
h=h_{ii}=h^{T}+2\partial_{i}h_{i}\,,$ (28)

where $h_{ij}^{TT}\,$, $h^{T}$ and $h_{i}$ are the T+L components of the ADM
decomposition. Using the $\mathcal{F}_{Diff}$ symmetry we can always fix
$W_{i}=0$. In fact, the gauge fixing condition is exactly (23). We then have

$h_{ij}=\frac{h_{ij}^{T\tau}}{1+\frac{1}{3}h}+\frac{1}{3}g_{ij}\frac{h}{1+\frac{1}{3}h}.$
(29)

In order to analyse the reduction to the Hamiltonian density in terms of the
physical degrees of freedom, the transverse traceless modes, we consider now
the covariant decomposition introduced by York york1973

$\pi^{ij}=\tilde{\pi}^{ijT\tau}+\tilde{\pi}^{ij\tau}+\frac{1}{3}g^{ij}\tilde{\pi},$
(30)

where $g_{ij}\tilde{\pi}^{ijT\tau}=g_{ij}\tilde{\pi}^{ij\tau}=0\,$,
$\nabla_{i}\tilde{\pi}^{ijT\tau}=0$ and $\tilde{\pi}^{ij\tau}$ is defined in
terms of a vector field $U^{j}$

$\tilde{\pi}^{ij\tau}=\nabla^{i}U^{j}+\nabla^{j}U^{i}-\frac{2}{3}g^{ij}\nabla_{k}U^{k}.$
(31)

It then follows from the constraints of the theory that

$\pi^{ij}=\tilde{\pi}^{ijT\tau}.$ (32)

We remark that $\tilde{\pi}^{ijT\tau}$ is a function of the $\pi^{ijTT}$
transverse-traceless ADM modes which are independent of the metric $g_{ij}$.
One way of making explicit this point is to consider the sequence
$\pi_{n}^{ij}\,$, $n=0,1,2,...,$ where $\pi_{0}^{ij}=\pi^{ijTT}$ and

$\partial_{i}\pi_{n+1}^{ij}+\Gamma_{ik}^{j}\pi_{n}^{ik}=0\,,\quad\delta_{ij}\pi_{n+1}^{ij}+h_{ij}\pi_{n}^{ij}=0.$
(33)

If this sequence is convergent, that is if there exists a fixed point of the
procedure, we obtain $\pi_{n}^{ij}\rightarrow\tilde{\pi}^{ijT\tau}$, which
depends on the metric $g_{ij}$ and $\pi_{ij}^{TT}$ the transverse-traceless
modes of the ADM decomposition. We argue now that the Hamiltonian density in
terms of physical degrees of freedom can be obtained in terms of
$\tilde{\pi}^{ijT\tau}$ and $h_{ij}^{TT}$

$\pi^{ijTT}\,,\,h_{kl}^{TT}\quad\longleftrightarrow\quad\tilde{\pi}^{ijT\tau}\left(h_{mn}^{TT},\pi_{pq}^{TT}\right)\,,\,\,h_{kl}^{TT}.$
(34)

The Lagrangian of the Hořava-Lifshitz theory at the kinetic-conformal point
evaluated on the constrained submanifold is given by (18), where the surface
terms $E$, (17), arise from the Fréchet differentiability of the action and
the expression of the Hamiltonian rewritten in terms of the constraints.

We then have, using that $g_{ij}\pi^{ij}=0$ is a constraint,

$\displaystyle\pi^{ij}\partial_{t}h_{ij}$ $\displaystyle=$
$\displaystyle\pi^{ij}\partial_{t}\left(\frac{h_{ij}^{T\tau}}{1+\frac{1}{3}h}+\frac{1}{3}g_{ij}\frac{h}{1+\frac{1}{3}h}\right)\,$
(35) $\displaystyle=$
$\displaystyle\pi^{ij}\partial_{t}\left(\frac{h_{ij}^{T\tau}}{1+\frac{1}{3}h}\right)+\pi^{ij}\partial_{t}g_{ij}\frac{h}{3\left(1+\frac{1}{3}h\right)}\quad.$

Also $\partial_{t}g_{ij}=\partial_{t}h_{ij}\,$, hence

$\pi^{ij}\partial_{t}h_{ij}=\left[\left(1+\frac{1}{3}h\right)\tilde{\pi}^{ijT\tau}\right]\partial_{t}\left[\frac{h_{ij}^{T\tau}}{1+\frac{1}{3}h}\right]\,,$
(36)

where, in the gauge $W_{i}=0\,$, $1+\frac{1}{3}h=1+\frac{1}{2}h^{T}\,$,
$h^{T}$ is the transverse component in the ADM decomposition and can be
obtained from the constraints of the theory in terms of $h_{ij}^{TT}$. We then
conclude that the canonical conjugate variables are

$\frac{h_{ij}^{TT}}{1+\frac{1}{2}h^{T}}\,,\quad\left(1+\frac{1}{2}h^{T}\right)\tilde{\pi}^{ijT\tau}.$
(37)

The value of the Hamiltonian on the constrained submanifold is then given by
the two surface terms in (17). The Lagrangian is

$L=\int
dtd^{3}x\left[(1+\frac{1}{2}h^{T})\tilde{\pi}^{ijT\tau}\right]\partial_{t}\left[\frac{h_{ij}^{TT}}{1+\frac{1}{2}h^{T}}\right]-E,$
(38)

where $E$ is the surface term (17). In contrast with the analysis in ADM2008 ,
in our argument there is no integration by parts.

##  4 The gravitational fields in the wave zone

We define the wave zone, as in Arnowitt1961 ; ADM2008 , by the following three
conditions: First, $kr\gg 1$ where $k$ is the wave number and $r$ is radial
distance. This condition is the same that defines the far-zone in linear
theories such as classical electrodynamics, and can be satisfied if the radial
distance is far enough from the sources. In addition to the previous
condition, for non-linear theories it is necessary to impose more restrictive
conditions in order to ensure that the self-interaction does not destroy the
free propagation of the dynamical modes. Then as a second condition, we demand
that the deviations of the fields from ”flat background” are of the order of
$\mathcal{O}(1/r)$, i.e.
$|g_{ij}-\delta_{ij}|\sim|N-1|\sim|N_{i}|=\mathcal{O}(A/r)\ll 1$, where
$A(t,\theta,\phi)$ represent generic functions of time, and angular
coordinates such that $A$ and all its derivatives are bounded. The third
condition is $|\partial g/\partial(kr)|^{2}\sim|\partial N/\partial(kr)|^{2}$
$\sim|\partial N_{i}/\partial(kr)|^{2}\ll|g-\delta|$, which is necessary to
guarantee that the interaction of sub-leading order modes can not interfere to
leading order $\mathcal{O}(1/r)$.

The wave zone in GR is a region of space far away from the source where the TT
modes of the metric, the physical degrees of freedom, propagate freely
according to a linear wave equation. On the near zone the interaction terms
are relevant and the propagation of the TT modes is not governed by the linear
wave equation. The T mode of the metric and lapse $N$ are nontrivial in the
wave zone and contribute to the gravitational energy and momentum. Beyond the
wave zone the dominant terms are the T mode and the lapse $N$ .

The solution of primary constraints (11) implies that the propagating parts of
the momenta in the wave zone behaves as

$\displaystyle\pi^{ijT}\sim\pi^{ijL}\sim\frac{B^{ij}}{r^{2}}+k\frac{\hat{A}^{ij}e^{ikr}}{r^{2}},$
(39)
$\displaystyle\pi^{ijTT}\sim\frac{B^{ij}}{r^{2}}+k\frac{\hat{A}^{ij}e^{ikr}}{r}\,,$
(40)

we use $A$ y $B$ as generic tensorial functions of angles such that they and
their derivatives are bounded.

If we use (39) and (40) in the secondary constrains (12) and (13) and the
transverse gauge $g_{ij,j}=0$, we obtain a coupled system of second-order
elliptical partial differential equations for variables $g^{T}$ and $N$. We
use here a different gauge condition than in the previous section. A more
adequate one for the present analysis. This is admissible since the
gravitational energy, the surface term $E$, obtained in the previous section
is independent of the gauge condition. If we multiply (12) by two and add the
result to (13) we decoupled the system and, since $\beta\neq\frac{\alpha}{2}$
for the values of $\beta$ and $\alpha$ determined from experimental data, we
can estimate $N$ and $g^{T}$ at low energies in the wave zone:

$N-1\sim g^{T}\sim\frac{B}{r}+\frac{\hat{A}e^{ikr}}{r^{2}}\,,$ (41)

these fields have a nontrivial Newtonian part $\mathcal{O}(1/r)$. They are
zero in a linearized version of the theory.

In the complete version of the theory, where higher order derivative terms are
included, there are oscillatory zero modes of the same order which have to be
eliminated to obtain consistency of the theory. We do not have this problem in
the low energy regime.

In the low energy regime, where only the contribution (8) to the potential is
considered, the term $\alpha a_{i}a^{i}$ breaks manifestly the relativistic
symmetry. Its contribution to the present analysis can be determined directly
from (41). We obtain

$\displaystyle
a_{i}\lesssim\frac{B_{i}}{r^{2}}+k\frac{\hat{A}_{i}e^{ikr}}{r^{2}},$ (42)
$\displaystyle
a_{i}a^{i}\lesssim\frac{B}{r^{4}}+k^{2}\frac{\hat{A}e^{ikr}}{r^{4}},$ (43)
$\displaystyle\nabla^{i}a_{i}\lesssim\frac{B}{r^{3}}+k^{2}\frac{\hat{A}e^{ikr}}{r^{2}}\,,$
(44)

although terms including the vector $a_{i}$ are not involved in the field
equations at order $\mathcal{O}\left(1/r\right)$, they do contribute to the
gravitational energy and to the Newtonian background.

From the transverse and longitudinal decomposition of the dynamical equations
(14) and (15), we obtain the canonical form of the wave equation

$\displaystyle\partial_{t}g_{ij}^{TT}=2\pi_{ij}^{TT}+\mathcal{O}(1/r^{2}),$
(45) $\displaystyle\partial_{t}\pi^{ijTT}=\frac{1}{2}\beta\Delta
g_{ij}^{TT}+\mathcal{O}(1/r^{2}),$ (46)

equivalently,

$\partial_{t}^{2}g^{ijTT}-\beta\Delta g^{ijTT}=0+\mathcal{O}(1/r^{2}).$ (47)

Then to the leading order $\mathcal{O}(1/r)$ the transverse traceless
components of the spatial metric satisfy a wave equation with speed of
propagation $\sqrt{\beta}$. Detection from gravitational waves arising from
the merge of the neutron star binary system GW170817 Abbott2017  and its
electromagnetic counterpart, ${\gamma}$-ray burst GBR170817A Abbott_2017GRB ,
restrict the space of $\beta$-parame-ter to $|1-\sqrt{\beta}|\leq 10^{-15}$
Abbott2019b . Then, if that is so, the prediction on gravitational waves of
Hořava-Lifshitz theory at the kinetic-conformal point, at low energies, is the
same as in GR.

##  5 Discussion and conclusions

We showed that in Hořava-Lifshitz theory at the kinetic-conformal point, in
the low energy regime, a wave zone can be consistently defined. In it the
physical degrees of freedom, which reduce to the transverse traceless
tensorial modes, satisfy a linear wave equation. The same one that arises from
a linear perturbative approach Blas2011 ; BellorinRestuccia2018 , but unlike
this there are Newtonian non-trivial contributions, that do not obstruct the
free propagation (radiation) of the physical degrees of freedom. Among these
contributions there are terms which manifestly break the relativistic
symmetry. These terms which determine a different physical behavior, for
example of the static, spherically symmetric solutions of the Hořava-Lifshitz
gravity theory compared to GR, do not contribute to the free propagation of
the physical degree of freedom to the dominant order in the wave zone.
However, they provide a relevant contribution to the gravitational energy and
the gravitational momentum.

The gravitational energy of the Hořava-Lifshitz gravity has been considered in
Blas2011  from an asymptotic analysis, however the canonical reduction of the
Hamiltonian to the physical degrees of freedom has not been analyzed. In GR
the reduction to a canonical formulation in terms of the TT modes can be
achieved by a suitable gauge fixing condition. The main point is that the
gauge fixing procedure for the Hořava-Lifshitz gravity reduces to a spacelike
coordinate election, since only reparametrization on time is allowed. We
obtained in section 3 the canonical reduction in terms of the true physical
degrees of freedom, which at linearized level reduce to the TT modes of the
metric. In general they are well defined functions of the TT modes of the
metric.

We showed, using the $\mathcal{F}_{Diff}$ symmetry, that the Hořava-Lifshitz
theory, at the kinetic-conformal point, can be canonically reduced to a
Hamiltonian formulation in terms of the physical degrees of freedom, in a
particular coordinate frame. They reduce to the TT tensorial modes in the wave
zone.

In Hořava-Lifshitz theory there does not exist a universal (scale invariant)
constant as the light velocity, however the energy-dependent coupling constant
$\sqrt{\beta}$, in the renormalization flow from the UV regime to the IR
point, should end up having a value very near or equal to the speed of light.
In that case, although the Hořava-Lifshitz theory breaks the relativistic
symmetry, the wave equation coincides with the relativistic one arising in GR.

We expect that the interaction terms with high spatial derivatives modify the
wave equation by introducing linear high order spatial derivatives in the wave
zone, i.e $\Delta^{2}$ and $\Delta^{3}$ operators in the propagating equation.
However, it is unknown if in this case the Newtonian background interacts in a
non-trivial way with the propagation of the physical degrees of freedom. In
particular, the resolution of the constraints is in this case a non-trivial
problem.

##  6 Aknowledgments

J. Mestra-Páez acknowledge financial support from Beca Doctorado Nacional 2019
CONICYT, Chile. N° BECA: 21191442. J. M. P is supported by the projects
ANT1956 and SEM18-02 of the University of Antofagasta, Chile.

## References

  * (1) M. C. Miller, N. Yunes, The new frontier of gravitational waves, Nature 568 (7753) (2019) 469–476.  doi:10.1038/s41586-019-1129-z.
  * (2) A. Nitz, C. Capano, S. Kumar, Y.-F. Wang, S. Kastha, M. Schäfer, R. Dhurkunde, M. Cabero, 3-ogc: Catalog of gravitational waves from compact-binary mergers (2021).  arXiv:2105.09151.
  * (3) P. Hořava, Quantum gravity at a Lifshitz point, Phys. Rev. D 79 (2009) 084008.  doi:10.1103/PhysRevD.79.084008.
  * (4) D. Blas, O. Pujolàs, S. Sibiryakov, Consistent extension of Hořava gravity, Phys. Rev. Lett. 104 (2010) 181302.  doi:10.1103/PhysRevLett.104.181302.
  * (5) A. Emir Gümrükçüoğlu, M. Saravani, T. P. Sotiriou, Hořava gravity after GW170817, Phys. Rev. D 97 (2018) 024032.  doi:10.1103/PhysRevD.97.024032.
  * (6) J. Bellorín, A. Restuccia, A. Sotomayor, Consistent Hořava gravity without extra modes and equivalent to general relativity at the linearized level, Phys. Rev. D 87 (2013) 084020.  doi:10.1103/PhysRevD.87.084020.
  * (7) D. Blas, O. Pujolas, S. Sibiryakov, Comment on “strong coupling in extended Hořava-Lifshitz gravity”[phys. lett. B 685 (2010) 197], Physics Letters B 688 (4-5) (2010) 350–355.  doi:10.1016/j.physletb.2010.03.073.
  * (8) R. Arnowitt, S. Deser, C. W. Misner, Republication of: The dynamics of general relativity, General Relativity and Gravitation 40 (9) (2008) 1997–2027.  doi:10.1007/s10714-008-0661-1.
  * (9) P. Hořava, C. M. Melby-Thompson, General covariance in quantum gravity at a lifshitz point, Phys. Rev. D 82 (2010) 064027\.  doi:10.1103/PhysRevD.82.064027.
  * (10) A. Wang, Y. Wu, Cosmology in nonrelativistic general covariant theory of gravity, Phys. Rev. D 83 (2011) 044031.  doi:10.1103/PhysRevD.83.044031.
  * (11) A. M. da Silva, An alternative approach for general covariant Hořava–Lifshitz gravity and matter coupling, Classical and Quantum Gravity 28 (5) (2011) 055011.  doi:10.1088/0264-9381/28/5/055011.
  * (12) J. Klusoň, Hamiltonian analysis of nonrelativistic covariant restricted-foliation-preserving diffeomorphism invariant Hořava-Lifshitz gravity, Phys. Rev. D 83 (2011) 044049.  doi:10.1103/PhysRevD.83.044049.
  * (13) K. Lin, A. Wang, Q. Wu, T. Zhu, Strong coupling in nonrelativistic general covariant theory of gravity, Phys. Rev. D 84 (2011) 044051.  doi:10.1103/PhysRevD.84.044051.
  * (14) T. Zhu, Q. Wu, A. Wang, F.-W. Shu, U(1) symmetry and elimination of spin-0 gravitons in Horava-Lifshitz gravity without the projectability condition, Phys. Rev. D 84 (2011) 101502.  doi:10.1103/PhysRevD.84.101502.
  * (15) J. Bellorín, A. Restuccia, Quantization of the Hořava theory at the kinetic-conformal point, Phys. Rev. D 94 (2016) 064041\.  doi:10.1103/PhysRevD.94.064041.
  * (16) C. Charmousis, G. Niz, A. Padilla, P. M. Saffin, Strong coupling in Hořava gravity, Journal of High Energy Physics 2009 (08) (2009) 070.  doi:10.1088/1126-6708/2009/08/070.
  * (17) M. Visser, Power-counting renormalizability of generalized Hořava gravity (2009).  arXiv:0912.4757.
  * (18) A. Papazoglou, T. P. Sotiriou, Strong coupling in extended Hořava-Lifshitz gravity, Physics Letters B 685 (2-3) (2010) 197–200.  doi:10.1016/j.physletb.2010.01.054.
  * (19) D. Orlando, S. Reffert, The renormalizability of Hořava-Lifshitz-type gravities, Classical and Quantum Gravity 26 (15) (2009) 155021.  doi:10.1088/0264-9381/26/15/155021.
  * (20) F.-W. Shu, Y.-S. Wu, Stochastic quantization of the Hořava gravity (2009).  arXiv:0906.1645.
  * (21) D. Benedetti, F. Guarnieri, One-loop renormalization in a toy model of Hořava-Lifshitz gravity, Journal of High Energy Physics 2014 (3) (2014) 78.  doi:10.1007/JHEP03(2014)078.
  * (22) A. Contillo, S. Rechenberger, F. Saueressig, Renormalization group flow of Hořava-Lifshitz gravity at low energies, Journal of High Energy Physics 2013 (12) (2013) 17.  doi:10.1007/JHEP12(2013)017.
  * (23) G. D’Odorico, F. Saueressig, M. Schutten, Asymptotic freedom in Hořava-Lifshitz gravity, Physical Review Letters 113 (17) (2014) 171101\.  doi:10.1103/PhysRevLett.113.171101.
  * (24) G. D’Odorico, J.-W. Goossens, F. Saueressig, Covariant computation of effective actions in Hořava-Lifshitz gravity, Journal of High Energy Physics 2015 (10) (2015) 1–30.  doi:10.1007/JHEP10(2015)126.
  * (25) A. O. Barvinsky, D. Blas, M. Herrero-Valea, S. M. Sibiryakov, C. F. Steinwachs, Renormalization of Hořava gravity, Physical Review D 93 (6) (2016) 064022\.  doi:10.1103/PhysRevD.93.064022.
  * (26) A. Wang, Hořava gravity at a Lifshitz point: a progress report, International Journal of Modern Physics D 26 (07) (2017) 1730014.  doi:10.1142/S0218271817300142.
  * (27) S. Shin, M.-I. Park, On gauge invariant cosmological perturbations in UV-modified Hořava gravity, Journal of Cosmology and Astroparticle Physics 2017 (12) (2017) 033.  doi:10.1088/1475-7516/2017/12/033.
  * (28) M. Pospelov, Y. Shang, Lorentz violation in Hořava-Lifshitz-type theories, Physical Review D 85 (10) (2012) 105001.  doi:10.1103/PhysRevD.85.105001.
  * (29) R. Arnowitt, S. Deser, C. W. Misner, Wave zone in General Relativity, Phys. Rev. 121 (1961) 1556–1566.  doi:10.1103/PhysRev.121.1556.
  * (30) D. Blas, H. Sanctuary, Gravitational radiation in Hořava gravity, Phys. Rev. D 84 (2011) 064004.  doi:10.1103/PhysRevD.84.064004.
  * (31) A. Restuccia, F. Tello-Ortiz, Pure electromagnetic-gravitational interaction in Hořava-Lifshitz theory at the kinetic conformal point, Eur. Phys. J. C (80) (2020) 86.  doi:10.1140/epjc/s10052-020-7674-7.
  * (32) J. Bellorín, A. Restuccia, F. Tello-Ortiz, Anisotropic coupling of gravity and electromagnetism in Hořava-Lifshitz theory, Phys. Rev. D 98 (2018) 104018.  doi:10.1103/PhysRevD.98.104018.
  * (33) A. Restuccia, F. Tello-Ortiz, Quantum aspects of the gravitational-gauge vector coupling in the Hořava-Lifshitz theory at the kinetic conformal point (2020).  arXiv:2012.09354.
  * (34) T. Regge, C. Teitelboim, Role of surface integrals in the Hamiltonian formulation of general relativity, Annals of Physics 88 (1) (1974) 286 – 318\.  doi:10.1016/0003-4916(74)90404-7.
  * (35) J. W. York Jr, Conformally invariant orthogonal decomposition of symmetric tensors on Riemannian manifolds and the initial-value problem of General Relativity, Journal of Mathematical Physics 14 (4) (1973) 456–464.  doi:10.1063/1.1666338.
  * (36) B. P. Abbott, et al., GW170817: Observation of gravitational waves from a binary neutron star inspiral, Phys. Rev. Lett. 119 (2017) 161101.  doi:10.1103/PhysRevLett.119.161101.
  * (37) B. P. Abbott, et al., Multi-messenger observations of a binary neutron star merger, The Astrophysical Journal 848 (2) (2017) L12.  doi:10.3847/2041-8213/aa91c9.
  * (38) B. P. Abbott, et al., Tests of General Relativity with GW170817, Phys. Rev. Lett. 123 (2019) 011102.  doi:10.1103/PhysRevLett.123.011102.
  * (39) J. Bellorín, A. Restuccia, Einstein’s quadrupole formula from the kinetic-conformal Hořava theory, Int. Mod. Phys. D 27 (2018) 1750174.  doi:10.1142/S0218271817501747.