# Comparative Analysis of Fluctuations in Viscoelastic Stress: A Comparison of
the Temporary Network and Dumbbell models

Arturo Winters  Hans Christian Öttinger  Department of Materials, ETH Zurich,
Vladimir-Prelog-Weg 5, Zurich 8093, Switzerland<EMAIL_ADDRESS>Jan
Vermant  Department of Materials, ETH Zurich, Vladimir-Prelog-Weg 5, Zurich
8093, Switzerland

###### Abstract

Abstract. Traditionally, stress fluctuations in flowing and deformed materials
are overlooked, with an obvious focus on average stresses in a continuum
mechanical approximation. However, these fluctuations, often dismissed as
“noise”, hold the potential to provide direct insights into the material
structure and its structure-stress coupling, uncovering detailed aspects of
fluid transport and relaxation behaviors. Despite advancements in experimental
techniques allowing for the visualization of these fluctuations, their
significance remains largely untapped, as modeling efforts continue to target
Newtonian fluids within the confines of Gaussian noise assumptions. In the
present work a comparative analysis of stress fluctuations in two distinct
microstructural models is carried out: the temporary network model and the
dumbbell model. Despite both models conforming to the Upper Convected Maxwell
Model at a macroscopic level, the temporary network model predicts non-
Gaussian fluctuations. We find that stress fluctuations within the temporary
network model exhibit more pronounced abruptness at local scale, with only an
enlargement of the control volume leading to a gradual Gaussian-like noise,
diminishing the differences between the two models. These findings underscore
the heightened sensitivity of _fluctuating rheology_ to microstructural
details and the microstructure-flow coupling, beyond what is captured by
macroscopically averaged stresses.

††preprint: AIP/123-QED

##  I Introduction

In 1827, BrownBrown (1828) discovered Brownian motion, a phenomenon that lies
at the heart of understanding transport phenomena in soft matter and motivated
Einstein’s doctoral thesis to look at fluctuations.Einstein (1905) Einstein
proposed a way to determine the size of molecules in a liquid which led to a
method for calculating Avogadro’s number, laying the groundwork for modern
physical chemistry and statistical physics. The present paper presents
insights into thermal fluctuations, but now in flowing complex materials,
using the viscoelastic Maxwell model as a fundamental rheological example. In
the field of material science, thermal fluctuations at equilibrium are crucial
for interpreting scattering experiments Zárate and Sengers (2006); Larsen and
Furst (2008); Venerus and Öttinger (2018) or and microrheological
experiments,Einstein (1905); Mason and Weitz (1995); Furst and Squires (2017)
providing a deeper understanding of material behavior at the microscopic
level. Furthermore, many models of complex materials such as bead-spring
modelsÖttinger (1996) are based on fluctuations at a microscopic scale. To
date, the role of Brownian motion in transport processes at nanoscale remains
an active field of research.Marbach, Dean, and Bocquet (2018); Marbach, Zheng,
and Holmes-Cerfon (2022)
The continuum mechanical modeling of fluctuations at equilibrium, rely on the
coarse grained description of Landau and Lipshitz,Landau and E.M. (1959) which
treats dissipative fluxes $\boldsymbol{J}$, e.g. the heat flux or the extra
stress tensor $\boldsymbol{\tau}$, as stochastic variables. The latter
quantity is additively split into an average value and a fluctuating part
$\delta\boldsymbol{J}$. The fluctuations at equilibrium fluctuation are
assumed to be random and uncorrelated at different points in space and moments
in time. The second moment of $\delta\boldsymbol{J}$ is hence delta-correlated
in space and time

$\langle\delta\boldsymbol{J}_{\alpha}\delta\boldsymbol{J}_{\beta}\rangle=2k_{\text{B}}TK_{\alpha\beta}\delta(\boldsymbol{r}-\boldsymbol{r^{\prime}})\delta(t-t^{\prime}),$
(1)

where $k_{\text{B}}$ is the Boltzmann constant, $T$ the temperature, and
$K_{\alpha\beta}$ are components of the so called Onsager dissipation matrix,
which are reflecting material properties such as the thermal conductivity or
the viscosity. The delta-correlation in time reflects the Gaussian noise, when
the system is at equilibrium. When we assume local but not global equilibrium,
we can develop fluctuating hydrodynamics, which successfully models
fluctuations in flows or with thermal gradients, in Newtonian fluids.Zárate
and Sengers (2006); Kirkpatrick, Ortiz de Zárate, and Sengers (2013) So far
previous work has mainly focused on Newtonian materials, in the present work
we turn to complex fluids, where the rheological properties of the fluids are
more complicated, and the fluids are described by either a structural variable
or “memory” effects. Extending this is motivated by recent experimental
developments as, for instance, by high speed confocal counterrotating rheo-
microscopyColombo  _et al._ (2019) where particle tracking methods can be
applied to flowing systems. An example is shown in fig. 1, where the tracks of
micrometer sized particles fluctuate during flow of a single relaxation time
(near Maxwellian) wormlike micellar solution characterized in Snijkers et
al.Snijkers  _et al._ (2009) The tracks show the convective motion of the
particles with superimposed noise in the displacement, resulting from the
combined effect of stress fluctuations on the particle.

Figure 1: Tracked 1-µm beads near the stagnation plane ($v,\,\omega$ for the
velocity and vorticity directions, respectively) in a shear flow at rate
$\dot{\gamma}=1\text{s}^{-1}$ wormlike micellar solution observed by high
speed rheo-confocal microscopy.Colombo  _et al._ (2019); Ershov  _et al._
(2022) The (near Maxwellian) solution at 25° characterized in Snijkers  _et
al._ (2009) has a single relaxation time $\lambda=1.5$ s and a modulus of
$G=33.4$ Pa. The tracks show a convective motion along the flow direction with
superimposed fluctuations.

The present works seeks to understand whether the fluctuations of the stress
in complex fluids can be used to infer more detailed information on
microscopic dynamics. The few works dealing with the topic,Hütter, Olmsted,
and Read (2020); Öttinger and Grmela (1997) as we explore further in sec.
II.2, confined themselves to model fluctuations of the Gaussian type as in eq.
(1). Our initial objective is to advance the concept of fluctuating rheology.
We begin by demonstrating that, even in the context of the simplest
viscoelastic model, the Upper Convected Maxwell Model (UCM), Gaussian
approximations are not universally applicable. Our second objective
specifically addresses the nuances of the UCM, adhering to the principle that
“the main value of this model is as a guide for qualitative thinking”.Ferry
(1980) Remarkably, both the Temporary Network Model (TNM) and the Hydrodynamic
Dumbbell Model (HDM), which serve as basic constitutive model systems for two
distinct groups of microscopic models, reproduce the Maxwellian
behavior.Larson (1999); Öttinger (1996) Till now, literature emphasized the
equality of the two models when describing the stress evolution,Ferry (1980);
Morrison (2001); Venerus and Öttinger (2018); Larson (1999) which is
intriguing given their fundamental differences on a microscopic scale. In this
paper we demonstrate how differences do come out when considering the stress
fluctuations and suggest methods to measure and use them.

The paper is organized as follows: First (sec. II.1), we describe the
specifics of the two approaches and how they are cross-grained to yield the
UCM. Thereafter, we introduce a thermodynamic framework to express
fluctuations (sec. II.2) and focus on the Fluctuation-Dissipation Theorem
(FDT) and its role in connecting the microscopic world with the macroscopic
one.Öttinger, Peletier, and Montefusco (2021) From sec. III onward we present
novel results. We start with the TNM and HDM dissipation potentials (sec.
III). The differences of the two models with respect to the probability
density $p$ and to $\boldsymbol{\tau}$ follow in sec. IV and V. In sec. VI we
present possible ways to measure them. Further, we underscore the significance
of fluctuating rheology to better understand micro structure - flow coupling.
In sec. VII, we synthesize our results (VII).

##  II Background

Figure 2: Linear viscoelastic moduli of an unentangled polymer melt of
polystyrene, using time temperature superposition with $a_{t}$ the horizontal
shift factor (data courtesy of Costanzo  _et al._ (2022)). A bead-spring Rouse
chain is used to fit the data, where the chain is composed of 4 connected
dumbbells.

Figure 3: Linear viscoelastic moduli of associative telechelic polymers,
replotted from Annable et al.Annable, Buscall, and Ettelaie (1996) The data is
fitted with a linear Maxwell model with one relaxation time. The dynamics of
this physical polymer network can be represented by a temporary network model.

Figure 4: In the dumbbell model, dumbbells vibrate and move through the
continuum volumes.

Figure 5: The temporary network model is composed of springs of a network that
are destroyed (and generated) with rate $\lambda^{-1}$.

###  II.1 Temporary network and dumbbell model

Maxwell introduced the first, most simple and widespread phenomenological
model for viscoelasticity in 1867, albeit for gases.Maxwell (1867) Later it
was enhanced by its observer-invariant formulation, known as the UCM,Oldroyd
(1950); Lodge (1974) which relates the stress tensor $\boldsymbol{\tau}$ to
relative deformations:

$\partial_{t}{\boldsymbol{\tau}}=\boldsymbol{\kappa}\cdot\boldsymbol{\tau}+\boldsymbol{\tau}\cdot\boldsymbol{\kappa}^{T}-G(\boldsymbol{\kappa}+\boldsymbol{\kappa}^{T})-\frac{4H}{\zeta}\boldsymbol{\tau},$
(2)

and its corresponding integral formulation, namely the Lodge equationLodge
(1974); Morrison (2001)

$\boldsymbol{\tau}=\frac{G}{\lambda}\int_{-\infty}^{t}e^{-(t-t^{\prime})/\lambda}\big{[}\boldsymbol{\delta}-\boldsymbol{C}^{-1}(t,t^{\prime})\big{]}\mathop{}\\!\mathrm{d}t^{\prime}.$
(3)

$G$ is the modulus, $\lambda$ the relaxation time,
$\boldsymbol{\kappa}=(\nabla\boldsymbol{u})^{T}$ is the transpose of the
velocity gradient tensor and $\boldsymbol{C}^{-1}$ the Finger strain tensor.
These two models can be found in rheology textbooks Macosko (1994); Larson
(1999); Venerus and Öttinger (2018); Phan-Thien (2002); Morrison (2001); Ferry
(1980) and reproduce the rate and time dependency of the stress evolution of
materials such as diluted polymer solutions (fig. 5) or networks of ssociative
polymers with weak physical cross links (fig. 5).

Two fundamentally different model systems, namely bead-spring dumbbells and
temporary networks, effectively represent the microscopic dynamics associated
with the macroscopic UCM.Rouse (1953); Öttinger (1996) In the HDM, a springs
with stiffness $H$ connects two beads at distance $\boldsymbol{Q}$ (fig. 5).
Under the effect of a flow field the spring is stretched while the rapid
interactions between the beads and the surrounding fluid are put into the
model as Gaussian noise, as described in equation (eq. (4)). In contrast, the
TNM is rooted in rubber elasticity.Green and Tobolsky (2004); Larson (1999)
The structural variable $\boldsymbol{Q}$ represents the end-to-end vector
between crosslinks (fig. 5). Strands within this network deform in an affine
manner and are subject to continuous cycles of breaking and reformation,
occurring at a consistent rate $\lambda_{\text{TNM}}$. The introduction of new
strands again ensures the adherence to a Gaussian equilibrium distribution
$p_{\text{eq}}=(2\pi\sigma^{2})^{-3/2}\exp{(-\boldsymbol{Q}^{2}/\sigma^{2})}$,
with $\sigma^{2}$ being the characteristic equilibrium length.

These models possess three levels of description. The first level describes
the evolution of $\boldsymbol{Q}$. The microscopic dynamics of HDM are
described by a stochastic differential equation:

$\mathop{}\\!\mathrm{d}\boldsymbol{Q}=(\boldsymbol{\kappa}\cdot\boldsymbol{Q}-\frac{2H}{\zeta}\boldsymbol{Q})\mathop{}\\!\mathrm{d}t+\sqrt{\frac{4k_{\text{B}}T}{\zeta}}\mathop{}\\!\mathrm{d}\boldsymbol{W},$
(4)

where the $\cdot$ represent a contraction of indices, $\zeta$ the friction
coefficient and $\boldsymbol{W}$ is a vector, where each entry is an
independent Wiener process.Öttinger (1996)
The TNM starts from a deterministic equation for the advection,

$\mathop{}\\!\mathrm{d}\boldsymbol{Q}=\boldsymbol{\kappa}\cdot\boldsymbol{Q}\mathop{}\\!\mathrm{d}t,$
(5)

and one for the probability density of a strand’s “lifetime” $t_{L}$Biller and
Petruccione (1990)

$\phi(t_{\text{L}})=\frac{\exp{(-t_{\text{L}}\lambda_{\text{TNM}}^{-1})}}{\lambda_{\text{TNM}}}.$
(6)

Second, we seek the evolution of the probability density $p(t,\boldsymbol{Q})$
(PDF). For the HDM we haveBird, Hassager, and Curtiss (2017)

