\section{Introduction}
\label{sec:intro}

In the last years flow networks have been defined and analyzed
to understand properties of fluid transport and mixing
\cite{sergiacomi2015flow,sergiacomi2015most,sergiacomi2015dominant,lindner2017spatio},
with particular emphasis on geophysical transport and its
biological implications
\cite{treml2008,jacobi2012identification,rossi2014hydrodynamic,dubois2016linking}.
The basis of the flow-network paradigm rests on addressing the
Lagrangian fluid dynamics of the system with the so-called
set-oriented methods
\cite{Froyland2003,Froyland2005,Froyland2007,Dellnitz2009,
Froyland2010,Santitissadeekorn2010,Froyland2012Rossi,Tallapragada2013,Bolltbook2013}
in which a discrete approximation to the motion of fluid is
characterized by a transfer or Perron-Frobenius operator
indicating which proportion of fluid is moved from one region
to another by the flow. Under the name of `mapping method'
these ideas have provided insight on mixing and its
optimization in the context of industrial flows
\cite{galaktionov2002mapping,singh2008optimizing,singh2008mapping}.
Besides the context of fluid flows, the framework is also
relevant to characterize the behavior of general dynamical
systems
\cite{froyland2003detecting,dellnitz2006graph,santitissadeekorn2007identifying,Bolltbook2013}.
The contact with network theory \cite{Newman2010} is made when
the matrix representing the transfer operator is interpreted as
the adjacency matrix of a network
\cite{dellnitz2006graph,santitissadeekorn2007identifying,Bolltbook2013,
sergiacomi2015flow,sergiacomi2015most,sergiacomi2015dominant,lindner2017spatio},
so that the weight of a link between two nodes is given by the
amount of flow between the corresponding spatial locations. An
alternative viewpoint in characterizing fluid flows with
network techniques assigns links between spatial regions
according to statistical correlations between their dynamical
variables, leading to correlation-based flow networks
\cite{molkenthin2014a,tupikina2016correlation,molkenthin2017edge}
and making a connection with the broad field of climate
networks
\cite{tsonis2004,yamasaki2008climate,donges2009a,barreiro2011}.

Flow networks have been analyzed with standard measures of
network topology, which have complemented the direct study with
specific fluid dynamics methods. For example, degree and
related properties of the different nodes have been related to
local mixing and dispersion, with explicit correspondence with
Lyapunov exponents and entropic measures
\cite{sergiacomi2015flow}. Betweenness centrality highlights
preferred transit nodes connecting distant regions
\cite{sergiacomi2015most}. Closeness and eigenvector centrality
distinguish regions dominated by laminar or by strong mixing,
and identify structures related to invariant manifolds
\cite{lindner2017spatio}. Spectral and other community
detection methods have been used to extract coherent regions
and basins of attractions in the fluid flow
\cite{froyland2007detection,froyland2010transport,jacobi2012identification,froyland2014how,
rossi2014hydrodynamic,sergiacomi2015flow}.

There is however a standard network measure, the
\emph{clustering coefficient}, which has not yet been
interpreted in the context of fluid flows. The aim of the
present paper is to provide such interpretation. We show that,
under some conditions, clustering characterizes in a useful way
periodic trajectories in fluids and in dynamical systems,
serving as a spectroscopic tool that detects the presence of
periodic or close-to-periodic fluid paths of given period, and
reveals their approximate location in space. Periodic orbits
are objects of fundamental importance in the theory of
dynamical systems\cite{ChaosBook} and in the context of fluid
flow they identify recurrent fluid-element motions. Even when
the conditions required to locate periodic trajectories are not
perfectly fulfilled, the clustering coefficient still reveals
interesting flow structures.

To achieve our aim we will consider examples representative of
flows with an increasing level of complexity. After revising
the methodology of network construction (Sect.
\ref{sec:flownet}) and discussing the general relationship
between clustering and closed paths in networks (Sect.
\ref{sec:clustering}), we address in Subsect.
\ref{subsec:steady2d} the simple situation of two-dimensional
(2d) steady flows. In this case the associated flow network is
also steady and the clustering measures have a clear meaning.
The interpretation of these quantities are discussed for 2d
periodic flows in Subsect. \ref{subsec:periodic2d} and for a
three-dimensional (3d) time-independent dynamical system in
Subsect. \ref{subsec:steady3d}. A situation of fully aperiodic
flow is briefly addressed in Subsect. \ref{subsec:aperiodic}.
The paper closes with a sections mentioning further extensions
(Sect. \ref{sec:extensions}) and the Conclusions section (Sect.
\ref{sec:conclusions}).

\section{Flow network construction from fluid motion}
\label{sec:flownet}

Fluid flow is a process occurring in continuous space. A
network representation of it requires some kind of
discretization or {\sl coarse-graining}. We briefly remind here
the basic steps needed to build discrete flow networks from the
Lagrangian fluid
dynamics\cite{sergiacomi2015flow,lindner2017spatio}.

First, the spatial domain occupied by the fluid is discretized
in a large number $N$ of boxes, $\{B_i, i=1,...,N\}$, so that
network node $j$ will represent the fluid box $B_j$. Although
this is not strictly necessary, we take here for simplicity all
boxes to have the same area or volume. Also, unless otherwise
stated, we consider an incompressible flow of constant density.
Then, all boxes contain the same amount of fluid.

To complete the definition of the transport network, we
establish a directional \emph{link} between two nodes when an
exchange of fluid occurred between the corresponding fluid
boxes during a given time interval. To do this, we use the
fluid {\sl flow map} $\Phi_{t_0}^\tau$:
\begin{equation}
{\bf x}(t_0+\tau)=\Phi_{t_0}^\tau ({\bf x}_0)
\label{map}
\end{equation}
which gives the position at time $t_0+\tau$ of the fluid
particle started at ${\bf x}_0$ at time $t_0$, and that is obtained
simply by integrating the equations of motion of the fluid
elements in the velocity field $\dot {\bf x}(t)={\bf v}({\bf x}(t),t)$. By
considering this Lagrangian motion of all fluid particles
inside a region $A$ we define the action of the fluid map on
whole regions: $A(t_0+\tau)=\Phi_{t_0}^\tau (A(t_0))$. By
applying the flow map to our discrete boxes, we obtain the
amount of flow among each pair of nodes, which will be taken as
the {\sl weight} of the corresponding link. More explicitly,
the fraction of the fluid that started at node $i$ at time
$t_0$ and ended up at node $j$ after a time $\tau$ is given by
\cite{Froyland2003,Froyland2005,Froyland2007,Dellnitz2009,Froyland2010,Santitissadeekorn2010}:
\begin{equation}
\mathbf{P}(t_0,\tau)_{ij} = \frac{m\left(B_i \cap
\Phi_{t_0+\tau}^{-\tau}(B_j)\right)}{m(B_i)} \ . \label{PF}
\end{equation}
The transfer or transport matrix $\mathbf{P}(t_0,\tau)_{ij}$ is a
discrete approximation to the Perron-Frobenius operator of the
flow. $m(A)$ is a measure assigned to the set $A$, which in our
case is its area or volume (and, for constant density, the
amount of fluid it contains). A probabilistic interpretation of
Eq. (\ref{PF}) is that $\mathbf{P}(t_0,\tau)_{ij}$ is the probability
for a particle to reach the box $B_j$, under the condition that
it started from a uniformly random position within box $B_i$.
We interpret $\mathbf{P}(t_0,\tau)$ as the adjacency matrix of a
weighted and directed network, so that $\mathbf{P}(t_0,\tau)_{ij}$ is
the weight of the link from node $i$ to node $j$, for given
$t_0$ and $\tau$. $\mathbf{P}(t_0,\tau)$ is row-stochastic, i.e. it
has non-negative elements and $\sum_{j=1}^{N}
\mathbf{P}(t_0,\tau)_{ij}=1$. For incompressible flows it is also
doubly stochastic.

A standard numerical estimation of (\ref{PF}) is obtained by
randomly placing at time $t_0$ a large number of particles
$N_i$ inside each box $B_i$ (with our equal area/volume choice,
all $N_i$ should be equal), and determine with the flow map in
which boxes are the final positions at time $t_0+\tau$. Then:
\begin{equation}
\mathbf{P}(t_0,\tau)_{ij} \approx \frac{\textrm{number of particles
from box $i$ to box $j$}}{N_i} \ . \label{PFparticles}
\end{equation}

Networks constructed from $\mathbf{P}(t_0,\tau)$ are \emph{weighted}
and \emph{directed} \cite{Newman2010}. In addition, the
dependence on $t_0$ gives in fact a different network for each
starting time $t_0$. We can consider the sequence
$\{\mathbf{P}(t,\tau), t=t_0,t_1,...\}$ as a \emph{temporal} network
\cite{holme2012temporal}, i.e. a network which is changing in
time. The case in which $\mathbf{P}(t_0,\tau)$ is independent on the
initial time $t_0$ gives a \emph{static} network. This occurs
when the flow velocity field is steady: ${\bf v}({\bf x},t)={\bf v}({\bf x})$.

It is
standard\cite{dellnitz2009seasonal,Froyland2012Rossi,froyland2014how}
to make a Markovian approximation for the dynamics at times
beyond $\tau$: we assume that the transport matrix at long
times is well approximated by a product of matrices at shorter
times, $\mathbf{P}(t_0,n\tau)=\mathbf{P}(t_0,\tau)\mathbf{P}(t_0+\tau,\tau) ...
\mathbf{P}(t_0+(n-1)\tau,\tau)$. This amounts to assuming that at the
beginning of each time interval the fluid particles are
reinitialized with uniform density in each box, losing memory
of their earlier positions. The assumption is equivalent to
introducing some artificial diffusion in the dynamics
\cite{froyland2013analytic}. In the limit of very small boxes
and very short time steps, the approximation becomes
correct\cite{froyland2013analytic}, as this computational
diffusion vanishes and we recover the exact Lagrangian dynamics
given by (\ref{map}). For finite boxes the added diffusion can
be thought as modeling unresolved scales of motion.

A path in a network is an ordered set of nodes $\{i,j,k,...\}$
such that there is a non-vanishing link between each
consecutive pair of nodes. A \emph{time-respecting
path}\cite{holme2012temporal} can be defined for a temporal
network $\{\mathbf{P}(t,\tau), t=t_0,t_1,...\}$ as a sequence of nodes
$\{i,j,k,...\}$ such that there is a non-vanishing link between
$i$ and $j$ for the first time interval $[t_0,t_1]$, a
non-vanishing link between $j$ and $k$ for the second interval
$[t_1,t_2]$, etc. Time-respecting paths characterize transport
in the flow network
\cite{lentz2013unfolding,sergiacomi2015most,sergiacomi2015dominant}.
Since nodes are finite regions of flow, time-respecting paths
can be though as `tubes' of fluid trajectories giving a
coarse-grained version of them. When no confusion could arise
we would refer to time-respecting paths also as
\emph{trajectories}. Of course, for static networks there is no
distinction between paths and time-respecting paths.

Many properties of networks defined by the transport matrix
$\mathbf{P}(t_0,\tau)$ have been already analyzed, as summarized in
the Introduction. We now consider the clustering coefficient.

\section{Clustering coefficient and closed paths}
\label{sec:clustering}

The clustering coefficient of a node measures the amount of
closed triangles in the network of which that node is a vertex.
\cite{newman2003structure,newman2009networks}. Depending on the
type of network (weighted, directed, ...) and of the kind of
triangles one is interested in, different clustering
coefficients can be
defined\cite{saramaki2007generalizations,fagiolo2007clustering}.
Triangles in a graph give closed paths and then it is natural
to look for relationships between these triangles and closed
trajectories of the fluid elements, the signature of periodic
orbits of the associated advection dynamical system. An
essential characteristic of the fluid trajectories is its
directionality. Thus, considering the \emph{directed} character
of the flow network would be an important ingredient to take
into account when defining a clustering coefficient of
relevance in the study of fluid flows. To focus more clearly on
this relevant aspect we will neglect in the following the
weighted character of flow networks, by considering the
unweighted network defined by the adjacency matrices
$\mathbf{A}(t_0,\tau)$, binary versions of the transport matrices
$\mathbf{P}(t_0,\tau)$:
\begin{eqnarray}
\mathbf{A}(t_0,\tau)_{ij}=1 \ \ \ &\textrm{if} \ &\mathbf{P}(t_0,\tau)_{ij}>0
 \  \textrm{and} \  i \neq j\nonumber \\
\mathbf{A}(t_0,\tau)_{ij}=0 \ \ \ &\textrm{if} \
&\mathbf{P}(t_0,\tau)_{ij}=0\
 \  \textrm{or} \  i=j \ . \label{binary}
\end{eqnarray}
$\mathbf{A}(t_0,\tau)$ still keeps the connectivity and directional
information, as well as the time-evolving character through the
time-dependence on $t_0$. Let us consider first the single
network defined by $\mathbf{A}(t_0,\tau)$ for fixed $t_0$. Not all
triangles present there are related to closed paths. Figure
\ref{fig:triangles} shows four possible configurations of the
directions in the edges of a triangle involving nodes $A$, $B$
and $C$. Only configuration a), which is called a
\textit{cyclic} triangle \cite{fagiolo2007clustering}, gives a
triangle that can be related to a closed path $A \rightarrow B
\rightarrow C \rightarrow A$. Given a node $i$ with in-degree
$k_i^{IN}=\sum_{j\ne i} \mathbf{A}_{ji}$, out-degree
$k_i^{OUT}=\sum_{j\ne i} \mathbf{A}_{ij}$ and with $k_i^B$ of these
links being bidirectional ($k_i^B=\sum_{j\ne i} \mathbf{A}_{ij}
\mathbf{A}_{ji}$), a \emph{cyclic clustering coefficient} $C^c_i$ is
defined as the ratio of all cyclic triangles involving node $i$
present in the network, divided by all possible cyclic
triangles that could have been constructed with these values of
$k_i^{IN}$, $k_i^{OUT}$, $k_i^B$. It can be
computed\cite{fagiolo2007clustering} from the diagonal elements
of the third power of the adjacency matrix $\mathbf{A}$ :
\begin{equation}
C^c_i = \frac{\left(\mathbf{A}^3\right)_{ii}}{k_i^{IN}
k_i^{OUT}-k_i^B} \ . \label{cc}
\end{equation}

\begin{figure}
\centering
\includegraphics[width=.8\columnwidth, clip=true]{triangles.pdf}
\caption{Four distinct configurations of the directions of the edges in a triangle involving
nodes $A$, $B$ and $C$. Case a) is a cyclic triangle. Node $A$ in configuration
c) acts as a source, whereas it is a sink in configuration d). $A$ is just an intermediate
node in configuration b). } \label{fig:triangles}
\end{figure}

With the particular normalization in Eq. (\ref{cc}), $C_i^c \in
[0,1]$. But, since we have already neglected weight
information, we will not pay much attention to the particular
value of $C_i^c$. Rather the interesting point is if it is zero
or not. In the last case there is at least one directed
triangle involving $i$ in the network. What is the meaning of
this when the network represents a fluid flow, with binary
adjacency matrix $\mathbf{A}(t_0,\tau)$? It means that there are at
least two intermediate nodes $m$, $n$ such that during the time
interval $[t_0,t_0+\tau]$, some fluid particles have been
transported from $i$ to $m$, some from $m$ to $n$, and some
from $n$ to $i$. Note that then, the values of $C_m^c$ and
$C_n^c$ are also non-zero.

In the case of a steady flow, the flow network will be
time-independent, and the dependence on the initial time $t_0$
in the transport and adjacency matrices will disappear:
$\mathbf{P}(t_0,\tau)=\mathbf{P}(\tau)$, $\mathbf{A}(t_0,\tau)=\mathbf{A}(\tau)$. In this
situation (under the Markovian assumption mentioned above),
there will be fluid particles actually following the cycle
$i\rightarrow m \rightarrow n$ during a time interval of
duration $3\tau$. Nodes $i$ with non-vanishing $C_i^c$ will be
part of time-respecting paths, i.e. coarse-grained fluid
trajectories, of period $3\tau$. Since fluid trajectories are
continuous and different points on them can be associated with
different values of $t_0$, the invariance of $\mathbf{A}$ with respect
to $t_0$ implies that the nodes with non-vanishing clustering
will arrange in (thick) lines providing approximations to
periodic orbits of period $3\tau$.

The situation is very similar if the flow is periodic with the
same period $\tau$ as the duration of the time-interval used to
construct the network (or a submultiple of it). In this case
$\mathbf{P}(t_0+\tau,\tau)=\mathbf{P}(t_0,\tau)$ (the same is valid for
$\mathbf{A}(t_0,\tau)$) and the repeated action at consecutive times
of the transport process is given by powers of a single matrix
$\mathbf{P}(t_0,\tau)$. At variance with the static case, the
dependence on $t_0$ remains (and it is periodic). Nodes of high
clustering will again give an approximation to the position at
time $t_0$ of periodic orbits of period $3\tau$. In general
these nodes will not form continuous lines, instead the nodes
with high clustering will move following the full trajectory by
changing $t_0$.

If the integration time $\tau$ defining the flow network does
not coincide with the period of the flow (or a multiple of it)
the nodes with a high clustering coefficient are more difficult
to interpret. It is still true that there is motion of fluid
particles between nodes $i\rightarrow n \rightarrow m$ during
the time interval $[t_0+\tau,\tau]$. But after this time the
new transport matrix $\mathbf{P}(t_0+\tau,\tau)$ will be different
from $\mathbf{P}(t_0,\tau)$ and then there is no guarantee that the
same motion will repeat three times to produce a periodic
time-respecting path of period $3\tau$. The situation is indeed
similar to that of a fully aperiodic or turbulent flow, which
keeps an arbitrary $t_0$ dependence on $\mathbf{P}(t_0,\tau)$ and
$\mathbf{A}(t_0,\tau)$.

\section{Application to specific flows}
\label{sec:specific}

\subsection{Steady two-dimensional flow} \label{subsec:steady2d}

\begin{figure*}
\centering
\includegraphics[width=\textwidth, clip=true]{gyre_tau2_3_4_6_eps0.pdf}
\caption{Clustering coefficient $C_i^c$ in the nodes of the
flow network constructed from the steady ($\epsilon=0$) double
gyre flow with $A=0.1$. a) $\tau=2$, b) $\tau=3$, c) $\tau=4$,
d) $\tau=5$. Panel c) also shows two examples of cyclic
triangles present in the system. } \label{fig:steady}
\end{figure*}

In this section and in the following we consider an analytic
model flow, the \emph{double-gyre}, which provides a convenient
workbench to relate fluid dynamics quantities to network
characteristics. See for example
\cite{shadden2005definition,farazmand2012computing} for basic
properties of this system and computations of its Lagrangian
coherent structures and Lyapunov fields. Flow networks have
already been constructed from it
\cite{sergiacomi2015dominant,lindner2017spatio}. The
double-gyre is a two-dimensional time-periodic flow defined in
the rectangular region of the plane, ${\bf x}=(x,y) \in
[0,2]\times[0,1]$. It is described by the streamfunction
\begin{equation}
\psi(x,y,t)= A \sin(\pi f(x,t))\sin(\pi y) \ ,
\label{dg-stream1}
\end{equation}
with
\begin{eqnarray}
f(x,t)&=& a(t)x^2+b(t)x \\
a(t)  &=& \epsilon \sin(\omega t) \ , \\
b(t)  &=& 1-2\epsilon \sin(\omega t) \ . \label{dg-stream2}
\end{eqnarray}
From these expressions, the velocity field is
\begin{eqnarray}
\dot x &=& -\frac{\partial\psi}{\partial y}=-\pi A \sin(\pi f(x,t))\cos(\pi y)\\
\dot y &=& \frac{\partial\psi}{\partial x}=\pi A \cos(\pi
f(x,t)) \sin(\pi y) \frac{\partial f(x,t)}{\partial x} \ .
\label{dg-v}
\end{eqnarray}
We take $A=0.1$. We consider in this subsection $\epsilon=0$,
for which the flow is steady, $f(x,t)=x$, and the terms
containing the periodic forcing are absent. This provides a
simple test case to check the behavior of the clustering
coefficient. Ideal fluid particles follow very simple
trajectories: they rotate following closed streamlines,
clockwise in the left half of the rectangle, and
counterclockwise in the right one. The central streamline
$x=1$, a heteroclinic connection between the hyperbolic point
at $(1,1)$ and the one at $(1,0)$, acts as a separatrix between
the two regions. There are hyperbolic fixed points at
$(x,y)=(0,0),(1,0),(2,0),(0,1),(1,1),(1,2)$ and elliptic fixed
points at $(x,y)=(0.5,0.5)$ and $(1.5,0.5)$. The frequency of
the orbits close to the elliptic points can be obtained by
linearization, giving $\omega_0=\pi^2 A$, from which the period
is $T_0=2/(\pi A)\approx 6.366$. The orbits arrange
concentrically around each elliptic point, with period
increasing from $T_0$ when increasing the distance to the
elliptic point and diverging for the nearly-square orbits that
are close to touch the domain walls and the central $x=1$
separatrix.

We discretize the fluid domain into $100\times 50 = 5000$
square boxes, so defining the nodes in our flow network (each
node thus represents a square region of size $0.02\times
0.02$). The transport matrix $\mathbf{P}(t_0,\tau)=\mathbf{P}(\tau)$ is
computed by releasing 400 particles from each of the boxes. We
analyze the resulting clustering of each node for different
values of $\tau$.