$\partial_{t}p=-\partial_{\boldsymbol{Q}}\cdot\left((\boldsymbol{\kappa}\cdot\boldsymbol{Q}-\frac{2H}{\zeta}\boldsymbol{Q})p\right)+\frac{2k_{\text{B}}T}{\zeta}\partial_{\boldsymbol{Q}}\cdot\partial_{\boldsymbol{Q}}p,$
(7)

where $\partial_{t},\partial_{\boldsymbol{Q}}$ represent the partial
differential operator with respect to time or the coordinates of the connector
vector, respectively. For the TNM it is given by :

$\partial_{t}p+\partial_{\boldsymbol{Q}}\cdot(\boldsymbol{\kappa}\cdot\boldsymbol{Q}p)=\lambda_{\text{TNM}}^{-1}[p_{\text{eq}}(\boldsymbol{Q})-p(\boldsymbol{Q},t)].$
(8)

Finally, to further coarse grain we extract the evolution equation of the 2nd
moment of $\boldsymbol{Q}$, namely the conformation tensor
$\boldsymbol{c}(\boldsymbol{x},t)=\langle\boldsymbol{Q}\boldsymbol{Q}\rangle=\int\boldsymbol{Q}\boldsymbol{Q}p\mathop{}\\!\mathrm{d}^{3}Q$,
from (7) and (8). The latter, deterministic equation is quantitatively equal
for both the TNM and the HDM for the identities
$\lambda=\zeta/4H=\lambda_{\text{TNM}}$, $\sigma^{2}=k_{\text{B}}T/H$ and
$G=nk_{\text{B}}T$, where $n$ represents the strand and dumbbell densities,
respectively. Kramers’ relation formalizes the momentum transmission over
continuum volumes performed by the springs as follows:

$\boldsymbol{c}(\boldsymbol{x},t)=nH\langle\boldsymbol{Q}\boldsymbol{Q}\rangle=nk_{\text{B}}T\boldsymbol{\delta}-\boldsymbol{\tau}(\boldsymbol{x},t).$
(9)

This connection links the evolution equation for the configuration tensor to
the stress and ultimately leads to the macroscopic version of the UCM (eq.
(2)), where $\boldsymbol{c}$’s dependence on the location coordinate
$\boldsymbol{x}$ is due to $\boldsymbol{\kappa}(\boldsymbol{x},t)$.

A key question we seek to address is: At which level of description and up to
which length scales do the two models differ? This deeper dive into the HDM
and the TNM shows that despite the different micro structural origin of the
evolution of the PDF, the second order moments evolve in the same manner, so
that the stress follows the UCM (or the integral version , the Lodge model).
However, a closer look at the stress fluctuations up to some lengthscale will
reveal the microscopic differences of the models. This requires us however, to
first introduce a suitable formulation of the FDT.

###  II.2 Relating fluctuations and dissipation

The FDT, first formulated in 1928 by Nyquist (1928), creates a connection
across length scales. It links the behavior of fluctuating microscopic world
or models to the dissipative part of associated macroscopic equations.

Models based on the FDT can capture important properties of complex fluids
even in non-equilibrium situations. For instance, the occurrence of $\zeta$ in
both terms for the noise and the deterministic motion (eq. (4)) mirrors the
FDT in Rouse bead-spring chains. Current research effortsIanniruberto,
Brasiello, and Marrucci (2019); Ianniruberto and Marrucci (2020) rely on the
FDT to replicate strongly nonlinear properties of polymer melts. Recently, the
FTD has been extended to active matter.Burkholder and Brady (2019)

A rigorous derivation of the FDT is restricted to thermodynamic equilibrium,
limiting its direct applicability to scenarios involving only linear
perturbations. This is also known as the FDT of the first kind. Conversely,
the FDT of the second kind extends to encompass more complex dynamics, such as
those observed in e.g. strong shear flows, where traditional linear approaches
fall short. It consists in linking the dissipative structure of the
macroscopic evolution to the distribution of their fluctuations.Montefusco,
Peletier, and Öttinger (2021) The credibility of this extended application is
supported by empirical successes and theoretical justifications based on
fundamental principles.Öttinger (2005) Yet its validity remains an open debate
(e.g. Matsumiya  _et al._ (2023)). Kim  _et al._ (2017) approximated reaction-
diffusion systems with the Landau-Lifshitz approach (eq. (1)) and found
inconsistencies in the results for small ensembles of molecules. A successive
work by Montefusco et al.Montefusco, Peletier, and Öttinger (2021) explained
the latter problems by revealing how not all processes can be consistently
described in terms of stochastic differential equations, since non-Gaussian
noise requires a different description.

In this paper we show that also the TNM requires a description beyond of non-
Gaussian noise. To facilitate this investigation, we introduce the framework
of GENERIC (General Equation for Non-Equilibrium Reversible-Irreversible
Coupling) formalism.Grmela and Öttinger (1997); Öttinger and Grmela (1997);
Öttinger (2005) This approach uniquely combines reversible dynamics with an
irreversible dissipative component, providing a cohesive structure for
articulating the FDT in a very natural manner. GENERIC expresses the evolution
of independent variables $\boldsymbol{x}$ in an nonequilibrium system isolated
from its surroundings, through the total energy $E$, the total entropy $S$ and
two linear operators $\boldsymbol{L},\,\boldsymbol{M}$, which capture the
physical constraints

$\frac{\mathop{}\\!\mathrm{d}\boldsymbol{x}}{\mathop{}\\!\mathrm{d}t}=\boldsymbol{L}(\boldsymbol{x})\cdot\partial_{\boldsymbol{x}}E(\boldsymbol{x})+\boldsymbol{M}(\boldsymbol{x})\cdot\partial_{\boldsymbol{x}}S(\boldsymbol{x}).$
(10)

As our focus lies on the dissipative characteristics, we omit the reversible
contribution ($L=0$) for the remainder of this work. The fluctuations link to
the irreversible component of eq. (10) through the FDT in its widespread
Green-Kubo formulationKubo, Toda, and Hashitsume (2023); Öttinger (2005)

$\boldsymbol{B}(\boldsymbol{x})\cdot\boldsymbol{B}(\boldsymbol{x})=2k_{\text{B}}\boldsymbol{M}(\boldsymbol{x}),$
(11)

where $\boldsymbol{B}(\boldsymbol{x})$ describes the Gaussian fluctuations
with the stochastic differential equation for the stochastic process
$\boldsymbol{X}_{\tau}$ given by:

$\mathop{}\\!\mathrm{d}\boldsymbol{X}_{t}=\bigr{[}\boldsymbol{M}(\boldsymbol{x})\cdot\partial_{\boldsymbol{x}}S(\boldsymbol{x})+k_{\text{B}}\partial_{\boldsymbol{x}}\cdot\boldsymbol{M}(\boldsymbol{x})\bigr{]}_{\boldsymbol{x}=\boldsymbol{X}_{t}}\mathop{}\\!\mathrm{d}t+\boldsymbol{B}(\boldsymbol{X}_{t})\cdot\mathop{}\\!\mathrm{d}\boldsymbol{W}_{t}.$
(12)

The necessity for this formulation arises from its foundational role in
modeling fluctuations within complex fluids. Of particular relevance to our
paper is the elucidation of noise enhancement dynamics. Specifically, for
Gaussian noise, one can readily isolate the dissipative component of the
macroscopic evolution equation given by the friction matrix $\boldsymbol{M}$,
as depicted in eq. 10. Subsequently, the noise can be readily enhanced on the
system with the formulation eq. (12). Conversely, given the noise component
($\boldsymbol{B}$), one can infer the corresponding macroscopic equation
$\boldsymbol{M}$. To date, this relationship has been perceived as a
comprehensive means of encapsulating noise in rheology.Hütter, Olmsted, and
Read (2020); Öttinger (2005)

A second formulation withing GENERIC uses the dissipation potential
$\Psi^{*}(\boldsymbol{x},\boldsymbol{\xi})$, instead of $\boldsymbol{M}$

$\frac{\mathop{}\\!\mathrm{d}\boldsymbol{x}}{\mathop{}\\!\mathrm{d}t}=\partial_{\boldsymbol{\xi}}\psi^{*}(\boldsymbol{x},\boldsymbol{\xi})\big{|}_{\boldsymbol{\xi}=\partial_{\boldsymbol{x}}S(\boldsymbol{x})}.$
(13)

Among the differences between eq. (10) and (13),Hütter and Svendsen (2013) we
highlight that eq. (10)-(11) is limited to to the description of Gaussian
noise, while $\psi^{*}$ can be used to model general Markow processes
according to a generalized FDTÖttinger, Peletier, and Montefusco (2021);
Montefusco, Peletier, and Öttinger (2021); Montefusco (2019)

$G_{\text{C}}=\frac{2k_{\text{B}}}{\tau}\ln{\langle
e^{\alpha(\boldsymbol{X}_{\tau}-\boldsymbol{x})}\rangle}\approx\psi^{*}(\boldsymbol{x},\boldsymbol{\xi})-\psi^{*}(\boldsymbol{x},\partial_{\boldsymbol{x}}S(\boldsymbol{x}))$
(14)

with $G_{\text{C}}$ as the cumulant generating function of the stochastic
process $\boldsymbol{X}_{\tau}$ in a time interval $\tau$. In presence of
fluctuations for $\alpha$ we have $\xi-\frac{\delta S}{\delta
p}=2k_{\text{B}}\alpha$. Consequently, we can think of $\xi$ as representing
the distance of the system form equilibrium. When $\xi=0$, the potential
$\psi^{*}$ equals its minimum at zero i.e. there is no dissipation in the
evolution of the macroscopic system. Eq. (14) captures how with growing
distance from equilibrium (larger values of $\xi$) the dissipative part of the
macroscopic evolution equation changes.

##  III Dissipation potentials

###  III.1 Origin of $\boldsymbol{c}$ fluctuations in the two models

To understand the difference in the fluctuations of the two models we first
focus on single trajectories of $\boldsymbol{Q}$’s components, for example in
flow direction ($Q_{1}$) (eq. (4) and (5-6)). Fig. 6 shows the different
origin of fluctuations in $\boldsymbol{c}$, namely the trajectories of
$\boldsymbol{Q}$. The latter and the following figure results from our
simulations. We solve the linear stochastich differential equation (4) of the
HDM with the Euler-Maruyama scheme.Öttinger (1996) For the TNM we apply the
algorithm presented by Biller and Petruccione (1990) to solve (5)-(6).
Basically the stochastic element in the trajectory is introduced through a
Gillespie algorithmGillespie (1977) for the strand’s lifetime. The simulations
are validated by comparison to analytical solution for
$\boldsymbol{c},p(\boldsymbol{Q})$.

Figure 6: Comparison of the individual trajectories in a simple shear flow of
$\boldsymbol{Q}$’s component in flow direction ($Q_{1}$). For the dumbbell,
the frequent, small collisions of the beads with the surrounding microscopic
particles enable the use of continuous sample paths with superimposed Gaussian
white noise $\mathop{}\\!\mathrm{d}\boldsymbol{W}_{t}$ (eq. (4)). In contrast,
the temporary network model exhibits deterministic trajectories, abruptly
interrupted by stochastic jumps when strands are replaced.

In the HDM noise is an intrinsic part of the dynamics (fig. 6). The frequent,
small interactions of microscopic particles allow the usage of continuous
sample paths with superposed Gaussian white noise
$\mathop{}\\!\mathrm{d}\boldsymbol{W}_{t}$. In contrast, the TNM exhibits
deterministic trajectories abruptly interrupted by stochastic jumps when
strands are replaced. Noise becomes apparent only at the coarse-grained level
of a discretized PDF (8), where the trajectories of a finite number of strands
are aggregated into a histogram. If the space of possible $\boldsymbol{Q}$
strands is partitioned in cells $j$ with volume $V_{j}$, the relative
frequency of ensemble members in the volume $V_{j}$ defines the probability
$P_{j}$. The piecewise constant function with value $p_{j}=P_{j}/V_{j}$ on the
cell $j$ provides a noisy histogram approximation to the probability density
$p(\bm{Q})$ (see fig. 7). Thus, the strand’s “jumps” cause fluctuations in the
columns of the histogram.

Figure 7: Typical histogram of the PDF for the temporary network model during
a step strain experiment at $t=\lambda\ln{2}$. Strands “jump” from one -column
$V_{j}$ to another when they are replaced. This stochastic changes in the
histogram cause subsequent fluctuations in the materials fluxes.

###  III.2 The dumbbells

Since the noise of the dumbbells is uncorrelated, the FDT results in the well
known Green-Kubo relation (11), which has been exploited to formulate noise in
$\boldsymbol{c}$ in past works.Hütter, Olmsted, and Read (2020); Öttinger
(2005) To express the PDF in terms of GENERIC and thus of $p(\boldsymbol{Q})$
functionals we haveÖttinger and Grmela (1997); Öttinger (2005)