Figure \ref{fig:steady} clearly shows how the clustering
coefficient computed from a flow network of interval $\tau$
highlights the periodic orbits of period close to $3\tau$:  For
small $\tau$ (not shown) clustering is essentially zero
everywhere. When $\tau=2$ clustering begins to be nonzero for
nodes close to the origin, revealing the presence there of
periodic orbits of period approximately $3\times 2=6$. The
exact linear period close to the origin is $T_0=6.366$,
increasing outwards. This kind of `resonance' at $3\tau=6$ is
broad, revealing the orbits for $\tau$ somehow larger and
smaller than the one corresponding to the exact value
$T_0/3=2.122$. When increasing $\tau$, orbits further from the
origin are highlighted, corresponding to the larger period. For
$\tau$ sufficiently large, orbits with period submultiple of
$3\tau$ are also detected (the origin again for $\tau\approx
4$, a second orbit with period close to $7.5$ for $\tau=5$,
etc.)

\subsection{Periodic two-dimensional flow}
\label{subsec:periodic2d}

\begin{figure}
\centering
\includegraphics[width=\columnwidth, clip=true]{CvsEps_tau5.pdf}
\caption{Mean clustering coefficient per node $\left< C^c \right>$ as a function of forcing intensity
$\epsilon$ for the flow network with $A=0.1$, $t_0=0$ and $\tau=5$ equal to the flow
period $2\pi/\omega$.}
\label{fig:SpectTau5}
\end{figure}

\begin{figure*}
\centering
\includegraphics[width=\textwidth, clip=true]{ClustPoincareTau5_eps0p12_0p09.pdf}
\caption{Clustering coefficient (panels a) and c)) and
Poincar\'{e} section at $t_0=0$ (panels b) and d)) for the flow
network constructed from the double gyre with $A=0.1$ and
$2\pi/\omega=\tau=5$. Panels a) and b): $\epsilon=0.12$. Panels
c) and d): $\epsilon=0.09$.}
\label{fig:periodic}
\end{figure*}

We consider now the periodically forced double gyre, i.e. Eqs.
(\ref{dg-stream1})-(\ref{dg-v}) with $\epsilon>0$. The
remaining parameters are the same as in the previous subsection
with the addition of the forcing frequency, which we take to be
$\omega=2\pi/5$ so that the flow period is $T=5$, and the
initial time that we always take as $t_0=0$. The procedure of
discretization and particle releasing needed to define the flow
network are also the same as in the previous Subsect.
\ref{subsec:steady2d}. When $\epsilon > 0$, complex behavior
including chaotic trajectories arises. The periodic
perturbation breaks the central separatrix and now some
interchange of fluid is possible between the left and the right
parts of the domain. The geometric structures involved in this
interchange have been studied with a variety of techniques
\cite{shadden2005definition,farazmand2012computing}.

The cleanest interpretation of the clustering coefficient is
obtained when the flow network is computed during one period of
the forcing, i.e. for $\tau=T=5$, so that
$\mathbf{P}(t_0+\tau,\tau)=\mathbf{P}(t_0,\tau)$ and a single matrix (and
thus a single network) describes the full dynamics. In this
case clustering is expected to be high in the regions
surrounding the locations of the orbits of period $3\tau=3T$,
i.e. the period-3 orbits. Submultiples of this period will also
lead to non-vanishing clustering, which includes the case of
period-1 orbits. The exact period-1 cycle is indeed excluded,
because self-loops into a node are eliminated by the definition
of the adjacency matrix $\mathbf{A}(t_0,\tau)$ in Eq. (\ref{binary}),
but still we expect the tube of trajectories surrounding the
periodic one, arising from the coarse-grained nature of the
network description, to contain paths that connect to
neighboring nodes and return some mass to the original one
after three iterations, giving non-vanishing clustering.

In order to monitor the changes in the clustering spatial
distributions as parameters are varied, we compute a global
quantity, $\left<C^c\right>$, the mean clustering per node. We
display it in Fig. \ref{fig:SpectTau5} as a function of
$\epsilon$. We see a marked peak occurring close to
$\epsilon\approx 0.13$. Figure \ref{fig:periodic} shows the
clustering coefficient of the different nodes in the network
(panel a)) for $\epsilon=0.12$, a value close to the maximum of
$\left<C^c\right>$, and compares it with the Poincar\'{e}
section (panel b)) of the same flow at the same time $t_0=0$.
We see that the dominant features in both the Poincar\'{e}
representation and the clustering spatial distribution are the
locations at $t_0=0$ of two period-1 orbits, each one
surrounded by the three elements of an elliptic period-3 orbit
(separated by the elements of the corresponding hyperbolic
orbits of period-3, which have also non-zero $C^c_i$). For
other parameters of the forcing these orbits are also present,
as shown in Fig. \ref{fig:periodic}c and d. Its footprint in
the coarse-grained phase space defined by the network is
however weaker as indicated by the smaller mean
$\left<C^c\right>$ in Fig. \ref{fig:SpectTau5}.

\begin{figure}
\centering
\includegraphics[width=\columnwidth, clip=true]{CvsTau_eps0p12.pdf}
\caption{Mean clustering coefficient per node $\left< C^c \right>$ as a function of $\tau$ for
the flow networks obtained from the double gyre with $A=0.1$, forced periodically with $2\pi/\omega=5$,
$t_0=0$ and $\epsilon=0.12$.} \label{fig:SpectEps012}
\end{figure}

When the period of the flow does not coincide with the period
$\tau$ used to define the network we are in a quasiperiodic
situation that we expect to be similar to the case of aperiodic
flows: we do not expect exact periodic orbits to be present.
Nevertheless, we can still explore the behavior of the
clustering with changing $\tau$. Again transport of particles
between regions $A$, $B$, $C$ following the pattern $A
\rightarrow B \rightarrow C \rightarrow A$ during the time
interval $[t_0,t_0+\tau]$ will lead to non-vanishing
clustering, although we can not be sure that they will give
periodic time-respecting paths after three iterations since
flow will not be the same in the next intervals
$[t_0+n\tau,t_0+(n+1)\tau]$. They will be approximate periodic
paths if the flow changes slowly on the time scale $\tau$.
Figure \ref{fig:SpectEps012} shows the mean clustering for
fixed $\epsilon=0.12$ and changing $\tau$. There is a sharp
peak at $\tau=5$, as expected, since as shown before there are
cycles $A\rightarrow B \rightarrow C \rightarrow A$ which
indeed lead to periodic trajectories. But we see also high
values of clustering for $2 \lesssim \tau \lesssim 2.6$
indicating that cyclic triangles also occur during the interval
$[t_0,t_0+\tau]$ for these values of $\tau$. Figure
\ref{fig:Eps0p12Tau2p2} shows that, certainly, these regions
containing triangles are well visible in the system for
$\tau=2.2$. Their locations indicate that they are remnants of
the $\epsilon=0$ periodic orbits. They are no longer periodic,
but between $[t_0,t_0+\tau]$ with $t_0=0$ there is cyclic
interchange of particles between these nodes.

\begin{figure}
\centering
\includegraphics[width=\columnwidth, clip=true]{dgyro_tau2p2_eps0p12.pdf}
\caption{Clustering coefficient of the nodes in the flow network of the double gyre with the parameters
of Fig. \ref{fig:SpectEps012} and $\tau=2.2$.}
\label{fig:Eps0p12Tau2p2}
\end{figure}

\subsection{Steady three-dimensional flow}
\label{subsec:steady3d}

\begin{figure}
\centering
\includegraphics[width=\columnwidth, clip=true]{CvsTauLorenz.pdf}
\caption{Mean clustering coefficient per node $\left< C^c \right>$ as a
function of $\tau$ for the Lorenz system.}
\label{fig:SpectLorenz}
\end{figure}

For steady three-dimensional flows (and in general for steady
$d$-dimensional flows), triangles $A\rightarrow B \rightarrow C
\rightarrow A$ identified during an interval $[t_0,t_0+\tau]$
will remain as triangles in successive intervals thus leading
to periodic orbits of period $3\tau$. This is similar to the
previously studied case of steady 2d flow but, at variance with
it, here the periodic orbits are generically isolated and they
could coexist in phase space with more complex chaotic
trajectories. In order to explore a situation markedly distinct
from the 2d steady double-gyre considered in Sect.
\ref{subsec:steady2d}, which is an incompressible flow, here we
consider a well known dissipative dynamical system, the Lorenz
flow \cite{lorenz1963deterministic,sparrow1982}. The
trajectories ${\bf x}(t)=(x(t),y(t),z(t))$ are the solutions of the
system of equations:
\begin{eqnarray}
\dot x &=& s (y-x)\nonumber \\
\dot y &=& r x - y -x z \nonumber \\
\dot z &=& -b z + x y\ .
\label{LorenzEqs}
\end{eqnarray}

We discretize the $(x,y,z) \in
[-20,20]\times[30,30]\times[0,50]$ portion of phase space into
$31616$ cubic boxes, so defining the nodes in our flow network
(each node thus represents a cubic region of lateral size
approximately $1.59$ units). The transport matrix
$\mathbf{P}(t_0,\tau)$ is computed by releasing 1000 particles from
each of the boxes. We take the well studied parameter set
$s=10$, $b=8/3$, $r=28$ for which the chaotic Lorenz attractor
is well developed and dominates the dynamics. In our initial
condition all boxes have the same density, but under evolution
the particles concentrate in the Lorenz attractor, leaving the
rest of phase space empty. The presence of cyclic triangles is
then expected only in the neighborhood of the attractor.
Embedded in it there is an infinity of (unstable) periodic
orbits of different periods
\cite{sparrow1982,viswanath2003symbolic}. These orbits have
been computed by a variety of techniques
\cite{sparrow1982,viswanath2003symbolic,galiasshort2008} and
are expected to give rise to non-vanishing clustering for the
flow network of $\tau$ given by one third of the corresponding
period. Figure \ref{fig:SpectLorenz} shows the mean clustering
per node for the networks constructed from the Lorenz system at
different values of $\tau$. This provides a kind of
`spectroscopy' by which the presence of orbits of different
periods is revealed. We see peaks at $\tau_1=0.52$,
$\tau_2=0.77$, $\tau_3=1.03$, etc. which correspond very well
(within the step sizes we have used to explore $\tau$) to
one-third of the period of the shortest periodic orbits of the
system for these parameter values
\cite{viswanath2003symbolic,galiasshort2008}: $3\tau_1=1.56$
resonates with the periodic orbit LR (meaning one turn to the
left region of the Lorenz attractor and another turn right,
repeated periodically) of period $T_1\approx 1.5587$. The
resonance at $3\tau_2=2.31$ reveals the LLR and RRL orbits,
both with period $T_2\approx 2.3059$. Next resonance, at
$3\tau_3=3.09$ locates the presence of the orbits LLLR and
RRRL, with $T_3\approx 3.0236$, and LLRR, with $T_4\approx
3.0843$ (our scan in $\tau$ and the node discretization has not
been fine enough to distinguish between them).

\begin{figure*}
\centering
\includegraphics[width=\textwidth, clip=true]{LorenzOrbits.pdf}
\caption{Clustering coefficient for the flow network
constructed for the Lorenz model at $\tau=0.52$ (left) and
$\tau=0.77$ (right). Also plotted (yellow thick lines partially
visible inside the red regions) are the periodic orbit LR
(left) and the two periodic orbits LLR and RRL (right). To make
visible all locations with nonvanishing $C_i^c$, nodes with
clustering coefficient larger than the maximum indicated in the
color bar are plotted with the same color.}
\label{fig:LorenzOrbits}
\end{figure*}

To show that certainly the peaks in Fig. \ref{fig:SpectLorenz}
correspond to `resonances' of the value of $\tau$ with these
periodic orbits of period $3\tau$ we show in Fig.
\ref{fig:LorenzOrbits} the clustering coefficient of the
different nodes together with the periodic orbits present at
these values of $\tau$. It is seen how the clustering
coefficient identifies in a coarse-grained way the location of
the periodic orbits.

\subsection{Aperiodic flow}
\label{subsec:aperiodic}

\begin{figure}
\centering
\includegraphics[width=\columnwidth, clip=true]{CvsTauOcean.pdf}
\caption{Mean clustering coefficient per node $\left< C^c \right>$
as a function of $\tau$ for the flow networks constructed from the MFS surface velocity field
in the Mediterranean starting at $t_0=$ January 1st 2010.}
\label{fig:SpectOcean}
\end{figure}

\begin{figure*}
\centering
\includegraphics[width=\textwidth, clip=true]{vfinal_tau6_clustering_startdate01012010_reds.pdf}
\caption{Clustering coefficient of the nodes of the flow
network constructed from the MFS surface velocity field in the
Mediterranean with starting date January 1st 2010 and $\tau=6$
days.} \label{fig:Ocean}
\end{figure*}

To explore the behavior of clustering in an aperiodic flow we
use a surface velocity field obtained from the Mediterranean
Forecasting System (MFS), an hydrodynamic model based on
NEMO-OPA (Nucleus for European Modelling of the
Ocean-PArallelis\cite{madec2008nemo}, version 3.2) at
$1/16^\circ$ horizontal resolution and 72 unevenly spaced
vertical levels\cite{oddo2009nested}. We consider the
horizontal motion on the third layer below the surface
corresponding to about 10 meters depth, which gives a good
description of surface transport not excessively affected by
wind stress.

The flow network is constructed
\cite{rossi2014hydrodynamic,sergiacomi2015flow,sergiacomi2015most}
by discretizing the continuous ocean into quasi-square boxes of
$1/8^\circ$ horizontal size using an area-preserving sinusoidal
projection. Each node is filled at $t_0$ with 100 ideal fluid
particles, and their trajectories integrated during a time
$\tau$ to get $\mathbf{P}(t_0,\tau)$.

Figure \ref{fig:SpectOcean} shows a very broad `resonance' for
$3 \leq \tau \leq 16$ days which identifies potential cyclic
motions of periods $9-48$ days, the range characteristic of
energetic mesoscale structures such as eddies. Figure
\ref{fig:Ocean} displays the spatial distribution of the
clustering coefficient for $\tau=6$ days. The structures and
their locations are reminiscent of ocean eddies, specially the
ones that are semipermanent\cite{millot2005circulation}, such
as the Alboran gyre (just east of Gibraltar strait), the Ionian
anticyclones (east of Sicily) or the Shikmona gyre (south of
Cyprus). Note that several of these objects have a concentric
structure reminiscent of the nested periodic orbits of Fig.
\ref{fig:steady}c,d or \ref{fig:Eps0p12Tau2p2}. For larger
$\tau$ the structures become much weaker (smaller $C^c_i$) but
also larger, suggestive of the big regional gyres in the
Adriatic, the Tyrrhenian , etc. We do not expect perfect
periodic orbits to be present in the turbulent particle motion
produced by this ocean velocity field, but the small and slow
displacements of the semipermanent
gyres\cite{millot2005circulation} detected by the clustering
coefficient (Fig. \ref{fig:Ocean}) makes very likely that they
contain nearly-periodic circulations with approximate period in
a range including $3\tau=18$ days. In fact, these structures
have\cite{millot2005circulation} radii of about 60-150 $km$ and
typical circulations speeds of 0.2-0.6 $m/s$ , leading to
circulation times in the range 7-40 days, consistent with the
periods associated to the broad maximum of Fig.
\ref{fig:SpectOcean}. In any case non-vanishing clustering
reveals cyclic motions in the interval $[t_0,t_0+\tau]$.

\section{Further extensions}
\label{sec:extensions}

Standard clustering coefficients count triangles. But one can
also count the number of rectangles, pentagons, etc. adjacent
to a given node in a network. We have checked explicitly in the
double-gyre flow example that, as expected, the generalized
cluster coefficients highlight orbits of period $4\tau$,
$5\tau$, etc. An example of this is shown in Fig. \ref{fig:C4}
for the double-gyre flow. We display a generalized clustering
coefficient $C_4$ which is the number of cyclic polygons of 4
sides passing through each node. Panel a), constructed with
$\tau=3$ for the steady case, highlights the same orbit of
period $4\tau=12$ as in Fig. \ref{fig:steady}c. Panel b)
displays $C_4$ for $\epsilon=0.07$ and flow period
$2\pi/\omega$ equal to the integration time $\tau=5$. The
locations of the central period-1 orbits are identified, as
well as four spots around each of them corresponding to the two
elements of elliptic and of hyperbolic period-2 orbits. There
are however other features of difficult interpretation,
apparently associated to the boundaries and to the central
vertical separatrix. Also the central period-1 orbit is not
apparent in Fig. \ref{fig:C4}a). We conclude that the
generalized clustering coefficients contain in principle
similar information to the standard triangular clustering, but
in a more imprecise way arising probably from the extension of
the Markov approximation to longer cycles.

\begin{figure}
\centering
\includegraphics[width=\columnwidth, clip=true]{C4lt50.pdf}
\caption{Generalized clustering coefficient, $C_4$, counting the number of cyclic polygons
of 4 sides at each node, for the double-gyre flow. a) Steady case, $\epsilon=0$, $\tau=3$;
an orbit with period $4 \tau=12$ is highlighted. b) Periodic case with $\epsilon=0.07$
and forcing period $2\pi/\omega=5=\tau$. }
\label{fig:C4}
\end{figure}

An alternative characterization of triangular flow patterns
which does not take into account the directional nature of flow
can be done with the \emph{undirected} cluster coefficient
\cite{newman2009networks,fagiolo2007clustering}. It is defined
by constructing a new \emph{undirected} network with adjacency
matrix $\mathbf{U}$ such that $\mathbf{U}_{ij}=1$ if either $\mathbf{A}_{ij}=1$ or
$\mathbf{A}_{ji}=1$, and else $\mathbf{U}_{ij}=0$. The degree of site $i$ is
now $k_i=\sum_{i\ne j} \mathbf{U}_{ij}$. The undirected clustering
coefficient $C_i^{u}$
is\cite{newman2009networks,fagiolo2007clustering} the number of
triangles in this new network involving site $i$, over the
total number of possible triangles at that node for such degree
($k_i (k_i-1)/2$):
\begin{equation}
C^u_i = \frac{\left(\mathbf{U}^3\right)_{ii}}{k_i (k_i-1)} \ .
\label{cu}
\end{equation}
We have computed $C^u_i$ for several of the flow examples
considered above. For the double gyre case non-vanishing
$C^u_i$ occurs essentially in the same nodes as non-vanishing
$C^c_i$. This indicates that nearly all triangles present in
double-gyre flow networks are of the cycle type (type a) in
Fig. \ref{fig:triangles}). This may be due to the strong
constraints imposed by the incompressibility of the flow. The
situation for the Lorenz system, however, is completely
different: the higher values of $C^{u}_i$ occur in wide tubes
that intersect the Lorenz attractor transversally, and that
seem to trace portions of the stable manifolds of the
non-vanishing unstable fixed points embedded in the attractor.
We interpret this preliminary observation as revealing an
abundance of triangles of type c) in Fig. \ref{fig:triangles}:
different particles released in node $A$, on a region close to
the stable manifold of one of the fixed points, fall onto the
attractor after a time $\tau$ at slightly different places, $B$
and $C$, which are connected by the dynamics in the attractor.
Further work is needed to establish if this observation is a
general feature of dissipative dynamical systems or
compressible flows, so that different clustering coefficients
of the associated flow network would be useful to locate
complex structures beyond periodic orbits.

We finally mention as another future development the full
consideration of the weighted nature of the flow network, i.e.
using matrix $\mathbf{P}$ instead of $\mathbf{A}$. We expect, for example,
that triangle strength would be higher around stable
trajectories as compared to unstable trajectories, since more
fluid particles will remain close to the stable ones.

\section{Conclusions}
\label{sec:conclusions}

In summary, for flow networks with integration time $\tau$
constructed for steady flows or for periodic ones of period
$\tau$ (or submultiple of it), the cyclic clustering
coefficient is non-vanishing in regions locating periodic
orbits of period $3\tau$. Both stable and unstable orbits are
detected. Computation of the transport or adjacency matrices on
which the network approach is based could be expensive (but
easily parallelizable) because of the large number of
Lagrangian particles that should be released. But once
obtained, the matrices can be used to obtain a variety of flow
diagnostics (Lyapunov exponents, community detection, and other
applications mentioned in the Introduction). The computation of
the clustering itself represents just a small computational
overhead that gives access to an important aspect,
periodicities, of flow transport.

Scanning a range of $\tau$ provides a kind of spectroscopic
tool that gives high mean clustering per node at $\tau$ when
periodic orbits of period $3\tau$ exist and have large impact
on the phase space of the system. Even for aperiodic networks
clustering still identifies cyclic triangular paths during the
interval $[0,\tau]$. They could again be interpreted as
periodic trajectories in an approximate way when the time scale
of change of the temporal network is slow compared to $\tau$.
It is likely that a more precise identification of periodic
paths in aperiodic networks could be achieved by replacing the
cyclic clustering coefficient by a \emph{temporal} version of
it, for example by replacing the third power of the single
matrix $\mathbf{A}$ in Eq. (\ref{cc}) by a product of three adjacency
matrices at consecutive times. This line of work is worth to be
explored in the future.

At variance with standard methods to compute periodic
trajectories in dynamical systems\cite{Nayfeh1995}, clustering
in networks highlights fluid nodes, which are finite regions,
and then it detects thick tubes of trajectories instead of
single pathways. Thus, this network approach can not compete
with standard dynamical systems methods when \emph{precision}
in the location of periodic trajectories is required.
Nevertheless the coarse-graining implied by the network
discretization and the Markov approximation is equivalent to
perturbing trajectories with diffusing
noise\cite{froyland2013analytic} and thus our thick network
paths would be very relevant objects when dealing with
transport and dispersion in the presence of diffusion and/or
noise. More in general, we expect that the connections among
dynamical systems and network theory here proposed would open
additional theoretical developments and fruitful transfer of
concepts and methods across these two fields.

\begin{acknowledgments}
We acknowledge financial support from Spanish Ministerio de
Econom\'{\i}a y Competitividad and Fondo Europeo de Desarrollo
Regional under grants LAOP, CTM2015-66407-P (MINECO/FEDER) and
ESOTECOS projects FIS2015-63628-C2-1-R (MINECO/FEDER), FIS2015-
63628-C2-2-R (MINECO/FEDER). E.S-G. acknowledges partial
support from the program ``Investissements d'Avenir'' launched
by the French Government  and implemented  by ANR under grants
ANR-10-LABX-54 MEMOLIFE  and ANR-11-IDEX-0001-02 PSL* Research
University.
\end{acknowledgments}


\section{Introduction}
In 1992, Lee Smolin posed the question\cite{smolin1992did}, ``Did the
Universe Evolve?''  His use of the word ``evolve'' was not in the
physicists' sense of dynamical time evolution, but in the biologists'
sense of speciation.  Smolin proposed that gravitational collapse
produces an offspring Universe with slightly different Standard Model
(SM) parameters from the progenitor: a microphysical analogy to the
microbiology of inheritance.  After many generations, he argued, we
should find ourselves at a local maximum over SM parameter space for
the formation of collapsed objects.  In this way, he proposed a
fascinating resolution to the Hierarchy Problem: it is our notion of
naturalness that we must adjust.

Of the criticisms (e.g. \cite{smolin2006status}) levelled against
Smolin's 1992 proposal, the only one relevant is that it
disagrees with observation (e.g. \cite{rothman1993smolin}).  We regard
this result as refutation only of Smolin's particular fitness
function: reproductive rate.  Note that this fitness function is
appropriate in the simplest of scenarios: e.g. bacteria with unlimited
food.  In any realistic setting, one cannot consider the organism in
isolation; competition and environmental pressures ultimately
determine the fate of \emph{populations}.  In fact, it is observation
of how interacting organisms within populations change over timescales
much greater than an individual's lifespan which led to Darwin's
conclusions\cite{darwin1872origin} on fitness, not the other way
around.

The purpose of this paper is to develop a notion of population and
ancestry via an extension to the Einstein-Hilbert action.  We restrict
ourselves to SM clones for simplicity, as we have yet no
\emph{experimental} basis for speculation on the microphysical origins
of speciation.  Our approach combines two well-motivated lines of
extension to General Relativity (GR): multiscalar-tensor theories and
multiple metric theories.  There is a vast and mature literature on
scalar-tensor theories.  Not only have their observational signatures
been rigorously characterised(e.g. \cite{damour1992tensor,
  damour1993nonperturbative}), but scalar fields are well known to
contribute negative pressure, and are a mainstay of primordial
inflation scenarios\cite{bassett2006inflation}.  Multiple metric
theories also have a long history (e.g. \cite{rosen1963flat}), with
well-known pathologies\cite{TEGP, boulanger2001inconsistency} having
motivated novel and potentially viable
(e.g. \cite{hassan2012bimetric}) approaches.  The result of our effort
will be a classical field theory which makes quantitative predictions
amenable to immediate confrontation with existing and future
astrophysical data.

The rest of this paper is organised as follows.  We first review the
necessary technical background and emphasise some common caveats when
working in multiple metric settings. We then motivate and introduce
our gravitational and matter actions, featuring $N$ metric degrees of
freedom.  Through conformal relations, we establish a consistent
notion of causality across the $N$ metrics by replacing $N-1$ metric
degrees of freedom with scalar degrees of freedom.  We perturbatively
expand the resulting gravitational action, which takes the form of a
multiscalar-tensor theory in the Jordan frame, and determine
conditions for dynamical stability.  We compute the equations of
motion, demonstrate the consistency of our branching process, and
examine the Newtonian limit.  We then discuss very many future
directions and conclude.  An appendix presents a conceptually
related model as a foil to emphasise our particular approach.  A
second appendix outlines an alternate route to a model with many of
the features developed in the main body, with differing aesthetic
criteria that may appeal to some researchers.  Throughout the paper,
we employ units where $c \equiv 1$.  Note that we will define notation
at its introduction, and will often abbreviate ``stress-energy'' as
simply ``stress.''  If context should require distinction, it will be
made explicitly.

\section{The phylogenic model}
We briefly review some of the mathematics involved when considering
multiple metrics defined on a common differentiable manifold
$\mathscr{M}$ of dimension $n$.  First recall any suitable almost
everywhere $C^\infty$ differentiable manifold admits the notion of
curves independently of any geometry.  These objects are maps $C(\tau)
: \mathscr{M} \leftarrow \tau \in [a,b] \subset \mathbb{R}$ which,
when composed with the coordinate charts $\xi^i : \mathbb{R}
\leftarrow \mathscr{M}$, become real-valued differentiable functions
$\xi^i \circ C$.  Derivatives of these compositions can then be used,
in a natural way\cite{o1983semi}, to construct tangent and cotangent
spaces on $\mathscr{M}$.

At one's discretion, notions like metric, connection, and volume form
(Jacobian) can be introduced to augment $\mathscr{M}$.  There is no
natural notion of uniqueness for any of these objects.  Given a
metric, however, one may uniquely construct a torsion-free connection
compatible with this metric.  This Levi-Civita connection is the only
type which we consider below.