$\begin{split}&S[p]=\int-p\frac{H\boldsymbol{Q}^{2}}{2T}-pk_{\text{B}}\ln{p}\mathop{}\\!\mathrm{d}Q^{3},\\\
&M[p]=-\partial_{\boldsymbol{Q}}\cdot\frac{2T}{\zeta}p(\boldsymbol{Q})\partial_{\boldsymbol{Q}}[\cdot].\\\
\end{split}$ (15)

It is then possible to construct a quadratic potential $\psi^{*}$

$\psi^{*}[p,\xi]=\frac{1}{2}M(\boldsymbol{Q})\xi^{2}=-\frac{1}{2}\partial_{\boldsymbol{Q}}\cdot\frac{2T}{\zeta}p(\boldsymbol{Q})\partial_{\boldsymbol{Q}}[\xi^{2}],$
(16)

so that both eq. (10) and eq. (13) reproduce the dissipative part of eq. (8):

$\begin{split}&\frac{\mathop{}\\!\mathrm{d}p}{\mathop{}\\!\mathrm{d}t}=\partial_{\boldsymbol{\xi}}\psi^{*}(\boldsymbol{x},\boldsymbol{\xi})\big{|}_{\boldsymbol{\xi}=\partial_{\boldsymbol{x}}S(\boldsymbol{x})}\\\
&=\partial_{\boldsymbol{Q}}\cdot\Big{[}\frac{2T}{\zeta}p\partial_{\boldsymbol{Q}}[\frac{H\boldsymbol{Q}^{2}}{2T}+k_{\text{B}}\ln{p}]\Big{]}\\\
&=\partial_{\boldsymbol{Q}}\cdot(\frac{2H}{\zeta}\boldsymbol{Q}p)+\frac{2k_{\text{B}}T}{\zeta}\partial_{\boldsymbol{Q}}\cdot\partial_{\boldsymbol{Q}}p.\end{split}$
(17)

Now, that eq. (16) allows us to express the irreversible part of the PDF (17),
we can derive $\psi^{*}$ for the TNM to compare the two dissipation
potentials.

###  III.3 The temporary network

One may verify that an entropy

$S[p]=-k_{\text{B}}\int
p(\boldsymbol{Q})\ln{\frac{p(\boldsymbol{Q})}{p_{\text{eq}}(\boldsymbol{Q})}}\mathop{}\\!\mathrm{d}^{3}Q,$
(18)

and a potential

$\psi^{*}[p,\xi]=\frac{4k_{\text{B}}}{\lambda}\int\sqrt{p(\boldsymbol{Q})p_{\text{eq}}(\boldsymbol{Q})}\big{(}\cosh{\frac{\xi(\boldsymbol{Q})}{2k_{\text{B}}}-1}\big{)}\mathop{}\\!\mathrm{d}^{3}Q$
(19)

reproduce the right-hand side of (8).

Now, armed with an enhanced understanding of fluctuations, as detailed in sec.
III.1 and illustrated in fig. 7, we proceed to apply the FDT as formulated in
eq. (14). This application involves a two-step limiting process - initially
expanding the ensemble sizes towards infinity, followed by the reduction of
the volumes $V_{j}$ of the cells of the observation cells to infinitesimally
small dimensions. Through this procedure, the resulting noisy histograms
gradually approximate the probability distribution $p(\boldsymbol{Q})$. The
theory of large deviations provides a quantitative framework for describing
this convergence,Montefusco (2019) offering precise insights into the behavior
of these distributions. In the realm of mathematics, this concept finds a more
rigorous formulation through the use of empirical measures. Here, we use the
physical insight of the transition from discrete $P_{j}$ to
$p(\boldsymbol{Q},t)$ to write

$\int\alpha_{N}(\bm{Q})\Delta
p(\bm{Q})d^{3}Q=\sum_{j}\alpha_{N}(\bm{Q}_{j})\Delta P_{j},$ (20)

where $\boldsymbol{Q}_{j}\in V_{j}$ and $\alpha_{N}(\bm{Q})$ indicates an
auxiliary quantity dependent on the $N$ strands constituting the ensemble. The
symbol $\Delta$ indicates the stochastic change in the time lapse $\tau$. For
small $\tau/\lambda$ we approximate the process by a discrete-time Markov jump
process. The transitions are the results of two independent Poisson processes
$L^{+}_{j},\,L^{-}_{j}$, where we justify the independence with two
assertions. First, $L^{+}_{j}$ depends on the entire ensemble, while
$L^{-}_{j}$ only depends on the members in cell $j$. Secondly, the constraint
of a constant ensemble size has no significant effect on the independence
property. A change in the ensemble composition can either take place, because
strands “die” with average $p_{-}=\frac{\tau}{\lambda}V_{j}p_{j}N$, or because
they form with average $p_{+}=\frac{\tau}{\lambda}V_{j}p_{\text{eq}}N$.
Given the non-Gaussian nature of the noise, we cannot use relation (11).
However, with the discrete representation of the stochastic process for small
ensembles, we express the left hand side of the generalized FDT (14)

$\begin{split}&G_{C}=\frac{2k_{\text{B}}}{\tau}\ln{\big{\langle}\exp{\sum_{j}\frac{\alpha_{N}(\boldsymbol{Q}_{j})}{N}(L_{j}^{+}-L_{j}^{-})}\big{\rangle}}\\\
&\overset{\text{ind.}}{=}\frac{2k_{\text{B}}}{\tau}\sum_{j}\ln{\big{\langle}\exp{\frac{\alpha_{N}(\boldsymbol{Q}_{j})}{N}L_{j}^{+}}\big{\rangle}}+\ln{\big{\langle}\exp{-\frac{\alpha_{N}(\boldsymbol{Q}_{j})}{N}L_{j}^{-}}\big{\rangle}}\\\
&\overset{\text{Poi.}}{=}\frac{2k_{\text{B}}N}{\lambda}\sum_{j}V_{j}\big{[}p_{\text{eq}}(\boldsymbol{Q}_{j})(e^{\alpha_{N}(\boldsymbol{Q}_{j})/N}-1)+p(\boldsymbol{Q}_{j})e^{\alpha_{N}(\boldsymbol{Q}_{j})/N}-1)\big{]}\\\
&=\frac{2k_{\text{B}}N}{\lambda}\Big{(}\int
p_{\text{eq}}(\boldsymbol{Q})e^{\alpha_{N}(\boldsymbol{Q})/N}+p(\boldsymbol{Q})e^{-\alpha_{N}(\boldsymbol{Q})/N}\mathop{}\\!\mathrm{d}^{3}Q-2\Big{)}.\end{split}$
(21)

In the second line we used the independence of the two processes, in the third
the cumulant generating function of a Poisson random variable, in the fourth
we assumed the limit of large ensembles and vanishing $V_{j}$. $G_{C}$ is
independent of the time interval $\tau$ and of the particular discretization
of the PDF. However, it depends on the ensemble size $N$. Large deviation
theory considers the effect of a single jump event, so that we are interested
in $N=1$ and write $\alpha_{1}(\bm{Q})=\alpha(\bm{Q})$. Then, the general FDT
(14) holds as $\psi^{*}$ from eq. (8) agrees with the microscopic fluctuations
in eq. (21).
This result reflects the distinct jump dynamics for one single strand depicted
in fig. 6. The implication for the fluctuations in $\boldsymbol{c}$ and thus
$\boldsymbol{\tau}$ are further explored in sec. V. The potential (19)
resembles more the one of chemical reactionsMontefusco, Peletier, and Öttinger
(2021) than the one of diffusion processes (16). This is no surprise as the
dynamic of the processes is close to the “jumps” of a chemical reaction. An
expansion of the $\cosh{}$ function shows that at higher values of the non-
equilibrium parameter $\xi$ the potential (19) differs more and more from the
quadratic potential of diffusive process as of the HDM.

##  IV Differences in the PDFs

At equilibrium both HDM and TNM exhibit the same distribution $p_{\text{eq}}$.
Under flow the stress $\boldsymbol{\tau}$ in both models is described by the
UCM (which also holds for the evolution of $\boldsymbol{c}$. However, sec. III
implies different fluctuations in $\boldsymbol{c}$ and thus underlying
different PDFs. To characterize this we calculate the evolution equation of
the 4th order moment, which for the Gaussian HDM model evolves according to
Wick’s theorem :

$\displaystyle\begin{split}\partial_{t}\mu_{abcd}&=\frac{\partial
u_{a}}{\partial x_{j}}\mu_{jbcd}+\mu_{ajcd}\frac{\partial u_{b}}{\partial
x_{j}}+\mu_{abjd}\frac{\partial u_{c}}{\partial
x_{j}}+\mu_{abcj}\frac{\partial u_{d}}{\partial x_{j}}\\\
&-\frac{8H}{\zeta}\mu_{abcd}+\frac{4k_{\text{B}}T}{\zeta}(\delta_{ab}\mu_{cd}+\delta_{bc}\mu_{ad}+\mu_{ab}\\\
&\delta_{cd}+\delta_{ac}\mu_{bd}+\delta_{bd}\mu_{ac}+\delta_{ad}\mu_{bc}),\end{split}$
(22)

For the non-Gaussian TNM PDF the 4th order moment evolves accoridng to:

$\displaystyle\begin{split}\partial_{t}\mu_{abcd}&=\frac{\partial
u_{a}}{\partial x_{j}}\mu_{jbcd}+\mu_{ajcd}\frac{\partial u_{b}}{\partial
x_{j}}+\mu_{abjd}\frac{\partial u_{c}}{\partial
x_{j}}+\mu_{abcj}\frac{\partial u_{d}}{\partial x_{j}}\\\
&+\lambda^{-1}(\mu_{abcd}^{\text{eq}}-\mu_{abcd}),\end{split}$ (23)

where $\mu_{abcd}^{\text{eq}}$ is the equilibrium 4th moment.

Analytical solutions to the step strain experiment can be used to illustrate
the differences between the evolution of the PDF’s. Before the instantaneous
deformation is imposed both model predict $p_{\text{eq}}$ at equilibrium. The
shear rate equals $\dot{\gamma}=\gamma\delta(t)$, so that integration of (4)
and (5)-(6) from $0^{-}$ to $0^{+}$ leads to
$\boldsymbol{Q}_{0}=\Gamma\cdot\boldsymbol{Q}_{\text{eq}}$ with
$\boldsymbol{\Gamma}=\exp{(\nabla\boldsymbol{u}^{T})}$, where the subscript eq
refers to equilibrium opposed to $0$ as the initial condition. Since the
configurations are convected, the PDFs remain Gaussian with first and second
order moment $\langle Q_{1}Q_{1}\rangle$ is proportional to the normal stress
$\tau_{11}$ according to

$\displaystyle\langle\boldsymbol{Q}\rangle_{0}=\boldsymbol{\Gamma}\cdot\langle\boldsymbol{Q}\rangle_{\text{eq}},$
(24)
$\displaystyle\langle\boldsymbol{Q}\boldsymbol{Q}\rangle_{0}=\boldsymbol{\Gamma}\cdot\langle\boldsymbol{Q}\boldsymbol{Q}\rangle_{\text{eq}}\cdot\boldsymbol{\Gamma}^{T}.$
(25)

To simplify the graphical representation we restrict ourselves to the $Q_{1}$
component so to observe the marginal distribution function $p(Q_{1})=\int\int
p(\boldsymbol{Q})\mathop{}\\!\mathrm{d}Q_{2}\mathop{}\\!\mathrm{d}Q_{3}$. Then
eq. (8) results in the superposition of two contributions

$p_{\text{TNM}}(t)=p_{\text{eq}}(1-\exp{(-\lambda^{-1}t)})+p_{0}\exp{(-\lambda^{-1}t)}.$
(26)

On the other hand, the PDF evolution eq. (7) is a Fokker-Planck equation with
linear coefficients. It consequently leads to a Gaussian distribution uniquely
determined by its first and second moments with known respective evolution
equations.Venerus and Öttinger (2018) Fig. 8 compares the two different
analytical solutions. Only at equilibrium (i.d. at $t<0$ or
$t\rightarrow\infty$) and at the initial condition $t=0$ are the two PDFs
equal. At all intermediate times, $p_{\text{TNM}}$ is non-Gaussian. The
distinction is best captured by fig. 9. It depicts the fourth cumulant
$C_{1111}$, which in 1-D corresponds to $C_{1111}=\mu_{1111}-3\mu_{11}^{2}$
and thus reflects the differences in eq. (22) and (23).
For flows strengths beyond our analytical framework the discrepancy is visible
in numerical simulations. For instance, in fig. 10, the impact of different
shear rates ($\dot{\gamma}$) on a start-up flow is depicted in relation to the
differences $\mu_{1212}^{\text{TNM}}-\mu_{1212}^{\text{HDM}}$. Given that the
distribution for the HDM is always Gaussian, this difference equates
$C_{1212}^{\text{TNM}}$, since $C_{1212}^{\text{HDM}}$ is constant $0$. In
agreement with the findings of sec. III, the stronger the perturbation from
equilibrium, the larger the difference.

Figure 8: The probability density function of the temporary network (crosses)
and the dumbbells (line) for different times after applying an instantaneous
step strain at $t=0$. The strain is of $1.5$ instead of $1$ to accentuate the
difference. Figure 9: 4th cumulant $C_{1212}$ as for a step strain
($\gamma=1$), an ensemble of $N=10^{5}$ and a time discretization of
$\mathop{}\\!\mathrm{d}t=0.03\lambda$ for both the temporary network and the
dumbbells. The maximum in the difference occurs at $t=\lambda\ln{2}$ and
relaxes slowly - in fact only for $t\rightarrow\infty$ the difference ceases
to exist. Figure 10: 4th cumulant $C_{1212}$ of the TNM for start up flows at
two different shear rates, an ensemble of $N=10^{5}$ and a time discretization
of $\mathop{}\\!\mathrm{d}t=0.03\lambda$. In agreement with the findings of
sec. III, the further the systems are from equilibrium (larger values of
$\xi$, or specifically for this example $\dot{\gamma}$), the more the TNM and
the HDM differ from eachother.

Until this point the discussion was confined to the configurational space of
$\boldsymbol{Q}$, its single trajectories or PDFs. In the following section we
elucidate how these result relate to the stress fluctuations in actual
volumes.

##  V The limit of Gaussian noise for stress fluctuations in complex fluids

In fluctuating hydrodynamics the relation for the fluctuations of a field
variable $\delta f(\boldsymbol{x},\boldsymbol{t})$ and the volume $\Delta V$
in which the latter is measured is

$\delta f_{V}(t)=\frac{1}{\Delta V}\int_{\Delta V}\delta
f(\boldsymbol{x},t)\mathop{}\\!\mathrm{d}^{3}x.$ (27)

The noise, normalized by $\Delta V$, is Gaussian, independent of the volume of
observation.Zárate and Sengers (2006) This same relation is recovered here for
the HDM,Hütter, Olmsted, and Read (2020); Öttinger (2005) a cartoon
representing the physical picture is given in fig. 11.

Figure 11: The difference in $\langle\delta\boldsymbol{c}^{2}\rangle$ is
maximal on the scale of _one_ structural variable and becomes
indistinguishable with increasing number of dumbbells/strands $N_{Q}$ scaling
as in eq. (28). The cartoon depicts how $N_{Q}$ i.e. the volume $V=N_{Q}/n$
relates to different characteristics lengths for the two systems.

In contrast, sec. III points out that for ensembles confined to $N=1$ the
TNM’s fluctuations are Poissonian. On small length scales they cannot be
described by Gaussian noise.

How do these findings relate to stress measurements? The origin of stress
fluctuations in these complex materials is the evolution of $\boldsymbol{c}$
(see eq. (9)). Thus, in homogeneous flow fluctuations
$\delta{\boldsymbol{c}}=\langle\boldsymbol{c}\rangle-\boldsymbol{c}$ directly
reflect the different statics of the HDM and the TNM. When the control volume
increases (right-hand side of fig. 11)), $\boldsymbol{c}$ results from the sum
of independent $\boldsymbol{Q}$’s contributions. Consequently, according to
the central limit theorem, larger volumes lead to indistinguishable Gaussian
noise. To quantify the convergence to the same noise, we define $N_{Q}$ as the
amount of $\boldsymbol{Q}_{k}$ strands/dumbbells in $\Delta V$, where $N_{Q}$
linearly relates to $\Delta V$ through the number density $n$. A local
measurement for a larger $\Delta V$, as in the right-hand picture of fig. 11,
consequently results in
$\boldsymbol{c}=\langle\boldsymbol{Q}_{k}^{2}\rangle_{N_{Q}}$. The averages of
$\boldsymbol{c}$ over $N_{c}$ measurements yield

$\begin{split}\langle\boldsymbol{c}^{2}\rangle_{N_{c}}&=\langle\langle\boldsymbol{Q}\boldsymbol{Q}\rangle_{N_{Q}}^{2}\rangle_{N_{c}}=\frac{1}{N_{c}N_{Q}^{2}}\sum_{j=1}^{N_{c}}\big{(}\sum_{k=1}^{N_{Q}}\boldsymbol{Q}_{k}\boldsymbol{Q}_{k}\big{)}^{2}\\\
&=\frac{1}{N_{Q}}\langle\boldsymbol{Q}\boldsymbol{Q}\boldsymbol{Q}\boldsymbol{Q}\rangle+\frac{N_{Q}-1}{N_{Q}}\langle\boldsymbol{Q}_{m}^{2}\boldsymbol{Q}_{l}^{2}\rangle\\\
&=\frac{1}{N_{Q}}\langle\boldsymbol{Q}\boldsymbol{Q}\boldsymbol{Q}\boldsymbol{Q}\rangle+\frac{N_{Q}-1}{N_{Q}}\boldsymbol{c}^{2},\end{split}$
(28)

where we applied binomial coefficients and the independence of
$\boldsymbol{Q}_{m}$ and $\boldsymbol{Q}_{l}$. Thence, the difference
$\langle\delta{\boldsymbol{c}}_{\text{TNM}}^{2}\rangle-\langle\delta{\boldsymbol{c}}_{\text{HDM}}^{2}\rangle$
scales with $1/N_{Q}$.

In agreement with the findings of sec. III, when one trajectory is considered
(left-hand side of fig. 11) the difference is maximal and corresponds to the
TNM’s fourth cumulant
$C_{ijkl}^{\text{TNM}}=\mu_{ijkl}^{\text{TNM}}-\mu_{ijkl}^{\text{HDM}}$. While
in the HDM the critical length scale associated to one trajectory is given by
the average distance $d$, in the case of the TNM we have $d=|\boldsymbol{Q}|$.

The measurement of the different $\delta\boldsymbol{c}$ in the two models,
requires a certain amount of data points. Numerical simulations of the two
models quantify exactly how many. The error in the 4th moment $\mu_{ijkl}$ is

$\epsilon_{ijklijkl}=\sqrt{\frac{\mu_{ijklijkl}-\mu_{ijkl}^{2}}{N_{c}-1}},$
(29)

so that for a detectable difference we need

${\epsilon}^{\text{HDM}}_{ijklijkl}+{\epsilon}^{\text{TNM}}_{ijklijkl}<<C_{ijkl}^{\text{TNM}}$
(30)

to hold. Interestingly, an application of these formulas to the analytical
solution for the step-strain example (8) shows $N_{c}$ to be independent from
the Maxwell parameters $G,\lambda$, but strongly dependent on $\gamma$ with a
leading order of $\mathcal{O}(\gamma^{-4})$. Differences in
$\delta\boldsymbol{c}$ are most easily to observe far from equilibrium
conditions as suggested by the potential in sec. III and fig. 8.

##  VI Discussion

###  VI.1 Relevance of stress fluctuations

Fluctuations are pivotal for _modeling coupled phenomena_ that a purely
continuum mechanical description cannot intuitively capture. Examples include
reaction-diffusion systems, where fluctuations explain phenomena such as
pattern formations,Wang  _et al._ (2007) fluctuation-induced
instabilities,Kessler and Levine (1998) front propagation dynamics,Khain, Lin,
and Sander (2011) and the emergence of new steady states.Togashi and Kaneko
(2004) Fluctuating hydrodynamics (1) with Gaussian noise proved able to
capture many of these effects in Newtonian fluids.Kim  _et al._ (2017) For
networks at specific length scales, however, we need the noise description
introduced in the present work. The noise enhancement of a UCM of Hütter,
Olmsted, and Read (2020) is only a _particular_ case reflecting the dynamics
of dumbbells.

Likewise, the results are relevant for _interpreting experiments_ involving
diffusion coupled with non-zero stresses. An illustrative example of this
application is observed in the TNM behavior of stretched fire ant
networks.Wagner  _et al._ (2023) In this scenario, the structural variable
$\boldsymbol{Q}$ extends as long as the ants’ legs. Hence, diffusive phenomena
in these stretched networks shall reflect the Poissonian fluctuations
dynamics. Likewise, microrheology measurements of a bead embedded in a network
under deformation Goren  _et al._ (2024); Hasnain and Donald (2006)
necessitate a deeper understanding of how Brownian forces change relative to
thermodynamic equilibrium. We demonstrated that even in a simple network model
like the TNM, the presence of stresses alters the distribution of stochastic
forces (see sec. VI.2 and fig. 12).

This study also underscores the importance of _stress fluctuations (
$\delta\boldsymbol{c},\,\delta\boldsymbol{\tau}$) as indicators of the
underlying structure_. Theoretical studies and simulations looking into
microscopic mechanisms governing the rheological behavior of colloidal
networks by analyzing spontaneous microscopic dynamics, confirming a relation
between microstructural dynamics and microscopic structure.Bantawa  _et al._
(2022) Similarly, stress correlations
$\langle\delta\boldsymbol{\tau}^{2}\rangle$, albeit in spatial dimensions
rather than temporal, were used to identify force chains in colloidal
networks, underscoring the significance of this metric in characterizing
microstructures.Vinutha  _et al._ (2023) Recent experiments and simulations
have also observed non-Gaussian fluctuations near the glass transition of
glass-forming liquids.Rusciano, Pastore, and Greco (2022) These studies
highlight the importance of changes of the temporal displacement distributions
in networks with increasing distance from equilibrium. We therefore suggest
the possible use of deviations from Gaussian Statistics in stress fluctuations
as a means to discern local loads in viscoelastic networks. Such an approach
may enable the detection of local stress concentrations. The present paper
identifies the general FDT formulation (Eq. (14)) as a suitable framework to
consistently coarse-grain such noise, offering a detailed examination of
fluctuations in a mesoscopic network representation (TNM) for microscopic
reference.

###  VI.2 Experimental determination of the difference in
$\delta\boldsymbol{c}$

To detect the differences in $\delta\boldsymbol{c}$ and, thus, in the stress
$\boldsymbol{\tau}$ (eq. (9)), several experimental methods may be :
_microrheology, Raman scattering, polarized fluorescence, polarimetry_.
In _microrheology_ the mean square displacement of small colloids ($0.1$ to
$10$ µm) suspended in a fluid is measured via microscopy or light scattering
and related to the linear viscoelastic fluid properties.Mason and Weitz
(1995); Savin and Doyle (2007); van Zanten, Amin, and Abdala (2004); He  _et
al._ (2010); Furst and Squires (2017) However, the techniques relies on an
equilibrium formulation. Fig. 1 shows that experimental acquisition of the
trajectories under flow is now possible. Fig. 12 illustrates the origin of the
bead’s fluctuations - namely the colloid’s interaction with the fluids
structure, as captured by the structural variable $\boldsymbol{Q}$. The bead’s
trajectories hence should reflect the distinct nature of fluctuations in the
TNM and the HDM in flow.

Figure 12: The fluctuations of a hard sphere immersed in a fluid during flow
reflect the fluctuations of the fluid’s structure through its structural
variable $\boldsymbol{Q}$.

_Raman scattering_ defines the frequency shift in a small portion of scattered
light due to vibrational and rotational motions within a sample.Fuller (1995)
The Raman tensor $\alpha_{ij}$ establishes a connection between the electric
vector of incident light and Raman scattered light. Since $\alpha_{ij}$ is a
function of segment vectors $\boldsymbol{Q}$, it enables both 2nd and 4th
moments to be determined (see, e.g., eq. (6)-(8) in Archer, Fuller, and
Nunnelley (1992)). The method has been successfully employed for the
quantitative description of $\langle\boldsymbol{Q}\boldsymbol{Q}\rangle$ and
$\langle\boldsymbol{Q}\boldsymbol{Q}\boldsymbol{Q}\boldsymbol{Q}\rangle$ for
LC systems,Huang, Archer, and Fuller (1996a) for elongated or sheared polymer
melts,Bower (1972); Purvis and Bower (1976); Archer, Huang, and Fuller (1994);
Huang, Archer, and Fuller (1996b) and for _in situ_ observation of
polymerization in a rheometer.Nadège  _et al._ (2013) In Archer, Huang, and
Fuller (1994), anisotropies were measured for an entangled polyisobutylene
melt during step uniaxial extension and compared with predictions from
different models, including the TNM. Consequently, Raman scattering has
already proven to be a valuable tool to confirm the difference in
$\delta\boldsymbol{c}$, but possibly time resolution can be improved.
Similarly, _polarized fluorescence_ provides the 2nd and 4th moments. It
relies on the different wavelengths of emitted and incident lightFuller
(1995); Larson (1999) and has been successfully applied to polymer
dynamics.Monnerie (1987); Lapersonne  _et al._ (1989); Viovy and Monnerie
(1985)
Linear birefringence using _polarimetry_ uses the dependency of refraction
from segmental orientation. The polarization by an electric field (such as a
light beam) results in a refractive index tensor $\boldsymbol{n}$. The latter
is then linearly related to $\boldsymbol{c}$ through the stress-optical
law.Larson (1999); Janeschitz-Kriegl (1983) Full tensor optical rheometry was
able to determine $\boldsymbol{c}$ in polymeric flowsKalogrianitis and van
Egmond (1997) and may also be extended to detect differences in
$\delta\boldsymbol{c}$ (28) through point-wise measurements (i.e. using cross
correlation techniques for beams sent at slightly different angles).