Observers using distinct metrics will arrive at distinct physical
conclusions.  For example, fix two events $x^\mu$ and ${x'}^\mu$ on
$\mathscr{M}$ and fix some curve $C$ such that $C(\tau_i)\equiv x^\mu$
and $C(\tau_f)\equiv {x'}^\mu$.  As stated above, this curve has a
field of velocity vectors $\rmd C^\mu/\rmd\tau$ that exist
independently of any geometry.  If $C$ is always time-like with
respect to metrics labelled $p$ and $q$, then we may consider the
proper separation as perceived by observers using metric $p$ versus
observers using metric $q$.  It may happen that
\begin{align}
\int_C \sqrt{\xo{p}{g}_{\mu\nu}\frac{\rmd x^\mu}{\rmd\tau}\frac{\rmd x^\nu}{\rmd\tau}}~\rmd\tau \neq \int_C \sqrt{\xo{q}{g}_{\mu\nu}\frac{\rmd x^\mu}{\rmd\tau}\frac{\rmd x^\nu}{\rmd\tau}}~\rmd\tau .
\label{eqn:proper_times}
\end{align}
Physically, the presence of multiple metrics requires an additional
means of assigning ``measurement apparatus'' to observers, and this
assignment partitions observers into distinct classes. For clarity,
``observers'' continues to mean ``frames which arrive at the same
physical conclusions under diffeomorphism.''  We regard this property
as a virtue of multiple metric theories, in that such theories
simultaneously acknowledge the mathematical non-uniqueness of metric
and remove the unique \emph{class} of observers from GR.  This
uniqueness in GR infamously led to Fock's somewhat amusing complaint
of, ``a widespread misinterpretation of the Einsteinian Gravitation
Theory as some kind of general relativity.''\cite{fock2015theory}.

In the multimetric setting, care must be exercised when executing
common tensor operations.  Raising and lowering indices on an
$q$-labelled object can only be accomplished with the $q$ metric.  For
example
\begin{align}
\xo{q}{R_{\mu\nu}}\xo{p}{g^{\mu\rho}} \neq \xo{q}{R^\rho}_\nu
\end{align}
and cannot be simplified further, without additional information relating
the $p$-metric to the $q$-metric.  Likewise there are now distinct
covariant derivatives, which are only compatible with their paired
metrics.  In other words
\begin{align}
\xo{p}{\nabla}_\mu\left(\xo{q}{g^{\mu\nu}}\xo{p}{T_{\nu\rho}}\right) \neq \xo{q}{g^{\mu\nu}}\xo{p}{\nabla}_\mu\xo{p}{T_{\nu\rho}}
\end{align}
because the $q$ metric is not compatible with the $p$ connection
and so cannot be commuted through this covariant differentiation.

\subsection{The phylogenic matter and gravitational actions}
\begin{figure}
\caption{Directed, acyclic graph representation of couplings,
  conformal relations, and relative orders for a model anchored at
  $p$.  Note terms $O(\kappa^2)$ are discarded in the
  nearest-neighbour model of \S\ref{sec:nearest_neighbor_model}.}
\includegraphics{tree_arrows.eps}
\label{fig:tree_schematic}
\end{figure}

The utility of GR and its metric foundation, as a quantitative model
for reality, has been rigorously verified over many decades of
spatiotemporal scale.  Consequently, any departures from GR must
either be very slight, subtly hidden, or both.  Consider a population
of $N$ metrics and $N$ copies of the SM.  Introduce a dimensionless
coupling constant $\kappa > 0$, and consider a directed, acyclic
graph (DAG) e.g. the rooted tree of Figure~\ref{fig:tree_schematic}.
Associate to each vertex $q$ of this tree a distinct metric
$g^q_{\mu\nu}$ and stress $T^q_{\mu\nu}$.  Given any two vertices $p$
and $q$ of this graph, let $r(p,q) \geq 0$ be the graph distance
between vertex $p$ and vertex $q$.  This distance is well-defined
because there is a unique simple path between any two vertices of a
tree.  We assign observers to metrics by augmenting a straightforward
generalisation of the canonical matter action, due to Hohmann and
Wohlfarth\cite{hohmann2010repulsive}, with the graph distance between
$p$ and $q$
\begin{align}
\delta \xo{p}{S}_M \equiv -{1\over 2}\int\rmd^4\xi~ \sum_{q}^N  \kappa^{r(p,q)}\xo{q}{T}_{\mu\nu}\delta\xo{q}{g^{\mu\nu}}\sqrt{q}.
\label{eqn:matter_actions}
\end{align}
Here $\sqrt{q}$ is shorthand for the natural volume form induced by
the $q$ metric.  Note that these $N$ matter actions, indexed by $p$,
each enforce\cite{fock2015theory} $N$ distinct statements of the
Equivalence Principle: observers using the $q$ metric will observe $q$
stress to follow $q$ metric geodesics.

Multiple metric gravitational actions are tightly constrained.  With
the exception of Hassan and Rosen's multiple metric
theory\cite{hassan2012bimetric}, Theorem 1.1 of
Boulanger~\emph{et. al.}  \cite{boulanger2001inconsistency} demands
that all metrics must be decoupled.  We thus consider the class of
models indexed by $p$
\begin{align}
\xo{p}{S}_G \equiv \frac{1}{16\pi G}\int\rmd^4\xi~ \sum_{q}^N  \kappa^{r(p,q)} \xo{q}{R}\sqrt{q}
\label{eqn:grav_actions}
\end{align}
where powers of $\kappa$ enter in the same manner as (\ref{eqn:matter_actions}).
Note that these actions suggest a natural perturbative treatment in
powers of $\kappa$.  If one varies the gravitational actions
(\ref{eqn:grav_actions}) with respect to the metric degrees of
freedom, combination with (\ref{eqn:matter_actions}) recovers $N$
distinct statements of \emph{exactly} Einstein's equations for each
$q$.  Thus, there is no way to experimentally distinguish any of these
actions from that of Hilbert and we have successfully hidden our new
degrees of freedom.

\subsection{Conformal relations and ancestry}
\label{sec:ancestry}
An undesirable feature of any model with $N$ possibly distinct metrics
is $N$ distinct notions of causal ordering.  We resolve this issue,
introduce coupling, and bring the model very close to already
well-established literature by demanding that each vertex's metric be
conformally related to its unique graph parent.  Let $p$ have a set of
children $\{a_1, \dots, a_{n_a}\}$.  Then, for the $i$-th child, there
exists a scalar field $\phi_{a_i}$ and constants $\sigma_{a_i}$ such that
\begin{align}
\xo{a_i}{g}_{\mu\nu} \equiv \exp(\sigma_{a_i}\phi_{a_i})\xo{p}{g}_{\mu\nu}
\label{eqn:conformal_relations}
\end{align}
where we have used the standard exponential parameterisation with the
factor of two removed for economy of notation.  Note that, from any
child's perspective, the parent metric acquires an inverted sign
\begin{align}
\xo{p}{g}_{\mu\nu} = \exp(-\sigma_{a_i}\phi_{a_i})\xo{a_i}{g}_{\mu\nu}.
\label{eqn:sign_inverse}
\end{align}
These $N-1$ scalar fields define the edges of the graph and replace
$N-1$ metric degrees of freedom.

Note carefully that the model is no longer permutation symmetric.  In
other words, for distinct $p$, the resultant models from application
of the extrema principle will yield non-equivalent equations of
motion.  Nevertheless, the \emph{qualitative} behaviour observed in
all such models will be unchanged.  This can be seen by considering a
model where each vertex has a fixed number of children $n_a$ and the
tree is countably infinite.  In this limit, one would have perfect
permutation symmetry.  For a discussion of an absolutely permutation
invariant construction, and why we have avoided this model, we refer
the interested reader to Appendix~\S\ref{sec:foil_model}.  In the
remainder of our discussion, we will focus on a specific model
anchored at $p$.  We will use the word ``native'' to refer to tensors
associated with vertex $p$, and the word ``foreign'' to refer to
tensors associated with other vertices.  Special attention must be
paid to the scalar fields, however, as they define the graph's edges,
and are a property defined between two vertices.

\subsection{Asexual reproduction and scalar potentials}
\label{sec:sex}
Our next task is to implement the Smolin-like branching process, which
he assumes to be associated with black hole (BH) production.  The
approach we will take is distinctly pragmatic, and motivated from
condensed matter systems.  Consider the bulk magnetisation of a
ferromagnetic material as it passes through the Curie Temperature from
above.  Here, a vector degree of freedom enters the dynamics at a
phase transition.  In complete analogy, we assume that a scalar degree
of freedom $\phi_a$ and a tensor degree of freedom $T^a_{\mu\nu}$
enter the dynamics during gravitational collapse.  In this way, we
seek to implement a proposal of Chapline and Laughlin who, a decade
earlier, argued\cite{chapline2003quantum, chapline2001quantum} that
one really should replace the singular collapse scenario with a phase
transition.  We will glue an $N$ vertex model to a separate $N+1$
vertex model, both anchored at $p$, using the final data of the
$N$-vertex model fields to ``bootstrap'' the initial data of the $N+1$
model in a consistent fashion.

Chapline and Laughlin assert that the simplest scenario resulting from
their phase transition should be an interior GR vacuum solution with
anomalously large cosmological constant.  Adapted to our context, in
its simplest form, this can be accomplished by a single scalar field
dominated by potential.  Yet unlike the matter actions, which are
constrained by the Equivalence Principle; and the gravitational
actions, which are constrained by Boulanger, \emph{et. al.}; in the
absence of a comprehensive theory, a potential is explicitly
phenomenological.  There is an extremely rich literature on
inflationary potentials, their observational signatures, and their use
in resolving outstanding problems in particle physics (for an
excellent review, see \cite{bassett2006inflation}).  Inflationary
theory, however, is dominated by particle and quantum-field theoretic
concepts, and any foray would be premature given the present limited
scope.  Thus, we only partially constrain the specific form of $V$,
and instead focus on how potentials should enter our actions.

Using Rosen's notation\cite{rosen1995discrete}, let $[A_{qw}]$ be the
adjacency matrix for the ancestral tree.  Let $d(p,q)$ be the \emph{signed}
relative graph depth (c.f. $r(p,q)$) between vertex $p$ and vertex
$q$, with ancestors negative consistent with (\ref{eqn:sign_inverse}).
We define
\begin{align}
\xo{p}{S}_V \equiv - \frac{1}{8 \pi G}\int\rmd^4\xi~ \sum_{q}^N  \sum_{w > q} \kappa^{r(p,w)}[A_{qw}] V(\phi_w)\sqrt{w}
\label{eqn:pot_actions}
\end{align}
where the notation $w > q$ means to consider pairings such that $w$
satisfies $d(p,w) > d(p,q)$.  Note that $V: \mathcal{F}(\mathscr{M})
\leftarrow \mathcal{F}(\mathscr{M})$ is a fixed map understood to have
dimension of inverse length squared.  A single map is consistent with
the Copernican Principle and our restriction to a population of SM
clones.  In reheating scenarios, energy within the scalar field is
transferred to SM degrees of freedom.  The Equivalence Principle then
suggests an absence of cross-couplings between scalar fields in $V$.
Note again that, assuming no active scalar degrees of freedom, one
just regenerates $N$ distinct cosmological constants.

\subsection{A possible observational signature}
An interesting and immediate consequence of
(\ref{eqn:conformal_relations}) is that all members of the ancestral
tree share \emph{the same} (tensor) gravitational radiation field in
vacuum.  This is because the conformally invariant Weyl tensor is the
only non-zero component of the Riemann tensor in vacuum.  This is
consistent with the presence of a single metric degree of freedom,
whose quantized linearization would supposedly provide the graviton.
Since the radiation field is shared, a significantly higher than
expected observation rate for BH mergers
(e.g. \cite{abbott2016gw151226}) could be the result of foreign
sources.  Given sufficient directionality, at distances near enough to
the Earth for appreciable amplitude, a slightly aspherical supernova
gravitational wave signature coupled with a lack of observation in the
EM and neutrino sectors would be a ``smoking gun'' independent of
gravitational wave polarisation.

\section{The nearest-neighbour approximation}
\label{sec:nearest_neighbor_model}
To begin exploration of the phylogenic model, we work to first order
in $\kappa$, as this simplifies the resultant equations of motion and
facilitates preliminary confrontation with experiment.  To bring the
gravitational Lagrangian into standard multiscalar-tensor form, it is
most natural to work with $p$'s metric.  Let $q$ index over $N-1$
vertices adjacent to $p$.  We define $x$ to be $p$'s parent, and $a$
to be a numeric index over $p$'s children, where we drop numeric
subscripts because no ambiguity can arise in this limited context.  By
the conformal relations (\ref{eqn:conformal_relations}), we have
that
\begin{align}
\xo{p}{g}_{\mu\nu} &\equiv \exp(\sigma_p\phi_p) \xo{x}{g}_{\mu\nu}  &\xo{a}{g}_{\mu\nu} &\equiv \exp(\sigma_a\phi_a)\xo{p}{g}_{\mu\nu} \\
\sqrt{p} &= \exp(2\sigma_p\phi_p)\sqrt{x}  & \sqrt{a} &= \exp(2\sigma_a\phi_a)\sqrt{p}
\end{align}
which is visualised in Figure~\ref{fig:tree_schematic}.
Recall\cite{carroll2004spacetime} that Ricci scalars for conformally
related metrics in $n$ spacetime dimensions satisfy
\begin{align}
\xo{a}{R} = \exp(-\sigma_a\phi_a)\left[\xo{p}{R} + (n-1)\left(\frac{1}{4}\sigma_a^2(2-n)\xo{p}{\nabla^\nu}\phi_a\xo{p}{\nabla_\nu}\phi_a - \sigma_a\xo{p}{\nabla^\nu}\xo{p}{\nabla_\nu} \phi_a\right)\right].
\label{eqn:conformal_ricci}
\end{align}
We then find terms of the form
\begin{align}
\xo{a}{R}\sqrt{a} = \exp(\sigma_a\phi_a)\left[\xo{p}{R} + (n-1)\left(\frac{1}{4}\sigma_a^2(2-n)\xo{p}{\nabla^\nu}\phi_a\xo{p}{\nabla_\nu}\phi_a - \sigma_a\xo{p}{\nabla^\nu}\xo{p}{\nabla_\nu} \phi_a\right)\right]\sqrt{p}
\end{align}
appearing in the gravitational action.  We may now integrate by parts, exploiting
\begin{align}
\exp(\sigma_a\phi_a)\xo{p}{\nabla^\nu}\xo{p}{\nabla_\nu}\phi_a  = \xo{p}{\nabla_\mu}\left[\xo{p}{g^{\mu\nu}}\exp(\sigma_a\phi_a)\xo{p}{\nabla_\nu}\phi_a\right] - \sigma_a\xo{p}{\nabla^\nu}\phi_a\xo{p}{\nabla_\nu}\phi_a\exp(\sigma_a\phi_a)
\end{align}
to remove the second derivatives.  This results in integrand terms of the form
\begin{align}
\xo{a}{R}\sqrt{a} = \exp(\sigma_a\phi_a)\left[\xo{p}{R} + \frac{\sigma_a^2}{4}(n-1)(n-6)\xo{p}{\nabla^\nu}\phi_a\xo{p}{\nabla_\nu}\phi_a\right]\sqrt{p} + \partial\sqrt{p}
\end{align}
where $\partial$ is a total divergence.  Substitution into
(\ref{eqn:grav_actions}), dropping $\partial$, gives the following
gravitational Lagrangian density
\begin{align}
\xo{p(1)}{\mathscr{L}_G} = &\left\{\xo{p}{R}\left[1 + \kappa\sum_q\xo{q}{\exp}(\phi\sigma)\right] - \sum_{q,w} \nabla^\nu\xo{q}{\phi}\nabla_\nu\xo{w}{\phi}\left[\frac{\kappa}{4}(n-1)(n-6)\xo{q}{\sigma^2} \xo{q}{\exp}(\phi\sigma)\delta^{qw}\right]\right\}\sqrt{p} + O(\kappa^2).
\label{eqn:grav_lagrangian_jordan}
\end{align}
We have omitted labels over the covariant derivatives, since they are
always with respect to $p$, and overset the exponential to indicate
that its arguments are regarded as particular to that specific label.
Note that (\ref{eqn:grav_lagrangian_jordan}) is expressed in the
Jordan frame.  Consistent with the notation of
\cite{kuusk2016invariant} we define
\begin{align}
A(\phi_q, \dots) &\equiv 1 + \kappa\sum_q\xo{q}{\exp}(\phi\sigma) + O(\kappa^2)\label{eqn:A}\\
B^{qw}(\phi_q, \dots) &\equiv \frac{\kappa}{4}(n-1)(n-6)\xo{q}{\sigma^2} \xo{q}{\exp}(\phi\sigma)\delta^{qw} + O(\kappa^2)\label{eqn:B}.
\end{align}
Note that terms of order $\kappa^2$ and higher would produce a
non-diagonal matrix in (\ref{eqn:B}).  This is because more distant
``relatives'' would be conformally related to $p$ by terms
$\exp(\sum_q \sigma_q\phi_q)$.  Upon differentiation, these would introduce
mixed products of covariant derivatives depending on the specific
relative.

\subsection{Classical instabilities}
A common pathology in classical field theories is unbounded energy
exchange, also called classical instability.  Independently of the
matter sector, a multiscalar-tensor theory is
unstable\cite{kuusk2016invariant} if
\begin{align}
\frac{B^{qw}}{2A} + \frac{3}{4A^2}\frac{\partial A}{\partial \phi_q}\frac{\partial A}{\partial \phi_w} \geq 0.
\end{align}
Note that the above expression becomes the coefficient on the kinetic
terms when the multiscalar-tensor theory is cast into the Einstein
frame.  Note that $A$ is always non-zero if
\begin{align}
\kappa \ll 1/N
\label{eqn:kappa_constraint}
\end{align}
and that this is equivalent to saying that there are not ``too many''
offspring for the given $\kappa$.  In this scenario, we may multiply
through by $2A$
\begin{align}
B^{qw} + \frac{3}{2A}\frac{\partial A}{\partial \phi_q}\frac{\partial A}{\partial \phi_w} \geq 0.
\label{eqn:ghostbuster}
\end{align}
At order $\kappa$, only the $B^{qw}$ term contributes provided that
the second term remains bounded.  Consider the diagonal where $q=w$.
Immediately, we must have $\forall q$ that
\begin{align}
\sigma_q \equiv
\left\{
\begin{array}{cc}
iy_q & 1 < n < 6 \\
y_q & n \geq 7
\end{array}\right.\qquad y_q \in \mathbb{R}
\end{align}
and we have arrived at a rather curious conclusion.  For spacetime
dimension $n=4$, each metric tensor in the population must be related
to every other by a pure phase.  Paying attention to the real portion
of (\ref{eqn:ghostbuster}), we also require that
\begin{equation}
\phi_q \in \left[-\pi/2, \pi/2\right] + 2\pi m \qquad m \in \mathbb{Z}
\label{eqn:phi_constraint}
\end{equation}
be dynamically enforced.  Investigation of the second term in
(\ref{eqn:ghostbuster}) gives
\begin{align}
-\frac{3\kappa^2\xo{q}{\exp}(2\sigma\phi)\left[1 + \kappa\sum_n\xo{n}{\exp}(-2iy\phi)\right]}{2 + 2\kappa^2 N + 4\kappa^2\sum_{n\neq m}\cos(y_n\phi_n - y_m\phi_m)}
\end{align}
where we have made the denominator real.  This term remains
$O(\kappa^2)$ provided that $\kappa \ll 1/\sqrt{N}$, which is
consistent with (\ref{eqn:kappa_constraint}).  The off-diagonal terms
remain similarly bounded, and can thus be discarded in the present
treatment.  We conclude the nearest-neighbour model is free of
instabilities in spacetime dimension $n=4$, provided that metrics are
conformally related by a pure phase and (\ref{eqn:phi_constraint})
remains satisfied.  The presence of this phase leads to many
stimulating questions of interpretation.  Unfortunately, these
questions would lead us too far from our present scope.  For
simplicity, we now focus on the model got by demanding the real
variation of the action to vanish identically, while permitting the
imaginary portion to float according to the solutions of the real
system.  For discussion of the $n\geq 7$ phylogenic model, which does
not require the introduction of any complex numbers, we refer the
interested reader to Appendix~\S\ref{sec:alternative}.

We must now define our approach to invariant spacetime lengths.  We
proceed by considering the effect of conformally scaling the usual
spatially flat Robertson-Walker (RW) metric
\begin{align}
\xo{p}{\rmd s^2} = -\rmd t^2 + p(t)^2\rmd \vec{x}^2
\end{align}
where $p(t)$ is vertex $p$'s scale factor, as reckoned from observers
within $p$.  Now consider a single child labelled $q$, and choose the
convention that $\sigma_q \equiv +i$.  By
(\ref{eqn:conformal_relations}), we have that
\begin{align}
\xo{q}{\rmd s^2} = -\exp(i\phi_q)\rmd t^2 + \exp(i\phi_q)p(t)^2\rmd \vec{x}^2
\end{align}
the real portion of which is
\begin{align}
\xo{q}{\rmd s^2} = -\cos(\phi_q)\rmd t^2 + \cos(\phi_q)p(t)^2\rmd \vec{x}^2.
\end{align}
Observers in $q$ will proceed to make physical conclusions with \emph{their} most natural time coordinate
\begin{align}
\tau(t) \equiv \int_c^t \sqrt{\cos(\phi_q(t'))}~\rmd t' \implies \rmd \tau = \sqrt{\cos(\phi_q(t))}~\rmd t
\label{eqn:foreign_time}
\end{align}
Thus, observers in $q$ simply perceive the usual RW metric
\begin{align}
\xo{q}{\rmd s^2} \equiv -\rmd \tau^2 + q(\tau)^2\rmd \vec{x}^2 \label{eqn:foreign_rw}
\end{align}
while observers in $p$ ``know better''
\begin{align}
q(\tau) \equiv p\sqrt{\cos(\phi_q)}\label{eqn:foreign_scale}.
\end{align}
Note that our definition (\ref{eqn:foreign_time}) again requires that
(\ref{eqn:phi_constraint}) remain satisfied at all times.  We thus see
that an exchange of the spacelike and timelike trajectories signals
the emergence of a classical instability.

By the Copernican Principle encoded into the gravitational action
(\ref{eqn:grav_actions}), we expect that offspring should begin from a
very tiny spatial point since this is how we appear to have begun.
Given (\ref{eqn:foreign_rw}) and (\ref{eqn:foreign_scale}), this
occurs at some branching time $t_*$, unique to each offspring, if
\begin{align}
\phi_{q}(t_*) \equiv \frac{(2m + 1)\pi}{2} \qquad m\in\mathbb{Z}.
\label{eqn:position_ic}
\end{align}
Near these points, the child's clock is significantly slowed from the
parent's perspective.  For the $n=4$ model, this effect is actually
\emph{symmetric} along the ancestral tree between any two vertices,
since cosine is even.  This is reminiscent of clock rates measured by
observers in relative inertial motion under Special Relativity.

\subsection{Metric equations of motion}
\label{sec:eom}
Having established that the nearest-neighbour model is free of
instabilities if (\ref{eqn:phi_constraint}) remains satisfied, we
proceed to compute the equations of motion in $n=4$ spacetime
dimensions.  Consequently, in all equations below, $\sigma \equiv \pm
i$.  Sign is of course determined by ancestry and $y \equiv 1$ has
been assumed without loss of generality.  Note that care must be taken
to not prematurely divide out exponential factors, as their real
portions may evaluate to zero.  Equations of motion will be found from
requiring that the real portion of
\begin{align}
\delta S \equiv \delta \xo{p}{S}_G + \delta \xo{p}{S}_M + \delta \xo{p}{S}_V
\end{align}
vanish.  Substitution of ancestry (\ref{eqn:conformal_relations}) into
the matter actions (\ref{eqn:matter_actions}) through order $\kappa$
gives
\begin{align}
\delta S_M = -\frac{1}{2}\int \rmd^4\xi~\left\{\left[\xo{p}{T}_{\mu\nu} + \kappa\sum_q\xo{q}{\exp}(\sigma\phi) \xo{q}{T}_{\mu\nu}\right]\delta\xo{p}{g^{\mu\nu}} - \kappa\sum_q\xo{q}{\sigma}\xo{q}{\exp}(2\sigma\phi)\xo{q}{T}\delta\xo{q}{\phi}\right\}\sqrt{p} + O(\kappa^2).
\end{align}
To determine the potential contribution, we first expand the adjacency
sum of (\ref{eqn:pot_actions}) to first order in $\kappa$
\begin{align}
\xo{p}{S}_V = -\frac{1}{8\pi G}\int \rmd^4\xi~\left[ \kappa \xo{p}{V}\sqrt{p} + \kappa\sum_a \xo{q}{V}\sqrt{a}\right] + O(\kappa^2),
\end{align}
substitute ancestry, and perform the variation
\begin{align}
\delta S_V = \frac{\kappa}{16\pi G}\int \rmd^4\xi~\left\{\delta\xo{p}{g^{\mu\nu}}\left[\xo{p}{V}\xo{p}{g_{\mu\nu}} + \xo{p}{g_{\mu\nu}}\sum_q\xo{q}{\exp}(2\sigma\phi)\xo{q}{V}\right] - 2\frac{\partial \xo{p}{V}}{\partial \phi}\delta\xo{p}{\phi} - 2\sum_q\xo{q}{\exp}(2\sigma\phi)\left[\frac{\partial \xo{q}{V}}{\partial \phi} + 2\sigma \xo{q}{V}\right]\delta\xo{q}{\phi}\right\}\sqrt{p}.
\end{align}

We leverage the results of \cite{kuusk2016invariant} to write down the
metric and scalar equations of motion using (\ref{eqn:A}) and
(\ref{eqn:B}).  We find for the mixed metric equation of motion
\begin{align}
\left(1+\kappa \sum_q \xo{q}{\exp}(\sigma\phi)\right)G^\mu_\nu - \kappa\sum_q \xo{q}{\exp}(\sigma\phi) \left[\left(\frac{\nabla^\alpha\xo{q}{\phi} \nabla_\alpha\xo{q}{\phi}}{4} - \xo{q}{\sigma} \nabla^\alpha\nabla_\alpha \xo{q}{\phi} \right)\delta^\mu_\nu  + \frac{\nabla^\mu\xo{q}{\phi}\nabla_\nu\xo{q}{\phi}}{2}  + \xo{q}{\sigma}\nabla^\mu\nabla_\nu\xo{q}{\phi} \right]& \nonumber \\
 - 8\pi G\xo{p}{T^\mu_\nu} +\kappa \xo{p}{V}\delta^\mu_\nu +\kappa\sum_q\xo{q}{\exp}(2\sigma\phi)\left[\xo{q}{V}\delta^\mu_\nu  - 8\pi G\xo{q}{T^\mu_\nu}\right] + O(\kappa^2)& = 0
\label{eqn:metric_eom}
\end{align}
where we are careful to raise the index on the foreign stress after we
have changed to the foreign metric, which introduces a conformal
factor.  At this point, recall the considerations of \S\ref{sec:sex},
where we note that the dynamics must not become acausal when a scalar
field enters during the collapse process.  In order that the effective
gravitational constant remain initially unchanged, we have immediately
that
\begin{align}
\xo{q}{\phi}(\vec{x}, t_*) \equiv \frac{\pi}{2} + 2\pi m \qquad m \in \mathbb{Z}
\label{eqn:birth_ic}
\end{align}
which is nothing more than the generalisation of
(\ref{eqn:position_ic}) to arbitrary metrics.  Since
(\ref{eqn:kappa_constraint}) guarantees that we do not have excessive
offspring, it is useful to divide the metric equations
(\ref{eqn:metric_eom}) by $A(\phi, \dots)$ and Taylor expand through
$O(\kappa)$
\begin{align}
G^\mu_\nu - \kappa\sum_q \xo{q}{\exp}(\sigma\phi) \left[\left(\frac{\nabla^\alpha\xo{q}{\phi} \nabla_\alpha\xo{q}{\phi}}{4} - \xo{q}{\sigma} \nabla^\alpha\nabla_\alpha \xo{q}{\phi} - \xo{q}{\exp}(\sigma\phi)\xo{q}{V}\right)\delta^\mu_\nu  + \frac{\nabla^\mu\xo{q}{\phi}\nabla_\nu\xo{q}{\phi}}{2}  + \xo{q}{\sigma}\nabla^\mu\nabla_\nu\xo{q}{\phi}\right]& \nonumber \\
 - 8\pi G\xo{p}{T^\mu_\nu} + \kappa \xo{p}{V}\delta^\mu_\nu- 8\pi G\kappa\sum_q\xo{q}{\exp}(\sigma\phi)\left[\xo{q}{\exp}(\sigma\phi)\xo{q}{T^\mu_\nu} - \xo{p}{T^\mu_\nu}\right] + O(\kappa^2)& = 0.
\label{eqn:metric_eom_einstein}
\end{align}
This equation is Einstein's equation with a position dependent
gravitational constant (somewhat obscured by grouping), augmented with
scalar kinetic and potential terms.  A new feature due to our matter
actions is the presence of foreign stress, which appears here as
purely dark matter.

\subsection{Scalar equations of motion and screening}
\label{sec:mimetic}
Suppressing $q$ labels on $\phi$ for clarity, since there is no
cross-coupling between the scalar fields at order $\kappa$, the scalar
equation of motion for each $q$ is found to be
\begin{align}
\kappa\xo{q}{\exp}(\sigma\phi)&\left(\xo{q}{\sigma}\xo{p}{R} + 3\nabla^\alpha\nabla_\alpha\phi + \frac{3\xo{q}{\sigma}}{2}\nabla^\alpha\phi\nabla_\alpha\phi \right) + 8\pi G\xo{q}{\sigma}\kappa\xo{q}{\exp}(2\sigma\phi)\xo{q}{T} \nonumber \\
& - 2\kappa\left[\frac{\xo{q}{\partial V}}{\partial \phi}\delta^p_q + (1-\delta^p_q)\xo{q}{\exp}(2\sigma\phi)\left(2\sigma \xo{q}{V} + \frac{\xo{q}{\partial V}}{\partial \phi}\right) \right] + O(\kappa^2) = 0.
\label{eqn:scalar_eom_ricci_child}
\end{align}
We now take the trace of the full metric EOM
(\ref{eqn:metric_eom})
\begin{align}
-\xo{p}{R}\left(1 + \kappa\sum_q\xo{q}{\exp}(\sigma\phi)\right) + &3\kappa\sum_q\xo{q}{\exp}(\sigma\phi)\left[\xo{q}{\sigma}\nabla^a\nabla_\alpha\xo{q}{\phi} - \frac{\nabla^\alpha\xo{q}{\phi}\nabla_\alpha\xo{q}{\phi}}{2}\right] \nonumber \\
&- 8\pi G\xo{p}{T} + 4\kappa\xo{p}{V} + \kappa\sum_q\xo{q}{\exp}(2\sigma\phi)\left[4V(\xo{q}{\phi}) - 8\pi G\xo{q}{T}\right] + O(\kappa^2) = 0.
\label{eqn:eom_trace}
\end{align}
and use it to remove the Ricci scalar to leading order
\begin{align}
\xo{q}{\exp}(\sigma\phi)\left[3\nabla^\alpha\nabla_\alpha\phi + \frac{3\sigma}{2}\nabla^\alpha\phi\nabla_\alpha\phi\right] &= \nonumber \\
8\pi G\xo{q}{\sigma}\xo{q}{\exp}(\sigma\phi)&\left[\xo{p}{T} - \xo{q}{\exp}(\sigma\phi)\xo{q}{T}\right] + 2\left[\frac{\xo{q}{\partial V}}{\partial \phi}\delta^p_q + (1-\delta^p_q) \xo{q}{\exp}(2\sigma\phi)\left(2\xo{q}{\sigma} \xo{q}{V} + \frac{\xo{q}{\partial V}}{\partial \phi}  \right)\right]+ O(\kappa).
\label{eqn:scalar_eom}
\end{align}
Note that the potential contribution distinguishes the scalar
equations, depending on relation within the ancestral tree.  In the
absence of any potential contribution, (\ref{eqn:scalar_eom}) takes
the familiar form of a curved-space wave equation with kinetic drag,
with the marked exception that the source term can \emph{zero-cross}.
Sign alteration of the source about
\begin{align}
\exp(\sigma \phi)\xo{p}{T} = \exp(2\sigma \phi)\xo{q}{T}
\end{align}
will drive the field to act as if it were free.  We are unaware of any
form of screening in the existing literature with this behaviour.
Note that screening is only possible for $\phi \in [0, \pi/4]$ if both
sources respect (or violate) the strong energy condition and
\begin{align}
|\xo{p}{T}| \leq |\xo{q}{T}|.
\end{align}
If $\phi \in [\pi/2, \pi/4)$ one, but not both, of the sources must
  violate the strong energy condition in order screen.  We digress
  momentarily from the $n=4$ theory to emphasise that screening
  persists in the purely real $n \geq 7$ theory.

\subsection{Branching consistency and initial conditions}
Since the equations of motion (\ref{eqn:metric_eom}) obtain from a
scalar action, they are automatically covariantly conserved.  This is
true for both an $N$ vertex model anchored at $p$ and an $N+1$ vertex
model anchored at $p$.  Since transitioning between these two models
at a collapse event must not violate local conservation, we require
that the covariant divergence of the newly introduced terms must
vanish.  Assign the label $q$ to the newly introduced vertex.  Then
consistency requires that
\begin{align}
\sigma\nabla_\mu\phi\Lambda^\mu_\nu + \nabla_\mu\Lambda^\mu_\nu =\Bigg|_{t_*} \exp(\phi \sigma)\left[8\pi G \xo{p}{\nabla}_\mu\xo{q}{T^\mu_\nu} - \nabla_\nu V + 2\sigma\nabla_\mu\phi\left(8\pi G\xo{q}{T^\mu_\nu} - V\delta^\mu_\nu\right)\right]
\label{eqn:cov_birth_n4}
\end{align}
where $\Lambda^\mu_\nu$ is defined in (\ref{eqn:Lambda}) and the $q$
label has been suppressed except in locations of possible confusion.
By inspection, (\ref{eqn:cov_birth_n4}) simplifies
considerably if the scalar field is initially gradient-free
\begin{align}
\left.\nabla_\mu \xo{q}{\phi}(x^\alpha)\right|_{t_*}  \equiv 0.
\label{eqn:velocity_ic}
\end{align}
Note that this initial condition leaves the initial value of the
potential and its derivative unconstrained (so long as they are
finite).  Substitution of (\ref{eqn:velocity_ic}) into
(\ref{eqn:cov_birth_n4}) gives
\begin{align}
\nabla_\mu\Lambda^\mu_\nu =\Bigg|_{t_*} 8\pi G \xo{p}{\nabla}_\mu\xo{q}{T^\mu_\nu}\exp(\phi \sigma)
\end{align}
where the gradient of $V$ has vanished because $V$ produces
real-valued $\mathcal{F}(\mathscr{M})$, whereas re-expression of the
covariant derivative will introduce metric factors.  Focusing on the
real portion removes the covariant accelerations, and using the
contracted Bianchi identity we find
\begin{align}
\frac{1}{4}\nabla_\nu\left[\nabla^\alpha\phi\nabla_\alpha\phi\right] + \frac{1}{2}\nabla_\mu\left[\nabla^\mu\phi\nabla_\nu\phi\right] =\Bigg|_{t_*} 4\pi G \left(\xo{q}{T^\alpha_\nu}\nabla_\alpha\phi  + \xo{q}{T}\nabla_\nu\phi\right)\sin(\phi).
\end{align}
Since each term on the left will contain a gradient of the field after
differentiation, we see covariant conservation is satisfied during
branching given the initial condition (\ref{eqn:velocity_ic}).  That
this initial condition is sufficient should not be surprising: energy
transfer associated with a collapse event should be (initially)
spatiotemporally localised to the collapse.

If we evaluate the parent scalar equation of motion at the initial
condition (\ref{eqn:birth_ic}), we find
\begin{align}
4\pi G\xo{p}{T} =\Bigg|_{t_*} \frac{\partial \xo{p}{V}}{\partial \phi}.
\end{align}
For all known matter and radiation, the left hand side is always less
than or equal to zero.  The analogous relation for an offspring scalar
equation at the initial condition introduces a sign.  Since $V$ is a
fixed map, this means that $V$ must have an extremum at the initial
condition.  Note that during the usual cold \emph{and} warm inflation
scenarios, the left hand side vanishes, which is consistent with the
extremum.  Further, this corresponds to a constant potential in the
neighbourhood of the initial condition, which suggests exponential
inflation.

\section{Newtonian behaviour}
\label{sec:consistency}
In this section, we begin the first investigations of whether the
nearest-neighbour model can reproduce the phenomenology of GR on solar
system scales.  First we will compute the Newtonian limit, and
highlight the promising screening behaviour.  Unfortunately, at this
point, we cannot use the well-known results of Esposito-Farese and
Damour\cite{damour1992tensor} to draw PPN conclusions because we have
multiple stress-energies.

We approximate solutions to scenarios featuring weak-fields
\begin{align}
V(\phi) \sim \epsilon ^2 &\qquad \frac{\partial V}{\partial \phi} \sim \epsilon ^2 \\
\xo{p}{\rho}(x^\alpha) \sim \epsilon ^2 &\qquad \xo{q}{\rho}(x^\alpha) \sim \epsilon^2\\
\xo{p}{g}_{\mu\nu} = \eta_{\mu\nu} + \epsilon^2 \xo{p}{h_{\mu\nu}}(x^\alpha) + O(\epsilon^4) &\qquad \xo{q}{\phi}(x^\alpha) = \xo{q}{\Phi} + \epsilon^2 \xo{q}{\phi}(x^\alpha) + O(\epsilon^4).
\end{align}
This is most appropriate for a parent $p$ considering contributions
from offspring $q$.  The Newtonian regime is additionally
characterised by slow motions
\begin{align}
\frac{\partial}{\partial t} \sim \epsilon \frac{\partial}{\partial x}.
\end{align}
We will regard the scalar contributions as a source, and must be
careful to consider their entire active gravitational mass
i.e. effective mass density plus spatial trace.  Using Will's
notation\cite{TEGP} for the gravitational potential, in this
approximation, combination of (\ref{eqn:metric_eom}) and
(\ref{eqn:scalar_eom}) gives the following Poisson equation after
taking real portions
\begin{align}
\nabla^2 U = 4\pi G \left[\xo{p}{\rho} - \sum_q\kappa\cos(\xo{q}{\Phi})\left(\xo{p}{\rho} - \xo{p}{T}/6\right) \right] + \kappa\sum_q\left\{4\pi G\cos(2\xo{q}{\Phi})\left[\xo{q}{\rho} - \xo{q}{T}/6\right] + \frac{1}{6}\left[\frac{\partial V}{\partial \phi}\sin(2\xo{q}{\Phi}) - \cos(2\xo{q}{\Phi})V\right]\right\}
\label{eqn:newtonian_limit}
\end{align}
Neglecting native pressures, this is Newton's equation with a weakened
gravitational constant, dark matter, and scalar potential.  We now
specialise this result to some relevant limits, and focus only on the
behaviour of a single contributing offspring $q$.

\subsection{Contribution to $p$ from inflating offspring}
We will assume the cold inflationary scenario for our single offspring
$q$.  Of course, if $q$ inflates, we expect $V$ to be of the same
order as the collapsed parent matter, so $V \sim O(1)$.  To get a feel
for what happens in this limit at a distance from the collapsed
object, we formally adjust $\kappa$ to guarantee the validity of the
Newtonian approximation from $p$'s perspective.  Removing all stress
from (\ref{eqn:newtonian_limit}), we find
\begin{align}
\nabla^2 U = \frac{\kappa}{6}\left[\frac{\partial V}{\partial \phi}\sin(2\xo{q}{\Phi}) - \cos(2\xo{q}{\Phi})V\right].
\end{align}
Since at $t_*$ the offspring is effectively
frozen from the parent's perspective, we find that
\begin{align}
\nabla^2 U = \frac{\kappa}{6} V\left(\frac{\pi}{2}\right).
\end{align}
 We emphasise strongly that we have not yet investigated the radial
 behaviour of the nearest-neighbour model in the static, spherically
 symmetric limit.  This limit should be the most relevant after
 collapse due to the relative clock rates between parent and
 offspring, and is required to confirm the expected localised
 quasi-static behaviour of a Chapline-Laughlin BH.

\subsection{Contribution to $p$ from late-stage offspring}
\label{sec:screening}
After the offspring has reheated, we expect the potential to no longer
contribute substantially to the equations of motion.  Since we
(present epoch observers) are stationed in $p$, it is reasonable that
$p$'s pressure be zero.  Due to the difference in clock rates, we
might expect vertex $q$ to be radiation dominated.  If we set $V$ and
$q$'s trace to zero in (\ref{eqn:newtonian_limit}) we find
\begin{align}
\nabla^2 U = 4\pi G \xo{p}{\rho}\left[1 - \frac{7\kappa}{6}\cos(\xo{q}{\Phi})\right] + \kappa4\pi G\cos(2\xo{q}{\Phi})\xo{q}{\rho}.
\label{eqn:newtonian_after_inflation}
\end{align}
Note that since we have assumed $\Phi \in [0, \pi/2]$, the dark matter
term can enter with \emph{repulsive} character if $\Phi > \pi/4$.
Such repulsive matter contributions on large scales are disfavoured by
the lack of deviations between positive mass simulations and
observation in the 3rd moment (and thus sensitive to sign) of the mass
convergence maps from weak-lensing surveys\cite{heymans2012cfhtlens}.
Yet, any means to realise ``negative energy'' without instability is
interesting.  For $\Phi \in [0, \pi/4]$, the scalar field will
oscillate about
\begin{align}
\frac{\cos(\Phi)}{\cos(2\Phi)} = \xo{q}{T}/\xo{p}{T}
\label{eqn:mimetic_condition}
\end{align}
given the constraints discussed in \S\ref{sec:mimetic}.  Since the
source term is essentially a driving force to a harmonic oscillator
with kinetic drag, this suggests damped oscillation about an
equilibrium point within this range.  More so, the equilibrium point
will track the evolving matter densities.  Remarkably, at exactly this
equilibrium point (\ref{eqn:newtonian_limit}) becomes
\begin{align}
\nabla^2 U = 4\pi G \xo{p}{\rho},
\end{align}
which is precisely Newton's equation.

\section{Discussion and future directions}
We have chosen to investigate the $n=4$ theory primarily because we
find ourselves within a 3+1 dimensional spacetime.  For the purpose of
model building, however, the pure phase model is awkward due to the
introduction of complex quantities.  Yet the pure phase model is
demonstrably perturbative in $\kappa$, and this permits calculation of
the condition under which the nearest-neighbour approximation remains
dynamically stable.  Whether this dynamical condition is honoured
provides the first target for subsequent cosmological investigations.

While it may be plausible that the stability condition be honoured by
a single contributing field, the superposition of many fields in the
$n=4$ theory, as would be required at higher order in $\kappa$, is
undesirable.  Since dynamical stability in the $n=4$ model is linked
to the preservation of spacetime interval sign, the
positive-semidefinite conformal factors required in models with $n\geq
7$ motivate a complete stability analysis of the full phylogeny.

An attractive feature for all models, independent of $n$, is that
consistency of the branching process requires that offspring universes
begin at a spatial point.  This ``big bang'' condition, though highly
dynamical from the perspective of the offspring, may appear initially
static from the perspective of the parent.  The behaviour of static,
spherically symmetric, solutions in the $N=2$ case from the parent's
perspective provides an attractive target for future investigations.

PPN constraint is a more demanding calculation.  Present frameworks
for PPN of multiscalar-tensor theories with potential are highly
non-trivial\cite{HohmannErik}, as is the PPN of theories with multiple
sources of stress (e.g. \cite{HohmannPPN}).  Given the highly
constrained nature of the model presented, encouraging cosmological
and spherically symmetric results should precede any PPN
investigation.  A reasonable first step would depart from the
brute-force perturbation pioneered by Will and Nordvedt\cite{TEGP}
toward a generalisation of Esposito-Farese's PPN technique to multiple
matter sources.  Inclusion of a potential would only be warranted if
cosmological and spherical investigations gave mandate.

It is very interesting that the model naturally features both
collisional and collisionless dark matter.  Collisional dark matter
(e.g. \cite{mahdavi2007dark}) would be inferred from dynamics of
source distributions within a single foreign contribution.
Collisionless dark matter (e.g. \cite{clowe2006direct}) would be
inferred from the dynamics of source distributions from multiple
foreign contributions.  These sources will not, in general, follow
geodesics of the native metric.  They will, however, alter native null
geodesics (i.e. lens) in the usual way.  The model also naturally
``tracks'' the native source.  Though at order $\kappa$ this tracking
is exact and hides the foreign source, at higher order this may no
longer be the case.  It should be noted that such ``tracking of the
luminous mass'' was an oft cited virtue of MOND-like phenomenologies
(e.g. \cite{milgrom1983modification, bekenstein2004relativistic}).
Cosmological constraints on Dark Matter, however, are severe; the
ratio of baryonic density to non-baryonic density cannot change.  This
would seem to exclude dominant parent and offspring contributions to
the Dark Matter density.  Sibling contributions at order $\kappa^2$,
however, provide an attractive candidate.  This is because the
Copernican Principle suggests that sibling growth and development
closely track our own.  Clearly, further exploration is justified.

\section{Conclusions}
We have implemented a Smolin-like branching multiverse as a
multiscalar extension to GR.  Our implementation seeks to produce a
more comprehensive model and remedy the lack of both population and
interaction within Smolin's original proposal.  We take the minimum
viable $N$ metric theory, in which metric degrees of freedom are
entirely decoupled, and replace $N-1$ metrics with scalar fields via
conformal relations.  The conformal relations enforce a directed,
acyclic graph structure upon the population of universes, i.e. a tree.
The result is a classical multiscalar-tensor field theory, and thus
amenable to experimental confrontation.  We analyse the model in $n$
spacetime dimensions and focus on a nearest-neighbour approximation
with $n=4$.  We determine the conditions for dynamical stability of
the this model and compute equations of motion.  The scalar equation
of motion exhibits a novel screening property: the field actively
seeks to decouple from stress under certain conditions.  We detail how
to consistently transition between an $N-1$ scalar field model to an
$N$ scalar field model, as would be required to guarantee well-defined
dynamics during reproduction events.  We compute the Newtonian limit
and show that, when applicable, the screening property reproduces
exactly Newton's equations.

\begin{acknowledgments}
The author would like to warmly thank Joel Weiner for guidance,
encouragement, and the crucial suggestion of conformal relation to
preserve causality.  The author further thanks Manuel Hohmann for
constructively critical discussions concerning early versions of this
work.  Many algebraic manipulations were verified using GNU Maxima.
Portions of this work were performed at the University of Tartu as a
Fulbright Fellow under the generous hospitality of the Laboratory of
Theoretical Physics, with financial support from the Fulbright
U.S. Student Program and the University of Hawai`i.
\end{acknowledgments}

\begin{appendix}
\section{Foil from absolute permutation symmetry}
\label{sec:foil_model}
In this appendix, we briefly discuss the following gravitational and
matter actions
\begin{align}
\delta \xo{p}{S}_M &\equiv -{1\over 2}\int\rmd^4\xi~ \sum_{q}^N  \kappa^{d(p,q)}\xo{q}{T}_{\mu\nu}\delta\xo{q}{g^{\mu\nu}}\sqrt{q} \\
\xo{p}{S}_G &\equiv \frac{1}{16\pi G}\int\rmd^4\xi~ \sum_{q}^N  \kappa^{d(p,q)} \xo{q}{R}\sqrt{q}
\end{align}
where signed relative graph depth \emph{replaces} the relative graph
distance $r(p,q)$.  Even with conformal constraint, the resultant
model is independent of the choice of $p$.  This can be seen by
noticing that a model anchored at $p$ and one anchored at any other
vertex are related by a scaling of $\kappa$ to some power.  This power
can be absorbed into the units of Newton's constant.

This absolute permutation symmetry might initially seem attractive.
We have avoided it for the following reasons.  Due to the symmetry,
the model can be analysed by anchoring within a leaf of the tree.
Then all couplings become $\kappa$ to some negative power, so let
$\kappa > 1$.  Unfortunately, all leaves enter at unit strength and
the notion of ``sibling'' is destroyed.  Similar relative couplings
across the tree grossly violate our intuitive notion of ``ancestry.''
For example, a newborn child in Australia should not influence a
similarly aged child born in Iceland more than either of their
parents, especially if this child were born two hundred years later.
Since the microevolutionary process works ``locally'' within ancestral
communities, the rigidly permutation symmetric model fails to capture
the essential aspects of actual biological populations.

\section{Conformal relation with the reals}
\label{sec:alternative}
In this appendix, we describe the essential differences between the
$n=4$ and $n \geq 7$ models.  We then briefly present the field
equations and the branching constraint in $n$ dimensions.

The essentially changed features are:
\begin{itemize}
\item{Spacetime dimension $\geq 7$: the actual dimension of spacetime
  must exceed 6.  This is not a constraint on the \emph{active}
  spacetime dimensions, and the usual compactification evasions may be
  employed.}
\item{Foreign time coordinates are always well-defined.}
\item{The partial order of the ancestral tree becomes explicitly
  physical.  An ancestor clock always runs faster, while a child clock
  always runs slower.  This can be used to define a conceptually
  distinct ``arrow of time.''}
\item{Perturbative treatment in $\kappa$ not guaranteed: the theory
  splits into two regimes, an Einstein regime similar to the pure
  phase model, and a ``nascent'' regime where the scalar field
  dominates the curvature.}
\item{Classical instabilities cannot be investigated perturbatively:
  since the perturbative expansion in $\kappa$ is no longer always
  valid, one must work to higher order in $\kappa$ to investigate
  whether terms remain positive.}
\end{itemize}
We proceed with little discussion.  The mixed metric equation of
motion is
\begin{align}
\left(1 + \kappa\sum_q\xo{q}{\exp}(\phi\sigma)\right)G^\mu_\nu + \kappa\sum_q\xo{q}{\sigma}^2\xo{q}{\exp}(\phi\sigma)\left[\delta^\mu_\nu\left(\frac{n^2 - 7n + 14}{8}\right)\nabla^\alpha\xo{q}{\phi}\nabla_\alpha\xo{q}{\phi} - \frac{(n-2)(n-5)}{4}\nabla^\mu\xo{q}{\phi}\nabla_\nu\xo{q}{\phi}\right] \nonumber \\
+\kappa\sum_q\xo{q}{\sigma}\xo{q}{\exp}(\phi\sigma)\left[\delta^\mu_\nu\nabla^\alpha\nabla_\alpha\xo{q}{\phi} - \nabla^\mu\nabla_\nu\xo{q}{\phi}\right] + \kappa\sum_q\xo{q}{\exp}(2\phi\sigma)\left[V(\xo{q}{\phi})\delta^\mu_\nu - 8\pi G\xo{q}{T^\mu_\nu}\right] - 8\pi G \xo{p}{T^\mu_\nu} + O(\kappa^2) = 0
\end{align}
while the scalar equation of motion is
\begin{align}
\xo{j}{\exp}(\sigma\phi)&\left[\xo{j}{\sigma}\xo{p}{R} + \frac{(n-1)(n-6)}{4}\left(2\xo{j}{\sigma^2}\nabla^\alpha\nabla_\alpha\xo{j}{\phi} + \xo{j}{\sigma^3}\nabla^\alpha\xo{j}{\phi}\nabla_\alpha\xo{j}{\phi}\right)\right] = \nonumber \\
&2\left[\frac{\xo{j}{\partial V}}{\partial \phi}\delta^p_j + (1-\delta^p_j) \xo{j}{\exp}(2\sigma\phi)\left(2\xo{j}{\sigma} \xo{j}{V} + \frac{\xo{j}{\partial V}}{\partial \phi}\right)\right] - \xo{j}{\exp}(2\sigma\phi)\left[8\pi G\xo{j}{T}\xo{j}{\sigma}\right] + O(\kappa) = 0
\end{align}
The branching conservation constraint in $n$-dimensions takes the following form
\begin{align}
\sigma\nabla_\mu\phi\Lambda^\mu_\nu + \nabla_\mu\Lambda^\mu_\nu =\Bigg|_{t_*} \exp(\phi \sigma)\left[8\pi G \xo{p}{\nabla}_\mu\xo{q}{T^\mu_\nu} - \nabla_\nu V + 2\sigma\nabla_\mu\phi\left(8\pi G\xo{q}{T^\mu_\nu} - V\delta^\mu_\nu\right)\right]
\label{eqn:branching_constraint}
\end{align}
where
\begin{align}
\Lambda^\mu_\nu \equiv G^\mu_\nu + \sigma^2\delta^\mu_\nu \left(\frac{n^2-7n+14}{8}\right)\nabla^\alpha\phi\nabla_\alpha\phi - \frac{(n-2)(n-5)}{4}\sigma^2\nabla^\mu\phi\nabla_\nu\phi + \sigma(\delta^\mu_\nu\nabla^\alpha\nabla_\alpha\phi - \nabla^\mu\nabla_\nu\phi).
\label{eqn:Lambda}
\end{align}
\end{appendix}