##  VII Summary

Fluctuations play a crucial role in understanding various processes, spanning
from phenomena such as diffusion to slip, and they must be accounted for in
order to accurately model coupled phenomena. Furthermore, even though often
perceived merely as noise, fluctuations directly probe the material structure
and in a unique manner, thereby uncovering the coupling of local dynamics and
structure. However, our understanding of fluctuations is largely based on
fluctuating hydrodynamics, which describes Newtonian fluids assuming Gaussian
noise.

The central message of this paper is that the Gaussian approximation does not
always hold in complex fluids. Gaussian noise, as assumed in fluctuating
hydrodynamics or in the widely used Fluctuation-Dissipation Theorem
formulation of Green-Kubo (eq. (11)), cannot describe the noise at the length
scale of $\boldsymbol{Q}$ in a temporary network model. These findings serve
as a caveat, that fluctuations in complex flows may differ fundamentally from
those in Newtonian fluids. We introduce the general formulation of the
Fluctuation-Dissipation Theorem (14) as an effective solution to this problem.
This formulation can account for the discrete jump dynamics occurring in
networks, making it a potent tool for developing a fluctuating rheology akin
to fluctuating hydrodynamics.

Second, Whereas the stress evolutions for dumbbell model and the temporary
network model are identical, the stress fluctuations do pick up the
differences in microstructure flow coupling. As measuring techniques advance,
there is a shift towards less coarse-grained descriptions in rheology. We
demonstrate how, within the context of local field theory, the most
fundamental quasi-linear viscoelastic model (UCM) exhibits varying
fluctuations depending on the underlying microscopic dynamics.

Finally, we provide quantification of these differences in stress
fluctuations. We suggest using the 4th moment to capture differences in
$\delta\boldsymbol{c}$, introduce an analytical solution to a step strain
example, and discuss the convergence behavior of $\delta\boldsymbol{c}$ with
increasing control volumes. Finally, we emphasize the relevance of these
findings for current research and list techniques to detect these differences.

###### Acknowledgements.

Thanks to Alberto Montefusco for the discussions on his seminal
work.Montefusco (2019) Thanks to Vincenzo Ianniello, Florence Müller and
Martin Kröger (ETH Zürich) for precious feedback about the work. At last,
thanks to Pierre Lehéricey (ETH Zürich) for helping with the experiments to
produce fig. 1.

## References

## References

  * Brown (1828) R. Brown, “A brief account of microscopical observations on the particles contained in the pollen of plants,” The Philosophical Magazine , 161–173 (1828).
  * Einstein (1905) A. Einstein, “Eine neue Bestimmung der Moleküldimensionen,” Ann. d. Physik 17, 549 (1905).
  * Zárate and Sengers (2006) J. Zárate and J. Sengers,  _Hydrodynamic Fluctuations in Fluids and Fluid Mixtures_ (Elsevier, 2006).
  * Larsen and Furst (2008) T. H. Larsen and E. M. Furst, “Microrheology of the liquid-solid transition during gelation,” Phys. Rev. Lett. 100, 146001 (2008).
  * Venerus and Öttinger (2018) D. Venerus and H. C. Öttinger,  _A Modern Course in Transport Phenomena_ (Cambridge University Press, 2018).
  * Mason and Weitz (1995) T. Mason and D. Weitz, “Optical measurements of frequency-dependent linear viscoelastic moduli of complex fluids,” Phys. Rev. Lett. 74, 1250 (1995).
  * Furst and Squires (2017) E. M. Furst and T. M. Squires,  _Microrheology_ (Oxford University Press, 2017).
  * Öttinger (1996) H. C. Öttinger,  _Stochastic Processes in Polymeric Fluids_ (Springer, Berlin, 1996).
  * Marbach, Dean, and Bocquet (2018) S. Marbach, D. Dean, and L. Bocquet, “Transport and dispersion across wiggling nanopores,” Nature Physics 14, 1108–1113 (2018).
  * Marbach, Zheng, and Holmes-Cerfon (2022) S. Marbach, J. Zheng, and M. Holmes-Cerfon, “The nanocaterpillar’s random walk: Diffusion with ligand-receptor contacts,” Soft Matter 18 (2022).
  * Landau and E.M. (1959) L. Landau and L. E.M.,  _Statistical Physics_ , Vol. 5 (Pergamon Press, 1959).
  * Kirkpatrick, Ortiz de Zárate, and Sengers (2013) T. R. Kirkpatrick, J. M. Ortiz de Zárate, and J. V. Sengers, “Giant casimir effect in fluids in nonequilibrium steady states,” Phys. Rev. Lett. 110, 235902 (2013).
  * Colombo  _et al._ (2019) G. Colombo, R. Massaro, S. Coleman, J. Läuger, P. V. Puyvelde, and J. Vermant, “Ultrafast imaging of soft materials during shear flow,” Korea-Australia Rheology Journal 31, 229–240 (2019).
  * Snijkers  _et al._ (2009) F. Snijkers, G. D’Avino, P. L. Maffattone, F. Greco, M. Hulsen, and J. Vermant, “Rotation of a sphere in a viscoelastic liquid subjected to shear flow. part ii: Experimental results,” J. Rheol. 53, 459 (2009).
  * Ershov  _et al._ (2022) D. Ershov, M.-S. Phan, J. Pylvänäinen, S. Rigaud, L. Blanc, A. Charles-Orszag, J. Conway, R. Laine, N. Roy, D. Bonazzi, G. Dumenil, G. Jacquemet, and J.-Y. Tinevez, “Trackmate 7: integrating state-of-the-art segmentation algorithms into tracking pipelines,” Nature Methods 19 (2022).
  * Hütter, Olmsted, and Read (2020) M. Hütter, P. Olmsted, and D. Read, “Fluctuating viscoelasticity based on a finite number of dumbbells,” Eur. Phys. J. E 43 (2020).
  * Öttinger and Grmela (1997) H. C. Öttinger and M. Grmela, “Dynamics and thermodynamics of complex fluids. ii. illustrations of a general formalism,” Phys. Rev. E 56, 6633–6655 (1997).
  * Ferry (1980) J. D. Ferry,  _Viscoelastic Properties of Polymers_ , 3rd ed. (John Wiley & Sons, Inc., 1980).
  * Larson (1999) R. G. Larson,  _The Structure of Rheology of Complex Fluids_ (Oxford University Press, 1999).
  * Morrison (2001) F. A. Morrison,  _Understanding Rheology_ (Oxford University Press, New york, 2001).
  * Öttinger, Peletier, and Montefusco (2021) H. C. Öttinger, M. A. Peletier, and A. Montefusco, “A Framework of Nonequilibrium Statistical Mechanics. I. Role and Types of Fluactuations,” J. Non-Equilib. Thermodyn. 46, 1–13 (2021).
  * Costanzo  _et al._ (2022) S. Costanzo, V. Ianniello, R. Pasquino, N. Grizzuti, G. Ianniruberto, and G. Marrucci, “Strain hardening of unentangled polystyrene solutions in fast shear flows,” Macromolecules 55, 9206–9219 (2022).
  * Annable, Buscall, and Ettelaie (1996) T. Annable, R. Buscall, and R. Ettelaie, “Network formation and its consequences for the physical behaviour of associating polymers in solution,” Colloids and Surfaces A: Physicochemical and Engineering Aspects 112, 97–116 (1996), a collection of Papers Presented at the Associating Polymers 1995 Meeting.
  * Maxwell (1867) Maxwell, “On the dynamicsl theory of gases,” Proc. Roy. Soc. A 15, 49 (1867).
  * Oldroyd (1950) J. G. Oldroyd, “On the formulation of rheological equations of state,” Proc. Roy. Soc. A 200, 523 (1950).
  * Lodge (1974) A. Lodge,  _Body tensor fields in continuum mechanics_ , 1st ed. (Academic Press Inc., New York, 1974).
  * Macosko (1994) C. W. Macosko,  _RHEOLOGY Principles, Measurements, and Applications_ , 1st ed. (Wiley-VCH, New York, 1994).
  * Phan-Thien (2002) N. Phan-Thien,  _Understanding Viscoelasticity: Basic of Rheology_ , 1st ed. (Springer, Heidelberg, 2002).
  * Rouse (1953) E. Rouse, J. Chem. Phys. 21, 1272 (1953).
  * Green and Tobolsky (2004) M. S. Green and A. V. Tobolsky, “A New Approach to the Theory of Relaxing Polymeric Media,” J. Chem. Phys. 14, 80–92 (2004).
  * Biller and Petruccione (1990) P. Biller and F. Petruccione, “Continuous time simulation of transient polymer network models,” J. Chem. Phys. 92, 6322 (1990).
  * Bird, Hassager, and Curtiss (2017) R. Bird, R. Hassager, and C. Curtiss,  _Dynamics of Polymeric Liquids_ , 2nd ed. (Wiley, New York, 2017).
  * Nyquist (1928) H. Nyquist, “Thermal agitation of electric charge in conductors,” Phys. Rev. 32, 110–113 (1928).
  * Ianniruberto, Brasiello, and Marrucci (2019) G. Ianniruberto, A. Brasiello, and G. Marrucci, “Modeling unentangled polystyrene melts in fast elongational flows,” Macromolecules 52, 4610–4616 (2019).
  * Ianniruberto and Marrucci (2020) G. Ianniruberto and G. Marrucci, “Origin of shear thinning in unentangled polystyrene melts,” Macromolecules 53 (2020).
  * Burkholder and Brady (2019) E. W. Burkholder and J. F. Brady, “Fluctuation-dissipation in active matter,” J. Chem. Phys. 150, 184901 (2019).
  * Montefusco, Peletier, and Öttinger (2021) A. Montefusco, M. A. Peletier, and H. C. Öttinger, “A Framework of Nonequilibrium Statistical Mechanics. II. Coarse-Graining,” J. Non-Equilib. Thermodyn. 46, 15–33 (2021).
  * Öttinger (2005) H. C. Öttinger,  _Beyond Equilibrium Thermodynamics_ (William and Sons, Inc., New Jersey, 2005).
  * Matsumiya  _et al._ (2023) Y. Matsumiya, T. Sato, Q. Chen, and H. Watanabe, “Rouse analysis of nonlinear rheology of unentangled polymer melts under fast shear: Viscoelastic response to superposed oscillatory strain,” Macromolecules 56 (2023).
  * Kim  _et al._ (2017) C. Kim, A. Nonaka, J. B. Bell, A. L. Garcia, and A. Donev, “Stochastic simulation of reaction-diffusion systems: A fluctuating-hydrodynamics approach,” J. Chem. Phys. 146, 124110 (2017).
  * Grmela and Öttinger (1997) M. Grmela and H. C. Öttinger, “Dynamics and thermodynamics of complex fluids. i. development of a general formalism,” Phys. Rev. E 56, 6620–6632 (1997).
  * Kubo, Toda, and Hashitsume (2023) R. Kubo, M. Toda, and N. Hashitsume,  _Statistical physics II. Nonequilibrium statistical mechanics. Rev. transl. of the rev. orig. Jap. ed. 2nd ed_ (Springer, 2023).
  * Hütter and Svendsen (2013) M. Hütter and B. Svendsen, “Quasi-linear versus potential-based formulations of force-flux relations and the GENERIC for irreversible processes: comparisons and examples,” Continuum Mech. Thermodyn. 25, 803–816 (2013).
  * Montefusco (2019) A. Montefusco,  _Dynamic Coarse-Graining via Large-Deviation Theory_ , Dissertation, ETHZ (2019).
  * Gillespie (1977) D. T. Gillespie, “Exact stochastic simulation of coupled chemical reactions,” J. Phys. Chem. 81, 2340–2361 (1977).
  * Wang  _et al._ (2007) H. Wang, Z. Fu, X. Xu, and Q. Ouyang, “Pattern formation induced by internal microscopic fluctuations,” J. Phys. Chem. A 111, 1265–1270 (2007).
  * Kessler and Levine (1998) D. A. Kessler and H. Levine, “Fluctuation-induced diffusive instabilities,” Nature 394, 556–558 (1998).
  * Khain, Lin, and Sander (2011) E. Khain, Y. Lin, and L. M. Sander, “Fluctuations and stability in front propagation,” Europhys. Lett. 93, 28001 (2011).
  * Togashi and Kaneko (2004) Y. Togashi and K. Kaneko, “Molecular discreteness in reaction-diffusion systems yields steady states not seen in the continuum limit,” Phys. Rev. E 70, 020901 (2004).
  * Wagner  _et al._ (2023) R. J. Wagner, S. Lamont, Z. T. White, and F. J. Vernerey, “The mechanical response of fire ant rafts,” arXiv: (2023).
  * Goren  _et al._ (2024) S. Goren, B. Ergaz, D. Barak, R. Sorkin, and A. Lesman, “Anisotropic micromechanics of fibrous gels stretched illuminate the cellular microenvironment,” Biophysical Journal 123, 557a (2024).
  * Hasnain and Donald (2006) I. Hasnain and A. M. Donald, “Microrheological characterization of anisotropic materials,” Physical Review 73, 031901 (2006).
  * Bantawa  _et al._ (2022) M. Bantawa, B. Keshavarz, M. Geri, M. Bouzid, T. Divoux, G. McKinley, and E. Del Gado, “The hidden hierarchical nature of soft particulate gels,” Nature Physics 8, 1178 (2022).
  * Vinutha  _et al._ (2023) H. Vinutha, F. D. Diaz Ruiz, X. Mao, B. Chakraborty, and E. Del Gado, “Stress–stress correlations reveal force chains in gels,” J. Chem. Phys. 158 (2023).
  * Rusciano, Pastore, and Greco (2022) F. Rusciano, R. Pastore, and F. Greco, “Fickian non-gaussian diffusion in glass-forming liquids,” Physical Review Letters 128, 168001 (2022).
  * Savin and Doyle (2007) T. Savin and P. S. Doyle, “Statistical and sampling issues when using multiple particle tracking,” Phys. Rev. E 76, 021501 (2007).
  * van Zanten, Amin, and Abdala (2004) J. van Zanten, S. Amin, and A. Abdala, “Brownian motion of colloidal spheres in aqueous peo solutions,” Macromolecules 37, 3874–3880 (2004).
  * He  _et al._ (2010) F. He, G. W. Becker, J. R. Litowski, L. O. Narhi, D. M. Brems, and V. I. Razinkov, “High-throughput dynamic light scattering method for measuring viscosity of concentrated protein solutions,” Analytical Biochemistry 399, 141–143 (2010).
  * Fuller (1995) G. G. Fuller,  _Optical Rheometry of Complex Fluids_ , 1st ed. (Oxford University Press, New York, 1995).
  * Archer, Fuller, and Nunnelley (1992) L. Archer, G. Fuller, and L. Nunnelley, “Dynamics of polymeric liquids using polarization-modulated laser raman scattering,” Polymer 33, 3574–3581 (1992).
  * Huang, Archer, and Fuller (1996a) K. Huang, L. A. Archer, and G. G. Fuller, “Polarization-modulated Raman scattering measurements of nematic liquid crystal orientation,” Review of Scientific Instruments 67, 3924–3930 (1996a).
  * Bower (1972) D. I. Bower, “Investigation of molecular orientation distributions by polarized raman scattering and polarized fluorescence,” J. Polym. Sci., Polym. Phys. Ed. 10, 2135–2153 (1972).
  * Purvis and Bower (1976) J. Purvis and D. I. Bower, “Molecular orientation in poly(ethylene terephthalate) by means of laser–raman spectroscopy,” J. Polym. Sci., Polym. Phys. Ed. 14, 1461–1484 (1976).
  * Archer, Huang, and Fuller (1994) L. A. Archer, K. Huang, and G. G. Fuller, “Orientation dynamics of a polymer melt studied by polarization-modulated laser Raman scattering,” J. Rheol. 38, 1101–1125 (1994).
  * Huang, Archer, and Fuller (1996b) K. Huang, L. A. Archer, and G. G. Fuller, “Microstructural dynamics of a homopolymer melt investigated using two-dimensional raman scattering,” Macromolecules 29, 966–972 (1996b).
  * Nadège  _et al._ (2013) B. Nadège, I. Youssef, M.-C. Chevrel, D. Chapron, C. Schrauwen, S. Hoppe, P. Bourson, and A. Durand, “In situ monitoring of styrene polymerization using raman spectroscopy. multi-scale approach of homogeneous and heterogeneous polymerization processes,” J. Raman Spectroscopy 44, 909–915 (2013).
  * Monnerie (1987) L. Monnerie,  _Developments in Oriented Polymers–2_ , 1st ed., edited by I. Ward (Elsevier, New York, 1987).
  * Lapersonne  _et al._ (1989) P. Lapersonne, J. F. Tassin, P. Sergot, L. Monnerie, and G. Le Bourvellec, “Fluorescence polarization characterization of biaxial orientation,” Polymer 30, 1558–1564 (1989).
  * Viovy and Monnerie (1985) J.-L. Viovy and L. Monnerie, “Fluorescence anisotropy technique using synchrotron radiation as a powerful means for studying the orientation correlation functions of polymer chains,” in  _Characterization of Polymers in the Solid State II: Synchrotron Radiation, X-ray Scattering and Electron Microscopy_ , edited by H. H. Kausch and H. G. Zachmann (Springer Berlin Heidelberg, Berlin, Heidelberg, 1985) pp. 99–122.
  * Janeschitz-Kriegl (1983) H. Janeschitz-Kriegl,  _Polymer Melt Rheology and Flow Birefringence_ (Delft University Press, 1983).
  * Kalogrianitis and van Egmond (1997) S. Kalogrianitis and J. W. van Egmond, “Full tensor optical rheometry of polymer fluids,” J. Rheol. 41, 343 (1997).


11institutetext: Department of Applied Mathematics and Computer Science
Technical University of Denmark, 2800 Kgs. Lyngby, Denmark
11email<EMAIL_ADDRESS>http://www.compute.dtu.dk/ 22institutetext: Institute of
Biological Psychiatry, Copenhagen University Hospital MHC Sct. Hans, Roskilde,
Denmark 33institutetext: Neurobiology Research Unit, Rigshospitalet,
Copenhagen, Denmark 44institutetext: Department of Computer Science,
Copenhagen University, Denmark

# Multi-view consensus CNN for 3D facial landmark placement

Rasmus R. Paulsen 11 0000-0003-0647-3215 Kristine Aavild Juhl 11
0000-0003-1434-697X Thilde Marie Haspang 11 Thomas Hansen 22
0000-0001-6703-7762 Melanie Ganz 3344 0000-0002-9120-8098 Gudmundur Einarsson
11 0000-0002-1326-2339

###### Abstract

The rapid increase in the availability of accurate 3D scanning devices has
moved facial recognition and analysis into the 3D domain. 3D facial landmarks
are often used as a simple measure of anatomy and it is crucial to have
accurate algorithms for automatic landmark placement. The current state-of-
the-art approaches have yet to gain from the dramatic increase in performance
reported in human pose tracking and 2D facial landmark placement due to the
use of deep convolutional neural networks (CNN). Development of deep learning
approaches for 3D meshes has given rise to the new subfield called geometric
deep learning, where one topic is the adaptation of meshes for the use of deep
CNNs. In this work, we demonstrate how methods derived from geometric deep
learning, namely multi-view CNNs, can be combined with recent advances in
human pose tracking. The method finds 2D landmark estimates and propagates
this information to 3D space, where a consensus method determines the accurate
3D face landmark position. We utilise the method on a standard 3D face dataset
and show that it outperforms current methods by a large margin. Further, we
demonstrate how models trained on 3D range scans can be used to accurately
place anatomical landmarks in magnetic resonance images.

###### Keywords:

3D Facial Landmarks Multi-View CNN Geometric Deep Learning

##  1 Introduction

3D face recognition and analysis has a long history with important efforts
dating back to work done in the early nineties [18] and with lots of work
published in the early 2000s [4]. Initially, 3D scanning devices were
expensive, complicated to use, and for laser scanning devices it required that
the subject had to be still and have closed eyes for a longer period of time.
However, the availability of 3D scanning devices ranging from highly accurate
clinical devices to consumer class products implemented in mobile devices has
dramatically increased in the last decade. This means that human face
recognition and analysis is moving from the 2D domain to the 3D domain. 3D
morphometric analysis of human faces is an established research topic within
human biology and medicine, where the applications range from 3D analysis of
facial morphology [20] to plastic surgery planning and evaluation [9]. Broadly
speaking, the analysis of facial 3D morphometry is based on either a sparse
set of landmarks that serves as a simple measure of facial anatomy or the
analysis of a dense set of points that are aligned to the 3D faces using a
template matching approach. The approach of matching a dense template can then
be used to solve the task of selecting a few anatomical landmarks as seen in
for example [19, 29].

While a substantial amount of work has been published on automated 3D landmark
placement, only a limited number of publications have drawn on the recent
drastic improvements in human pose estimation based on deep learning. In this
paper, we describe a framework for automated 3D landmarking of facial surfaces
using deep learning techniques from the field of human pose estimation.

A classic approach to finding facial landmarks and face parameterisations
consists of fitting a 3D statistical shape model to either one or several
views of the 3D surface or directly to the surface [24, 2, 42]. These
approaches use a learned statistical deformation model based on both geometry
and texture variations that is able to synthesise faces within a learned low-
dimensional manifold. By computing the residuals between the actual views and
the synthesised face rendering, the optimal parameterisation of the face model
can be found. In the seminal paper [2] this is done in a standard penalised
optimisation framework, while newer approaches cast the parameter optimisation
into a deep learning framework [42]. These methods work well if the learned
model is broad enough to fit all new faces, and can successfully recover the
face pose. Unusual facial expressions or pathologies might confuse the
methods, since they fall outside the learned appearance manifold.

Placing landmarks purely on the basis of surface geometry is described in
[16], where landmarks are placed by computing the correspondence between a
template face and a given face. The correspondence is computed by minimising
bending energy between surface patches of the reference face and the target
face. In [10] a machine learning approach is described where a set of local
geometrical descriptors are extracted from facial scans and used to locate
landmarks. The descriptors include surface curvature similar to what we
propose. The concept of using local 3D shape descriptors to locate landmarks
is also exploited in [30], where a facial landmark model is fitted to
candidate locations found using curvature and local shape derivatives.

In this work, we use facial surfaces from a range scanner, where the data
consist of surfaces with associated textures. We also use surfaces extracted
as iso-surfaces from magnetic resonance (MR) images of the human head, where
the face is only represented with its geometry. From a geometric point of
view, the face is then a 2D surface embedded in 3D space and it is
topologically equivalent to a disc. Recently, the analysis of this type of
data using deep learning has seen a drastic increase under the term _geometric
deep learning_ [5], where one focus is the transformation of the
representation of, for example, triangulated surfaces into a canonical
representation that is suited for convolutional neural networks [17]. Spectral
analysis of surfaces is one approach to the required domain adaptation as
described in, for example, [3]. Another approach is to use a volumetric
representation where the surface is embedded in a 3D volume [31] and the data
is transformed into a volumetric occupancy grid. The goal of the method
described in [31] is to classify entire objects into a set of predefined
classes. Volumetric methods are still hampered by the drastic memory
requirements for true 3D processing. In [31] the volume size is restricted to
30x30x30, thus severely limiting the spatial resolution. An alternative
approach is to render the surface or scene from multiple views and use a
standard image based convolutional neural network (CNN) on the rendered views.
This is described and analysed in detail in [31, 37], the conclusion being
that with the current memory limits and CNN architectures, multi-view
approaches outperform volumetric methods. Multi-view approaches are
conceptually similar to picking up an object and turning it around while
looking (with one eye) for features or to identify the object. In [38], multi-
view CNNs are used for silhouette prediction with convincing results. However,
there is rapid development in all approaches such as the voxel based method in
[33], which can also estimate the orientation of the object in question.

Figure 1: Overall system overview. A 3D facial scan is rendered from several
views that are fed to a CNN. After a 3D view ray voting process, the result is
accurately placed 3D landmarks on the face surface.

In this paper, we propose a multi-view approach to identify feature points on
facial surfaces. From one direction, a feature point can be seen as a ray in
3D space. By combining results from several views, a 3D landmark can be
estimated as a consensus between feature rays, similarly to classic approaches
from multi-view geometry [21]. A complete overview of the system can be seen
in Fig. 1. A similar idea is presented in [15], where a depth image from a
hand is synthesised into three orthogonal views. A CNN is then trained to
generate 2D heatmaps of hand landmark locations for each of these three views.
The final 3D landmark locations are found by fusing the heatmap probability
distributions as a set of 3D Gaussian distributions and applying a learned
parameterised hand-pose subspace. We propose using more 3D projections (in
this paper, 100) and an outlier-robust method to fuse the individual heatmap
results. Compared to the network used in [15], we propose a deeper network
based on the stacked hourglass network that currently gives state-of-the-art
results for human pose tracking and human face landmark detection [26, 8, 39,
7].

The inspiration for this work comes from recent advances in human pose
tracking, where very deep convolutional networks have been trained to identify
feature points on humans in a variety of poses and environments [26]. The
method is based on heatmap regression, where each individual landmark is coded
as a heatmap. The idea behind the used hourglass network is based on an
aggregation of local evidence [6, 8], where heatmaps of individual landmarks
are created in the first part of the network and fed into the next layer, thus
enabling the network to refine its knowledge of the spatial coherence of
landmark patterns. A recent paper [23], estimates 3D landmarks from 2D photos
using the joint correspondence between frontal and profile landmarks and that
method gets state-of-the-art results on standard test sets used for facial
tasks. We propose a network architecture resembling the architecture described
in [23] but our end goal is different.

In our work, the metric is the landmark localisation error measured in
physical units (in this case, mm) since the end application is often related
to physical estimation of facial morphology. However, in most state-of-the-art
methods in landmark placement based on facial 2D photos, the accuracy is
measured as Normalised Mean Error (NME), which is the average of landmark
errors normalised by the bounding box size [23], which is not a particular
good measure when working with 3D surfaces. In the work of [8, 23], 3D
landmarks are estimated purely from 2D photos by regressing the unknown
z-coordinate using a CNN network. The landmark accuracy in these cases is
highly dependent on the heatmap resolution. They also do not utilise the
availability of a true underlying 3D surface. In our work, we demonstrate a
novel way to combine the output from state-of-the-art networks in a geometric
based consensus to produce highly accurate 3D landmark predictions.

##  2 Methods

A summary of the method can be seen in Fig. 1. The overall idea is that the 3D
surface is rendered from multiple views and for each view, landmark candidates
are estimated. Each landmark estimate is now considered as a ray in 3D space
and the final landmark 3D position is estimated using an outlier robust least
squares estimate between landmark rays.

###  2.1 Multi-view rendering

The scans are rendered using an OpenGL rendering pipeline. The scan is placed
approximately at the origin of the coordinate system and the camera is placed
in 100 random positions around the face, and with the focal point at the
origin. The parameters of the camera are determined so the entire scan is in
view. We use an orthographic projection. For each view, we render a view of
the scan and store the 2D coordinates of the projected ground-truth 3D
landmarks. A simple ambient white light source is used. All surfaces are
rendered in a monochrome non-texture setup and surfaces with texture are also
rendered with full RGB texture colours.

It is an ill-posed problem to pre-align a general surface scan to a canonical
direction as also described in [37]. However, for facial surfaces it is common
to start the pipeline by identifying feature points such as the eyes and the
nose and using them to pre-align the scan or to crop an area of interest as in
[19]. In this work, we do not rely on pre-alignment and simply have the loose
assumption that the scanned surface contains a face among other information,
such as shoulders, and that an approximate direction of the face is given by
the scanner. It also means that we do not assume that a given camera position
can be in any fixed position with relation to the facial anatomy. This makes
the algorithm general and not specific to facial anatomies. The rendering
pipeline is used in both generating 2D training images and computing the
landmarks on an unseen 3D surface.

###  2.2 Geometric Derivatives

To enable the implicit use of geometry, a set of geometric representations are
also rendered. The first is a distance map representation that is computed as
being the OpenGL z-buffer, where the precision has been optimised by setting
the near and far clipping planes as close as possible to the scans bounding
box. Using a depth map for human feature recognition was popularised in the
seminal articles on pose recognition using depth sensors [35] and later
applied in, for example, 3D estimation of face geometry from 2D photos [34].
The standard geometry-only view is also rendered by disabling the texture.
This is the representation used in the multi-view papers [31, 37].

While curvature is implicitly represented in the depth map, we also render a
view where the surface is grayscale coded according to the local mean
curvature. Our aim is to use the method on surfaces containing surface noise,
and therefore a robust curvature estimation is needed. Traditional methods
that estimate curvature based on 1-ring neighbours of a vertex are too noisy.
We use an approach where for each vertex $P$, the algorithm finds the
curvature of the sphere that passes through $P$ and that best approximates the
set of neighbours $P_{N}$ of $P$. The curvature is estimated by first doing an
inverse projection of neighbouring points, so they are lying on an approximate
plane, and then performing an eigenvector analysis of this point cloud. In
this work, neighbour points $P_{N}$ are found by a mesh based region growing
algorithm that includes points connected to $P$ at a maximal distance
empirically chosen to be 10 mm. Proofs and implementation details can be found
in [12, 28].

###  2.3 Network Architecture and Loss Function

We use the two-stack hourglass model described in [26], which is based on the
residual blocks described in [22]. We focus on having a higher resolution of
the predicted heatmap than in previous work. After the heatmap prediction in
the hourglass block, the heatmap is upsampled twice using nearest neighbours
followed by a 3x3 learnable convolutional layers as suggested in [27]. The
input to the network for a single view is 2D renderings of the 3D surface. The
network is flexible with the number of input layers. Using renderings of the
textured surface (RGB) adds three layers, while the depth rendering, the
geometry rendering and the curvature rendering each add one layer. The input
images are 256 x 256, the dimensions throughout the hourglass stacks are 128 x
128 and the heatmap is upsampled to 256 x 256. The entire network can be seen
in Fig. 2. The ground truth is one heatmap per landmark, with a Gaussian
kernel placed at the projected 2D position of the landmark. A cross entropy
loss function is used. The heatmap estimates from the first and second
hourglass modules are concatenated together to form a combined loss function,
as demonstrated in previous work [26]. This ensures intermediate supervision
of the network. Only one network is used and is able to recognise landmarks
from all view directions.

Figure 2: Network architecture. The input is an image with a varying number of
channels (NC) and the output is one heatmap per landmark. NL is the number of
landmarks. Blue boxes are feature maps and the number on top is the number of
feature maps. The spatial size of the layer is written in the lower left
corner (where necessary).

The network is implemented with Tensorflow and trained using RMSPROP with an
initial learning rate of 0.00025, a decay of 0.96 and a decay step of 2000.
The batch size is between 4 and 8, the drop-out rate is 0.02. The network
converged after around 60-100 epochs depending on the used rendering input.
The network was trained and evaluated on one NVIDIA Titan X GPU card. The
network is not optimised with respect to processing time and currently, it
takes approximately 20 seconds to process 100 renderings, including overhead
for reading and storing intermediate results.

###  2.4 Landmark Detection and Consensus Estimation

For a given input image, the output of the network are $NL$ heatmaps, where
$NL$ is the number of landmarks. A 2D landmark is found as the position of the
maximum value in the associated heatmap. This is illustrated in the top box of
Fig. 1. However, to avoid using landmarks that are obviously not located
correctly, only landmarks belonging to heatmap maxima over a certain threshold
are considered. We have experimentally chosen a threshold of $0.5$. When a 2D
landmark candidate has been found, the parameters of the corresponding 3D ray
are computed using the inverse camera matrix used (and stored) for that
rendering. Details of camera geometry and view rays can be found in [21]. This
is illustrated in the bottom part of Fig. 1 (3D view rays). For a given
landmark, this results in up to $NV$ rays in 3D space, where $NV$ is the
number of rendered views.

When the set of potential landmark rays for a given landmark has been
computed, the landmark can be found as the crossing of these rays. In
practice, the rays will not meet in a single point and some of the rays will
be outliers due to incorrect 2D landmark detections. In order to robustly
estimate a 3D point from several potentially noisy rays, we use a least
squares (LSQ) fit combined with RANSAC selection [14]. When each ray is
defined by an origin $\mathbf{a}_{i}$ and a unit direction vector
$\mathbf{n}_{i}$, the sum of squared distances from a point $\mathbf{p}$ to
all rays is:

$\sum_{i}\mathbf{d}^{2}_{i}=\sum_{i}\Big{[}(\mathbf{p}-\mathbf{a}_{i})^{T}(\mathbf{p}-\mathbf{a}_{i})-\big{[}(\mathbf{p}-\mathbf{a}_{i})^{T}\mathbf{n}_{i}\big{]}^{2}\Big{]}$
(1)

Differentiating with respect to $\mathbf{p}$ results in a solution
$\mathbf{p}=\mathbf{S}^{+}\mathbf{C}$, where $\mathbf{S}^{+}$ is the pseudo-
inverse [1] of $\mathbf{S}$. Here
$\mathbf{S}=\sum_{i}(\mathbf{n}_{i}\mathbf{n}^{T}_{i}-\mathbf{I})$ and
$\mathbf{C}=\sum_{i}(\mathbf{n}_{i}\mathbf{n}^{T}_{i}-\mathbf{I})\mathbf{a}_{i}$

In the RANSAC procedure, the initial estimate of $\mathbf{p}$ is based on
three random rays and the residual is computed as the sum of squared distances
from the included rays to the estimated $\mathbf{p}$. The result is a robust
estimate of a 3D point based on a consensus between rays where outlier rays
are not included (Fig. 1, LSQ+RANSAC). The final 3D landmark estimation is
done by projecting the found point to the closest point on the target surface
using an octree based space-division search algorithm. The result can be seen
in Fig. 1(right).

##  3 Data

####  DTU-3D

The dataset consists of facial 3D scans of 601 subjects acquired using a
Canfield Vectra M3 surface scanner. The scanner accuracy is specified to be in
the order of 1.2 mm (triangle edge length). Each face has been annotated with
73 landmarks using the scheme described in [13] and seen in Fig. 3. The faces
in this dataset are all captured with a neutral expression. The data are used
without being cropped, so both ears, neck and shoulders are partly present
(see Fig. 1 for an example scan from the database).

#### BU-3DFE

The database contains 100 subjects (56 female, 44 male). Each subject
performed seven expressions in front of the 3D face scanner, resulting in
2,500 3D textured facial expression models in the database [40]. Each scan is
annotated with 83 landmarks and the faces have been cropped to only contain
the facial region.

#### MR

One Magnetic Resonance (MR) volume of a human head acquired using a 3T Siemens
Trio scanner with a T1-weighted MEMPRAGE sequence with an isotropic voxel size
of 1 mm. We use the N3 algorithm [36] for bias field correction and the
intensity normalisation approach from [11]. The outer skin surface is
extracted using marching cubes [25] with an experimentally chosen iso-surface
value of 20. The surface can be seen in Fig. 4(right); the skin surface is
noisy due to inconsistences in MR values around the skin-air interface.

Figure 3: Left) The 73 landmark scheme used on the DTU-3D database [13].
Right) Landmark localisation errors on the DTU-3D database with different
rendering configurations.

##  4 Results and Discussion

The DTU-3D data were divided into a training set with 541 faces and a
validation set of 60 faces. For both the training and the validation set, the
faces were rendered from 100 camera positions using RGB, geometry, curvature,
and depth rendering. The network is trained and validated using different
render combinations (as seen in Fig. 3). The BU-3DFE was divided into a
training set with all scans belonging to 46 females and 34 males and a
validation set with all scans belonging to 10 females and 10 males. This
results in a training set of 2000 faces and a validation set of 500 faces. The
network using BU-3DFE was trained and evaluated using RGB renderings.

For a given validation face, the 3D landmarks are computed using the proposed
method and compared with the ground truth landmarks. The landmark localisation
error is computed as the Euclidean distance between an annotated landmark and
an estimated landmark. We report the mean and standard deviation of the
landmark localisation error per landmark in millimeters as seen in Tab. 1 and
Fig. 3.

| Salazar [32] | Gilani [16] | Grewe [19] | This paper
---|---|---|---|---
images | 350 | 2500 | 2500 | 500
| Mean | impr. | Mean | impr. | Mean | impr. | Mean | SD
Ex(L) (17) | 9.63 | 73% | 4.43 | 41% | 2.95 | 11% | 2.59 | 1.53
En(L) (21) | 6.75 | 71% | 4.75 | 64% | 3.04 | 37% | 1.89 | 0.98
Ex(R)(25) | 8.49 | 66% | 4.34 | 33% | 3.22 | 11% | 2.85 | 1.5
En(R) (29) | 6.14 | 70% | 3.29 | 33% | 3.23 | 44% | 1.8 | 0.89
Ac(L) (42) | 6.47 | 59% | 4.3 | 38% |  |  | 2.61 | 1.41
Ac(R ) (38) | 7.17 | 58% | 4.28 | 29% |  |  | 2.96 | 1.56
Sn (40) |  |  | 3.9 | 31% | 1.97 | -29% | 2.52 | 1.69
Ch(L) (47) |  |  | 6 | 86% |  |  | 2.182 | 1.44
Ch(R) (53) |  |  | 5.45 | 68% |  |  | 2.42 | 1.44
Ls (50) |  |  | 3.2 | 19% |  |  | 2.33 | 1.31
Li (55) |  |  | 6.9 | 99% |  |  | 2.5 | 1.41
Mean | 7.44 |  | 4.62 |  | 2.88 |  | 2.42 |
Table 1: Results on the BU-3DFE database. Error is in mm. Improvements
compared to previous work are in percentages. The number in parentheses after
the landmark id is the corresponding landmark on Fig. 3.

As can be seen in Fig. 3, the method locates a large set of landmarks with a
localisation error in the range of 2 mm. This is in the limit of what an
experienced operator can achieve [13]. The landmarks that have high errors are
also the landmarks that are typically very difficult for a human to place.
Landmarks placed on the chin and on the eyebrows are very hard to place
manually, due to the weak anatomical cues. An extensive analysis of landmark
errors can be found in [13], where the inter-observer variance is reported for
the landmark set used in the DTU-3D. It was found that the landmark around the
chin has a manual inter-observer in the range of $6mm$. The errors in Fig. 3
for the DTU-3D base is probably more due to inconsistency in the manual
annotation than from the presented method. The method handles uncropped 3D
scans where both partial hair, ears, neck and shoulders are present as seen in
Fig. 1, meaning that the only pre-processing step needed is a rough estimate
of the overall direction of the head.

In Fig. 3 it can be seen that rendering using RGB textures generally performs
very well. This is not surprising since textured surfaces contain many visual
cues. The mode where depth is used in combination with the curvature also
performs very well. This result enables the method to be used on data where
RGB texture is not naturally present, such as iso-surfaces from modalities
like computed tomography (CT) or MR imaging. Using the depth layer alone
yields reasonable results on many landmarks, but in particular around the
chin, where few visual cues are seen in the depth image, the errors are large.
Using geometry rendering yields slightly worse results than the RGB and
depth+curvature renderings, but still on par with the state-of-the-art
algorithms. We also tested the RGB + depth + curvature + geometry
configuration but did not achieve superior results.

The generalisability of the method was tested on the MR scan using depth +
curvature rendering. The method successfully locates landmarks on the MR iso-
surface despite significant noise. The landmarks around the eyes and nose have
an error level in the range of 2-3 mm, while the errors on the chin are
larger, as also seen in the results on the pure range scan. It can be seen
that the estimation of the curvature yields reasonable visual results despite
the very high noise levels on the surface. Placing landmarks in 3D volumes is
notoriously difficult. While deep learning approaches are being used for true
3D landmark localisation [41], the methods are still limited by the
prohibitive GPU memory requirements for true 3D processing. We believe that
our approach offers an alternative way to handle complicated landmark
placement in 3D volumes.

Figure 4: Finding landmarks on iso-surfaces extracted from an MR image using
combined curvature and depth rendering. Left) examples of the rendered
curvature and depth images. Right) the 3D landmarks identified on the skin
iso-surface. The colour coding is the localisation error compared to manually
annotated landmarks.

We have experimented with the number of views, using 25, 50, 75, and 100 view
renderings. Going from 25 to 100 views decreases the landmark localisation
error with less than a millimeter, meaning that in future applications the
number of renderings could be significantly lowered. We have chosen to use
random camera positions instead of selecting a fixed set of pre-defined
positions. We believe that using random positions works as an extra data
augmentation step. Finally, the results on the benchmark data set BU-3DFE as
seen in Tab. 1 demonstrate that the proposed method outperforms state-of-the-
art methods by a large margin.

Compared to the methods that are dominant for 2D photos faces-in-the-wild as
for example [8, 23, 7], our metric is a landmark localisation error in mm,
while they report the error as given as a percentage compared to either the
eye-to-eye distance [7] or the bounding box diagonal length [8, 23]. The
average adult eye-to-eye distance is 62 mm. In [7] they report an error rate
of $4.5\%$ that is equivalent to 2.8 mm. This is a very impressive result but
their measured error is still above what we report.

Since the 3D geometry of the face is available, it is possible to do a view-
ray versus surface intersection test to for example determine if a landmark is
placed on an occluded part of the surface. We have chosen not do that, since
the method in practise is good at predicting the positions of landmarks on
occluded parts of the face due to the high spatial correlation of landmark
positions as for also demonstrated in [8, 23, 7]. Furthermore, the landmarks
with low confidence are removed in the RANSAC procedure and this is often
landmarks that from a given viewpoint are hard to see due to for example
occlusion.

Future works includes improved maxima selection in the heatmap based on more
rigorous statistical assumptions and using a fitting function to locate maxima
with sub-pixel precision. In addition, different renderings and surface
signature techniques to further enhance surface structure can be exploited.
Furthermore, with increased GPU memory in sight, it will be interesting to
design novel network architectures with higher spatial resolutions in mind and
also test other blocks like the inception-resnet blocks used in [23].

##  5 Conclusion

In this paper, we demonstrated a multi-view pipeline to accurately locate 3D
landmarks on facial surfaces. The method outperforms previous methods by a
large margin with respect to landmark localisation error. Furthermore, the
effect of using different rendering types was demonstrated, suggesting that a
model trained on high-resolution 3D face scans could be used directly to
accurately predict landmarks on surfaces from a completely different scanning
modality. It was demonstrated that using a combination of surface curvature
rendering and depth maps performs on par with using RGB texture rendering,
proving that implicit geometric information can be very valuable even when
observed in 2D.

While geometric deep learning is a rapidly growing field and volumetric
methods are gaining foothold, this paper shows that the concept of multi-view
rendering of 3D surfaces currently produces state-of-the-art results with
regard to feature point location.

Code, trained models and test data can be found here:
http://ShapeML.compute.dtu.dk/

## Acknowledgements

We gratefully acknowledge the support of NVIDIA Corporation with the donation
of the Titan Xp GPU used for this research.

## References

  * [1] Ben-Israel, A., Greville, T.N.: Generalized inverses: theory and applications. Springer (2003)
  * [2] Blanz, V., Vetter, T.: A morphable model for the synthesis of 3d faces. In: Proc. Computer graphics and interactive techniques. pp. 187–194 (1999)
  * [3] Boscaini, D., Masci, J., Rodolà, E., Bronstein, M.: Learning shape correspondence with anisotropic convolutional neural networks. In: Proc. NIPS. pp. 3189–3197 (2016)
  * [4] Bowyer, K.W., Chang, K., Flynn, P.: A survey of approaches and challenges in 3d and multi-modal 3d+ 2d face recognition. Computer vision and image understanding 101(1), 1–15 (2006)
  * [5] Bronstein, M.M., Bruna, J., LeCun, Y., Szlam, A., Vandergheynst, P.: Geometric deep learning: going beyond Euclidean data. IEEE Signal Processing Magazine 34(4), 18–42 (2017)
  * [6] Bulat, A., Tzimiropoulos, G.: Convolutional aggregation of local evidence for large pose face alignment. In: Proc. BMVC (2016)
  * [7] Bulat, A., Tzimiropoulos, G.: Two-stage convolutional part heatmap regression for the 1st 3d face alignment in the wild (3dfaw) challenge. In: Proc. ECCV. pp. 616–624 (2016)
  * [8] Bulat, A., Tzimiropoulos, G.: How far are we from solving the 2d & 3d face alignment problem?(and a dataset of 230,000 3d facial landmarks). arXiv preprint arXiv:1703.07332 (2017)
  * [9] Chang, J.B., Small, K.H., Choi, M., Karp, N.S.: Three-dimensional surface imaging in plastic surgery: foundation, practical applications, and beyond. Plastic and reconstructive surgery 135(5), 1295–1304 (2015)
  * [10] Creusot, C., Pears, N., Austin, J.: A machine-learning approach to keypoint detection and landmarking on 3d meshes. International journal of computer vision 102(1-3), 146–179 (2013)
  * [11] Dale, A.M., Fischl, B., Sereno, M.I.: Cortical surface-based analysis: I. segmentation and surface reconstruction. Neuroimage 9(2), 179–194 (1999)
  * [12] Delingette, H.: Modélisation, Déformation et Reconnaissance d’Objets Tridimensionnels á l’Aide de Maillages Simplexes. Ph.D. thesis, L’École Centrale de Paris (1994)
  * [13] Fagertun, J., Harder, S., Rosengren, A., Moeller, C., Werge, T., Paulsen, R.R., Hansen, T.F.: 3d facial landmarks: Inter-operator variability of manual annotation. BMC medical imaging 14(1), 35 (2014)
  * [14] Fischler, M.A., Bolles, R.C.: Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography. In: Readings in computer vision, pp. 726–740. Elsevier (1987)
  * [15] Ge, L., Liang, H., Yuan, J., Thalmann, D.: Robust 3d hand pose estimation in single depth images: From single-view cnn to multi-view cnns. In: Proc. CVPR. pp. 3593–3601 (2016)
  * [16] Gilani, S.Z., Shafait, F., Mian, A.: Shape-based automatic detection of a large number of 3d facial landmarks. In: Proc. CVPR. pp. 4639–4648. IEEE (2015)
  * [17] Goodfellow, I., Bengio, Y., Courville, A., Bengio, Y.: Deep learning. MIT press (2016)
  * [18] Gordon, G.G.: Face recognition based on depth and curvature features. In: Proceedings 1992 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. pp. 808–810. IEEE (1992)
  * [19] Grewe, C.M., Zachow, S.: Fully automated and highly accurate dense correspondence for facial surfaces. In: Proc. ECCV. pp. 552–568. Springer (2016)
  * [20] Hammond, P., et al.: 3d analysis of facial morphology. American journal of medical genetics Part A 126(4), 339–348 (2004)
  * [21] Hartley, R., Zisserman, A.: Multiple view geometry in computer vision. Cambridge university press (2003)
  * [22] He, K., Zhang, X., Ren, S., Sun, J.: Deep residual learning for image recognition. In: Proc. CVPR. pp. 770–778 (2016)
  * [23] J. Deng, Y. Zhou, S.C., Zafeiriou, S.: Cascade multi-view hourglass model for robust 3d face alignment. In: FG (2018)
  * [24] Jourabloo, A., Liu, X.: Pose-invariant 3d face alignment. In: Proc. ICCV. pp. 3694–3702 (2015)
  * [25] Lorensen, W.E., Cline, H.E.: Marching cubes: A high resolution 3d surface construction algorithm. In: ACM siggraph computer graphics. vol. 21, pp. 163–169. ACM (1987)
  * [26] Newell, A., Yang, K., Deng, J.: Stacked hourglass networks for human pose estimation. In: Proc. ECCV. pp. 483–499. Springer (2016)
  * [27] Odena, A., Dumoulin, V., Olah, C.: Deconvolution and checkerboard artifacts. Distill (2016). https://doi.org/10.23915/distill.00003
  * [28] Paulsen, R.R.: Statistical shape analysis of the human ear canal with application to in-the-ear hearing aid design. Ph.D. thesis, Technical University of Denmark (2004)
  * [29] Paulsen, R.R., Marstal, K.K., Laugesen, S., Harder, S.: Creating ultra dense point correspondence over the entire human head. In: Proc. SCIA. pp. 438–447. Springer (2017)
  * [30] Perakis, P., Passalis, G., Theoharis, T., Kakadiaris, I.A.: 3d facial landmark detection under large yaw and expression variations. IEEE transactions on pattern analysis and machine intelligence 35(7), 1552–1564 (2013)
  * [31] Qi, C.R., Su, H., Nießner, M., Dai, A., Yan, M., Guibas, L.J.: Volumetric and multi-view cnns for object classification on 3d data. In: Proc. CVPR. pp. 5648–5656 (2016)
  * [32] Salazar, A., Wuhrer, S., Shu, C., Prieto, F.: Fully automatic expression-invariant face correspondence. Machine Vision and Applications 25(4), 859–879 (2014)
  * [33] Sedaghat, N., Zolfaghari, M., Amiri, E., Brox, T.: Orientation-boosted voxel nets for 3d object recognition. In: British Machine Vision Conference (BMVC) (2017)
  * [34] Sela, M., Richardson, E., Kimmel, R.: Unrestricted facial geometry reconstruction using image-to-image translation. arXiv (2017)
  * [35] Shotton, J., Fitzgibbon, A., Cook, M., Sharp, T., Finocchio, M., Moore, R., Kipman, A., Blake, A.: Real-time human pose recognition in parts from single depth images. In: Proc. CVPR. pp. 1297–1304 (2011)
  * [36] Sled, J.G., Zijdenbos, A.P., Evans, A.C.: A nonparametric method for automatic correction of intensity nonuniformity in mri data. IEEE transactions on medical imaging 17(1), 87–97 (1998)
  * [37] Su, H., Maji, S., Kalogerakis, E., Learned-Miller, E.: Multi-view convolutional neural networks for 3d shape recognition. In: Proceedings of the IEEE international conference on computer vision. pp. 945–953 (2015)
  * [38] Wiles, O., Zisserman, A.: Silnet: Single-and multi-view reconstruction by learning from silhouettes. In: Proc. BMVC (2017)
  * [39] Yang, J., Liu, Q., Zhang, K.: Stacked hourglass network for robust facial landmark localisation. In: Proc. CVPR. pp. 2025–2033. IEEE (2017)
  * [40] Yin, L., Wei, X., Sun, Y., Wang, J., Rosato, M.J.: A 3d facial expression database for facial behavior research. In: Proc. FGR. pp. 211–216. IEEE (2006)
  * [41] Zheng, Y., Liu, D., Georgescu, B., Nguyen, H., Comaniciu, D.: 3d deep learning for efficient and robust landmark detection in volumetric data. In: Proc. MICCAI. pp. 565–572. Springer (2015)
  * [42] Zhu, X., Lei, Z., Liu, X., Shi, H., Li, S.Z.: Face alignment across large poses: A 3d solution. In: Proc. CVPR. pp. 146–155 (2016)