\section{Introduction}\label{sec1}

There is no doubt that Roger Penrose is one of the precursors of the mathematical foundations of causal  theory in General Relativity. His famous notes on causality \cite{penrose0} have become a classic textbook and a solid starting point for graduate students and researchers alike.  In this regard, one of Penrose's great achievements in Mathematical Relativity is the formal study of the asymptotic structure of spacetime, which arose with the introduction of the conformal compactification and Penrose diagrams \cite{zero}. It is through the notion of conformal infinity that the notion of a black hole is abstracted, thus allowing the development of  geometric and causal theoretic methods in their analysis. In spite of its success and widespread use, the conformal compactification approach to the asymptotic structure of spacetime has some downsides. Most notably, given a  spacetime, there is not a straightforward way to decide whether it admits a conformal compactification or not, or even in the affirmative case, a canonical way of constructing such a compactification. This issue was also tackled by Penrose. In their seminal work,  Geroch, Kronheimer and Penrose \cite{gkp} provided an alternative way to deal with the structure at infinity of a spacetime that relies exclusively in the causal structure of a distinguishing spacetime. Their construction follows in spirit the classical constructions in elementary geometry, where an ideal point (or point at infinity) is attached to a family of curves having a common end. Instead of using parallel rays (as is the case say, in hyperbolic geometry) they considered the causal structure of spacetime and declared that two curves have a common point at future (past) infinity if their corresponding chronological past (future) sets agree. Though elegant and simple at first sight,  the construction of the so called \emph{causal boundary} of a spacetime involves many subtleties, specially when trying to endow it with some additional structure, like a topology or causal relations (see \cite{flores:final,waveS} and references there in for an up to date account). Recently, significant results have been accomplished in this regard, most notably the introduction of a notion of black hole based on the causal boundary \cite{bhcausal}.

On the other hand, we have witnessed in the past few years a surge in the use of non-smooth geometric methods in Mathematical Relativity. Diverse settings as cone structures \cite{samann03,mingucone}, $C^0$ metrics \cite{galloling,grantc0, hevelin02,ling01,sbierski} and Lorentzian length spaces --to mention just a few-- have proven useful in exploring scenarios where (metric) smoothness can not be guaranteed, as the ones linked to recent observations \cite{ebh,ligo}. As a matter of fact, the use of non-smooth methods is not new. In the context of causality, Penrose and Kronheimer were the first to provide an abstract framework that does not require a metric structure at all \cite{kpcs}. Their notion of \emph{causal space} lays at the foundations of the theory of Lorentzian pre-length spaces first introduced by Kunzinger and S\"amman \cite{KSlls}. The purpose of this work is to present the future (or past) causal completion of a globally hyperbolic spacetime as a Lorentzian pre-length space, thus adding an interesting source of examples to this rapidly growing field \cite{samann02,hevelin,beran02,beran01,solismontes,hendike,kunzinger02}.

This work is organized as follows. In section \ref{sec:prelim}  we establish the basic facts about the causal completion and Lorentzian pre-lenth spaces, as well as the notation that will be used throughout this work. In section \ref{sec:gen} we prove that the future (past) causal completion admits a natural Lorentzian pre-length structure. Finally, in section \ref{sec:warped} we exhibit this structure in a class of warped product spacetimes.

\section{Preliminaries}\label{sec:prelim}

\subsection{Causal completions}

Let $(M,g)$ be a strongly causal spacetime and $\ll$, $\leq$ its usual chronological and causal relations, that is,  $p \ll q$ if and only if there exists a smooth future-directed timelike curve that joins $p$ with $q$, whereas $p \leq q$ if and only if there exists a smooth future-directed causal curve between these points. We define the \emph{chronological} (\emph{causal}) \emph{future} and \emph{past sets} in the standard way:
\begin{eqnarray*}
I^{+}(p) =\{q \in M \mid p \ll q\}, &\qquad& J^{+}(p)=\{q \in M \mid p \leq q\},\\
I^{-}(p) =\{q \in M \mid q \ll p\}, &\qquad& J^{-}(p)=\{q \in M \mid q \leq p\}.
\end{eqnarray*}
A sequence of points $\{x_{n}\}$ is called a \emph{future-directed chain} if $x_{n} \ll x_{n+1}$ for all $n$ and \emph{past directed} if $x_{n+1} \ll x_{n}$ for all $n \in \mathbb{N}$. Moreover, it will be called \emph{inextensible} if $\{x_{n}\}$ is not convergent. A subset $P \subseteq M$ is called \emph{past set} if $P=I^{-}(P)$ and for a subset $S \subset M$ the \emph{common past} is defined as
\[
\downarrow S=I^{-}(\{q \in M \mid q \ll s \ \forall s \in S\}).
\]
We can define the notions of \emph{future sets} and the \emph{common future} $\uparrow S$ in a time dual way. We will say that $P$ is an \emph{indecomposable past set} (or IP for short)  if $P$ cannot be written as $P=P_{1} \cup P_{2}$ where $P_{1}, P_{2} \subset P$ are disjoint proper past subsets of $P$.  As it turns out, there are only two classes of indecomposable past sets: $(1)$ the chronological past of points $I^{-}(p)$, which will be called \emph{proper idecomposable past sets}  (PIP) and $(2)$ the chronological past of inextensible future-directed chains $I^{-}(\{x_{n}\})$ which will be called \emph{terminal idecomposable past sets} (TIP).

The \emph{future causal completion} $\hat{M}$ of $(M,g)$ is the set of all idecomposable past sets (IPs). Observe that being $(M,g)$ strongly causal, then it is past distinguishing. Hence if $I^-(p)=I^-(q)$, then $p=q$. Thus, we have that any point $p \in M$ determines a unique PIP. The \emph{future causal boundary} $\hat{\partial}M$ is then identified with the TIPs. Thus we have

\begin{equation*}
	\begin{split}
IPs \equiv & PIPs \cup TIPs, \\
\hat{M}\equiv & M \cup \hat{\partial}M.
\end{split}
\end{equation*}

In a similar way we can define the \emph{past causal completion} as $\check{M} \equiv IFs$ and the \emph{past causal boundary} as  $\check{\partial}M=TIFs$. That is,

\begin{equation*}
	\begin{split}
		IFs \equiv & PIFs \cup TIFs \\
		\check{M}\equiv & M \cup \check{\partial}M
	\end{split}
\end{equation*}

Since both $\hat{M}$ and $\check{M}$ include a copy of the spacetime $(M,g)$, it is natural trying to define the total causal completion  of $(M,g)$ as $\hat{M}\cup\check{M}$, where the PIP $I^-(p)$ is identified with the PIF $I^+(p)$ ---since both represent the point $p\in M$. However, it was observed early on that such a construction may lead to inconsistencies, as in some cases further identifications on the boundaries $\hat{\partial}M$ and $\check{\partial}M$ ought to take place \cite{gkp,marolf,waveS,Sz}.  Thus finding an approach to constructing the causal boundary that works in full generality proved to be a delicate task, which was completed only recently \cite{flores:final}.   In order to avoid such intricacies, we will focus only on the future causal completion $\hat{M}$.

Notice that, as the notation suggest, the  relation $\hat{\ll}$ on $\hat{M}$ given by
\begin{equation}
	\label{eq2}
	\begin{split}
	P \hat{\ll} Q &\Leftrightarrow \exists q \in Q\setminus P {\text{ such that }} P \subset I^{-}(q)  \\
	\end{split}
\end{equation}
is indeed transitive. In fact, if $P \hat{\ll} Q \hat{\ll} R$ there exist $q \in Q \setminus P$ and $r \in R \setminus Q$ such that $P \subset I^{-} (q) \subset Q \subset I^{-}(r) \subset R$, thus $P \hat{\ll} R$.

Hence, we can think of $\hat{\ll}$ as providing a chronological structure on $\hat{M}$ \cite{harrischro}. The chronological future and past sets so induced will be denoted by $\hat{I}^+(P)$ and $\hat{I}^-(P)$.

Further, we can endow $\hat{M}$ with a sequential topology which is compatible with the chronology $\hat{\ll}$ just defined.\footnote{Recall that on any spacetime $(M,g)$ the chronological sets $I^+(p)$, $I^-(p)$ are open in the manifold topology.} Thus, consider the  limit operator $\hat{L}_{chr}$ over the sequences of past sets given by
\begin{equation}
\label{eq3}
P \in \hat{L}_{chr}(\{P_{n}\}\}) \Leftrightarrow P \subset LI(\{P_{n}\}) \text{ and it is maximal in } LS(\{P_{n}\}).\footnote{The theoretical inferior and superior limits of subsets are defined as $LI(\{P_{n}\})=\bigcup_{n=1}^{\infty}\bigcap_{m=n}^{\infty} P_{m}$ and $LS(P_{n})=\bigcap_{n=1}^{\infty} \bigcup_{m=n}^{\infty} P_{m} $, respectively.}
\end{equation}
and define the \emph{future chronological topology} $\hat{\mathcal{T}}_{chr}$ by its closed subsets as follows: a subset $C \subset \hat{M}$ is closed if and only if for any sequence $\sigma \subset C$ the contention $\hat{L}_{chr}(\sigma) \subset C$ holds.  Thus we have the following result (see the proof of Thrm. 3.27 in \cite{flores:final})
\begin{theorem}
\label{teo:propcomple}
Let $(M,g)$ be a strongly causal spacetime and $\hat{M}$ its future causal completion endowed with the chronological structure induced by \eqref{eq2} and the topology induced from the chronological limit \eqref{eq3}. Then:
	\begin{itemize}
		\item[(i)] The inclusion $M \hookrightarrow \hat{M}$ is continuous. Moreover, the restriction of the chronological topology to $M$ is the manifold topology
		\item[(ii)]  The future causal completion is complete: any future directed chain $\{P_n\}$ in $\hat{M}$ converges in $\hat{\mathcal{T}}_{chr}$. In particular, any future inextensible timelike curve $\gamma$ on $M$
		has an endpoint in $\hat{M}$.

		\item[(iii)] The sets $\hat{I}^+(P)$ and  $\hat{I}^-(P)$ are open for all $P \in \hat{M}$.

		\item[(iv)] $(\hat{M},\hat{\mathcal{T}}_{chr})$ is a $T_1$ topological space.

	\end{itemize}
\end{theorem}

Since there are examples in which $(\hat{M},\hat{\mathcal{T}}_{chr})$ is not Hausdorff, we can not aim to furnish a manifold structure on it for such cases. However, some extra structure can be added in particular cases. For instance, a linear connection can be defined on $\hat{M}$ provided that $(M,g)$ is static and spherically symmetric  \cite{harriscon}.

\subsection{Lorentzian pre-length spaces}

In their remarkable paper \cite{kpcs}, Kronheimer and Penrose developed the notion of a causal space by abstracting the fundamental properties of the chronological and causal relations. Such an axiomatic approach to causality has proven useful in many circumstances, for instance in the causet approach to quantum gravity \cite{surya}. In recent times, the search for applying synthetic geometrical methods to mathematical relativity sparkled a renewed interest in developing causality in an abstract setting. Lorentizan pre-length spaces are a refinement of the notion of causal spaces that incorporates a time distance function, thus providing an analog for the length structure that serves as a building block for the well established synthetic theory in metric spaces \cite{burago}.

As in the seminal work by Kunzinger and S\"amann \cite{KSlls} we define a Lorentzian pre-length space  $(X,d,\ll,\leq, \tau)$ as a metric space $(X,d)$ along with two relations $\ll$, $\leq$,  ---named \emph{chronological} and \emph{causal}--- and a \emph{time separation function} $\tau: X\times X\to [0,\infty]$ that satisfies the following axioms:
\begin{enumerate}
\item $\leq$ is a pre-order,
\item  $\ll$ is a transitive relation contained in $\leq$,
\item $\tau$ is a lower semi-continuous function ---with respect to $d$--- satisfying
\begin{itemize}
    \item $\tau(x,z)\geq \tau(x,y) + \tau(y,z)$ for all $x\leq y \leq z$,
    \item $\tau(x,y)>0$ if and only if $x\ll y$.
\end{itemize}
\end{enumerate}

 As an immediate consequence of the definition we have two of the most important features of the causal structure of a spacetime: (i) the chronological sets $I^+(p)$, $I^-(p)$ are open, and (ii) either $x\leq y\ll z$ or $x\ll y\leq z$ implies $x\ll z$.\footnote{This is commonly known as the \emph{push up property}.}  These properties, along some additional structure enable us to build a causality ---an even a causal hierarchy--- that closely resembles the usual causal structure of a smooth spacetime \cite{ACS,KSlls}.  Examples of Lorentzian pre-length spaces include a wide variety of structures, like cones \cite{samann03,mingucone}, causally plain $C^0$ spacetimes \cite{grantc0}, contact structures \cite{hendike} and Lorentzian taxicab-type spaces \cite{solismontes}.

\section{Lorentzian pre-length space structure on $\hat{M}$}\label{sec:gen}

As we discuss in the previous section, the topological space $(\hat{M},\hat{\mathcal{T}}_{chr})$ might not be metrizable in general. Thus, in order to construct a Lorentzian pre-length space structure on it, a refinement on the topology is in order.  In a recent work \cite{costaflores:newtopology}, Costa, Flores and Herrera studied the so called \emph{closed limit topology} $\hat{\mathcal{T}}_{c}$ (or CLT for short). The Hausdorff limit operator
\begin{equation}
\hat{L}_{H}(\{P_{n}\})=\{P \in \hat{M} \mid P=\liminf(\{P_n\})=\limsup(\{P_{n}\})\}
\end{equation}
generates the closed sets in $\hat{\mathcal{T}}_{c}$ as follows: a subset $C \subset \hat{M}$ is closed if and only if $\hat{L}_{H}(\{P_{n}\}) \subset C$ for every sequence $\{P_{n}\} \subset C$.\footnote{Observe that  $\# L_{H}(\sigma) \in \{0,1\}$ and thus $\mathcal{T}_{c}$ is a {\it first order topology} (see \cite{covering}). In other words, the following equivalence holds
\[
\{P_{n}\} \rightarrow_{\hat{\mathcal{T}}_{c}} P \Leftrightarrow \hat{L}_{H}(\{P_{n}\})=\{P\}.
\] }
We now summarize the relevant features of the CLT topology in globally hyperbolic spacetimes (refer to Thrms 4.1 and 4.2 in \cite{costaflores:newtopology})
\begin{theorem}\label{teo:clt}
If $(M,g)$ is globally hyperbolic, then, the following statements hold for the topological space $(\hat{M},\hat{\mathcal{T}}_{c})$:
\begin{itemize}
	\item[$(i)$] The natural inclusion $i: M \rightarrow \hat{M}$ given by $i(p)=I^{-}(p)$ is an open continuous map. Moreover, $i(M)$ is an open dense subset of $\hat{M}$, the induced topology on $M$ is the manifold topology, $\hat{\partial} M$ is closed and $(\hat{M},\hat{\mathcal{T}}_{c})$ is second countable.

	\item[$(ii)$] The chronological sets $\hat{I}^{\pm}(P)$ are open subsets for all $P \in \hat{M}$.

	\item[$(iii)$] Any future directed chain in $\{P_{n}\} \subset \hat{M}$ converges in $\hat{\mathcal{T}}_{c}$
		\item[$(iv)$] The topological space $(\hat{M},\hat{\mathcal{T}}_{c})$ is metrizable.
\end{itemize}
\end{theorem}

Notice that in general, the chronological topology is coarser than $CLT$. Moreover, they coincide when their corresponding limit operators agree. The next result provides necessary and sufficient conditions for that (see \cite[Thrm. 5.3]{costaflores:newtopology})

\begin{theorem}\label{teo:opagree}
Let $\hat{M}$ be a future completion endowed with the limit operators $L_{H}$ and $\hat{L}_{chr}$.
Both limit operators coincide if, and only if, the chronological topology is Hausdorff.
\end{theorem}

Recall $\hat{M}$ is endowed with the chronological relation $\hat{\ll}$ given by \eqref{eq2}. Thus, $\hat{\ll}$ naturally induces an associated causal relation ${\leq}_\circ$ in $\hat{M}$  (see \cite[Definition 2.22]{sanchezminguzzi}) by setting
\begin{equation*}
P {\leq}_\circ Q \qquad \Leftrightarrow \qquad \hat{I}^{-} (P) \subset \hat{I}^{-} (Q) \quad \text{and} \quad \hat{I}^{+} (Q) \subset \hat{I}^{+} (P).
\end{equation*}
By construction, this relation is reflexive, transitive and contains $\hat{\ll}$, thus it is a causal relation according to the definition of a pre-length space.

The causal relation $\leq_\circ$ though constructed by a standard procedure, it may seem artificial at first sight. The simple contention of past sets provides right away a pre-order in $\hat{M}$. That is,  the relation $\hat{\leq}$ on $\hat{M}$ given by
$$P \hat{\leq} Q \Leftrightarrow P \subset Q.$$
might be considered as well \cite{outline}.

The following well known result will be used extensively from now on. We include its proof here for completeness (see for example \cite[Proposition 2.11]{costaflores:newtopology})

\begin{proposition}\label{lemadesigualdad}
Let $(M,g)$ is a globally hyperbolic spacetime. If $P \in \hat{\partial} M$ then for any $p\in M$ we have $P\hat{\nleq} I^-(p)$.  \footnote{Equivalently, if $P,Q \in \hat{M}$ are IPs such that $P \in \hat{\partial} M$ and $P \hat{\leq} Q$, then $Q \in \hat{\partial} M$. }
\end{proposition}

\begin{proof}
Suppose by contradiction that $P \subset I^{-}(p)$ and consider $\gamma:[a,b) \rightarrow M$ an inextensible future directed timelike curve such that $P=I^{-}(\gamma)$. Since $P$ is a TIP
then we must have that $P$ is a proper subset of $I^{-}(p)$, thus, if we take $\gamma(a)$ then we must have that $\gamma(s) \in J^{+}(\gamma(a)) \cap I^{-}(p) \subset J^{+}(\gamma(a)) \cap J^{-}(p)$
for all $s \in (a,b)$ and this gives an inextensible timelike curve imprisoned in a compact subset thus contradicting strong causality.
\end{proof}

We now show that the causal relations $\leq_\circ$ and $\hat{\leq}$ coincide when $(M,g)$ is globally hyperbolic.

\begin{proposition}\label{relacionescausalescoinciden}
Let $(M,g)$ be a globally hyperbolic spacetime. For all $P,Q \in \hat{M}$, $P\hat{\leq} Q $ if and only if $P {\leq}_\circ Q$.
\end{proposition}

\begin{proof}
Suppose first that $P \hat{\leq} Q$, i.e. $P \subset Q$. If $Q$ is a TIP then by Proposition \ref{lemadesigualdad} $\hat{I}^{+} (Q) = \emptyset$ and we have $\hat{I}^{+} (Q) \subset \hat{I}^{+} (P)$. If $Q = I^{-} (q)$ is a PIP, consider $R  \in \hat{I}^{+} (Q)$, then $Q \hat{\ll} R$, which means there exists $r \in R \setminus Q$ such that $Q \subset I^{-} (r)$ and since $ P \subset Q \subset I^{-} (r)$ with $r \in R \setminus P$ we have $R \in \hat{I}^{+} (P)$. Thus $\hat{I}^{+} (Q) \subset \hat{I}^{+} (P)$.

Similarly, if $R \in \hat{I}^{-} (P)$ then $R \hat{\ll} P$, and there exists $r \in P \setminus R$ with $R \subset I^{-}(p)$. Since $P \subset Q$, $p$ also belongs to $Q$ and we get $R \hat{\ll} Q$, that is, $R \in \hat{I}^{-} (Q)$. This means that $\hat{I}^{-} (P) \subset \hat{I}^{-} (Q)$ and we have $P {\leq}_\circ Q$.

On the other hand, let us suppose $P {\leq}_\circ Q$, that is,
$$\hat{I}^{+} (Q) \subset \hat{I}^{+} (P) \qquad \text{ and } \qquad \hat{I}^{-} (P) \subset \hat{I}^{-} (Q).$$
Take $p_0 \in P$, then there exists $q_0 \in P \setminus I^{-} (p_0)$ with $I^{-} (p_0) \subset I^{-} (q_0)$ therefore $I^{-} (p_0) \hat{\ll} P$, that is, $I^{-} (p_0) \in \hat{I}^{-}(P) \subset \hat{I}^{-}(Q)$. Thus, there exists $q \in Q\setminus I^{-}(p_0)$ con $I^{-} (p_0) \subset I^{-} (q) \subset Q$. If $\{p_n\}$ is a future directed chain generating  $I^{-}(p_0)$ we have that $p_n \in I^{-}(q)$ for $n$ large enough and then $p_0 \in \overline{I^{-}(q)} = J^{-}(q)$, since $M$ is causally simple. In consequence, $p_0 \leq q \ll q_m,$ for some $m$ large enough, where $\{q_m\}$ is a future-directed chain that generates $Q = I^{-}(\{q_m\})$. Then, $p_0 \in Q$ and $P \subset Q$, which by definition means $P \hat{\leq} Q$.
\end{proof}

As an immediate consequence, we have the following

\begin{coro}\label{hatMescausal}
Let $(M,g)$ be a globally hyperbolic spacetime. Then the relation $\hat{\leq}$ is a partial order on $\hat{M}$.\footnote{This amounts to saying that the causal space $(\hat{M},\hat{\ll},\hat{\leq})$ satisfies the causality axiom. }
\end{coro}

As expected, the relations $\hat{\ll}$ and $\hat{\leq}$ naturally extend those of $M$.

\begin{proposition}\label{relacionesextendidas}
Let $(M,g)$ be globally hyperbolic. For all $p,q \in M$, $I^{-} (p) \hat{\ll} I^{-} (q)$ if and only if $p \ll q$. Moreover, $I^{-} (p) \hat{\leq} I^{-} (q)$ if and only if $p \leq q$.
\end{proposition}

\begin{proof}
Assume first that $I^{-} (p) \hat{\ll} I^{-}(q)$. This implies, by definition, that there is a $r \in I^{-}(q) \setminus I^{-} (p)$ such that $I^{-} (p) \subset I^{-} (r)$. Then $p \ll r \ll q$ which implies $p \ll q$.

Conversely, if $p \ll q$ then $p \neq q$ because $(M,g)$ is distinguishing. Consider $r \in M$ with $p \ll r \ll q$. Thus, $r \in I^{-}(q) \setminus I^{-} (p)$ with $I^{-} (p) \subset I^{-} (r)$, that is, $I^{-} (p) \hat{\ll} I^{-} (q)$. This concludes the first part of the proof.

Now let us suppose that $I^{-} (p) \hat{\leq} I^{-} (q)$, that is, $I^{-} (p) \subset I^{-} (q)$. For $p \in M$, we take $\{p_n\}$ a future-directed chain that generates $I^{-} (p)$, that is, a sequence with $p_n \ll p_{n+1}$ for all $n \in \mathbb{N}$ and $p_n \to p$. Notice that for all $n \in \mathbb{N}$ we have that $p_n \ll p$, that is, $p_n \in I^{-} (p)$ and so $p_n \in I^{-}(q)$. Therefore $p_n \to p \in \overline{I^{-}(q)} = J^{-}(q)$ because $(M,g)$ is causally simple. Therefore, $p \leq q$.

On the other hand, assume $p \leq q$. It suffices to prove that $I^{-} (p) \subset I^{-} (q)$. If $x \ll p$, by the push-up property we get $x \ll q$. Therefore, $I^{-} (p) \subset I^{-}(q)$.
\end{proof}

So far, we have exhibited a metric topology, as well as chronological and causal relations in the future causal completion $\hat{M}$. In the remaining of this section we will be dealing with the construction of an adequate time separation $\hat{\tau}$  for $\hat{M}$. As a first step, notice that if $Q \in \hat{M}$ is generated by a future-directed chain $\{q_n\}$  and $p \in Q$, then there exists $N \in \mathbb{N}$ such that for all $n\geq N$, $p \ll q_n \ll q_{n+1}$. By the reverse triangle inequality of $\tau$ we have
$$\tau (p, q_{n+1}) \geq \tau (p,q_n) + \tau (q_n, q_{n+1}) > \tau (p,q_n).$$
Hence, the sequence of real numbers $\{\tau (p,q_n) : p \ll q_n\}$ is strictly increasing. Consequently, in view of Proposition \ref{lemadesigualdad} we define $\hat{\tau}$ as follows: for $p \in M$, $P,Q \in \hat{M}$
\begin{itemize}
\item[(i)] $\hat{\tau} (P,Q) := 0$ if $P \in \hat{\partial} M$,
\item[(ii)] $\hat{\tau} (I^{-} (p), Q) := \sup \{\tau (p,q_n): Q = I^{-} (\{ q_n \})\}.$
\end{itemize}

Since the chain $\{q_n\}$ that generates a past set $Q$ might not be unique, we first have to show that $\hat{\tau}$ is well defined. First note that if $\hat{\tau} (P,Q) = 0$ there is nothing to prove.

Alternatively, If $\{q_n\}$ and $\{r_m\}$ are two different future-directed chains generating $Q$. It is enough to show that
$$\sup \{\tau(p,q_n) : Q = I^{-} (\{q_n\})\} = \sup \{\tau(p,r_m) : Q = I^{-} (\{r_m\})\}.$$
For any $p \in Q$ there exists $q_n \in Q$ with $p \ll q_n$, and hence there exists $r_m \in Q$ with $q_n \ll r_m$. Thus $p \ll q_n \ll r_m$ and by the reverse triangle inequality
$$\sup \{ \tau (p,r_m) \} \geq \tau (p,r_m) \geq \tau (p,q_n) + \tau (q_n, r_m) > \tau (p,q_n).$$
Hence $\sup \{ \tau (p,r_m) \}$ is an upper bound for $\{\tau (p, q_n)\}$ and thus
$$\sup \{\tau(p,q_n) : Q = I^{-} (\{q_n\})\} \leq \sup \{\tau(p,r_m) : Q = I^{-} (\{r_m\})\}.$$
Similarly,
$$\sup \{\tau(p,r_m) : Q = I^{-} (\{r_m\})\} \leq \sup \{\tau(p,q_n) : Q = I^{-} (\{q_n\})\}.$$

We now notice that $\hat{\tau}$ extends $\tau$ to the future causal completion $\hat{M}$.

\begin{proposition}\label{htauextiende}
If $P=I^-(p)$ and $Q=I^-(q)$, then $\hat{\tau} (P,Q) = \tau (p,q)$.
\end{proposition}

\begin{proof}
By definition
$$\hat{\tau} (P,Q) = \sup \{\tau (p,q_n) : Q = I^{-} (\{q_n\})\}.$$
Since $Q = I^{-} (q)$ is a PIP, the chain $\{q_n\}$ must converge to $q$. Moreover, the sequence $\{\tau (p,q_n)\}$ is non decreasing and since $(M,g)$ is globally hyperbolic, $\tau$ is a continuous function. Therefore
$$\hat{\tau} (P,Q) = \sup \{ \tau (p,q_n): Q = I^{-}(\{q_n\}) \} = \lim_{n \to \infty} \tau (p,q_n) = \tau (p,q).$$
\end{proof}

The following series of lemmas are intended to show that $\hat{\tau}$ defines a time separation function on $(\hat{M},\hat{\mathcal{T}}_c)$ when $(M,g)$ is a globally hyperbolic spacetime.

\begin{lemma}[Positivity]\label{clasificacausalidad}
If $P,Q \in \hat{M}$, then $\hat{\tau} (P,Q) > 0$ if and only if $P \hat{\ll} Q$.
\end{lemma}

\begin{proof}
We begin assuming that $P \hat{\ll} Q$. Then $P = I^{-}(p)$ for some $p\in M$ by Proposition \ref{lemadesigualdad}.  By definition, there exists $q \in Q\setminus I^{-}(p)$ such that $I^{-}(p) \subset I^{-} (q)$, then $p \leq q$. Moreover, if $\{q_n\}$ is a future directed chain generating $Q$ then there exists $q_m \in \{q_n\}$ such that $p \leq q \ll q_n$ for all $n\ge m$ which implies  by the reverse triangle inequality
$$\tau (p,q_n) \geq \tau (p,q) + \tau(q,q_n) > 0.$$
Therefore
$$\hat{\tau}(I^{-}(p), Q ) = \sup \{\tau(p,q_n) :  Q = I^{-} (\{q_n\})\} > 0.$$

Conversely, assume that $\hat{\tau} (P,Q) > 0$. Which by the definition of $\hat{\tau}$ means that $P$ is a PIP. Suppose that $P = I^{-}(p)$. If $\hat{\tau} (I^{-}(p) , Q) < +\infty$, for any $\epsilon > 0$ there exists $q_n$ with
$$\tau (p,q_n) > \hat{\tau}(I^{-}(p),Q) - \epsilon.$$
If we take $\epsilon = \frac{\hat{\tau}(I^{-}(p),Q)}{2}$, we have that $\tau(p,q_n) > 0$ which implies $p \ll q_n \ll q_{n+1}$. Hence, $q_n \in Q \setminus I^{-}(p)$ with $I^{-}(p) \subset I^{-}(q_n)$ and then $I^{-}(p) \hat{\ll} Q$. Now, if $\hat{\tau}(I^{-} (p), Q) = +\infty$, for any $N \in \mathbb{N}$ there exists $q_m \in \{q_n\}$ such that
$$\tau (p,q_m) \geq N > 0,$$
which implies $p \ll q_n$ for all $n\geq m$. Then $I^{-}(p) \hat{\ll} Q$.
\end{proof}

The following result is needed in the proof of the reverse triangle inequality.

\begin{lemma}\label{lema:rev}
If $p \notin Q$ or $p \in \partial Q$, then $\hat{\tau} (I^{-} (p), Q) = 0$.
\end{lemma}

\begin{proof}
Let $\{q_n\}$ be a future-directed chain generating $Q \in \hat{M}$. Proceeding by contradiction, let us suppose that $\hat{\tau} (I^{-} (p), Q ) > 0$ then there exists $q_k \in \{ q_n \}$ with $\tau (p,q_k) > 0$, but this only occurs if $p \ll q_k$, which is a contradiction because $p \notin Q$ or $p \in \partial Q$, and in either case, there can be no such $q_k$.
\end{proof}

\begin{lemma}[Reverse triangle inequality]\label{desigualdadtrianguloinvertida}
If $P,Q,R \in \hat{M}$ are such that $P \hat{\leq} Q \hat{\leq} R$, then
$$\hat{\tau} (P,R) \geq \hat{\tau} (P,Q) + \hat{\tau} (Q,R).$$
\end{lemma}

\begin{proof}
We proceed  case by case. If $P,Q,R$ are all TIPs then
$$\hat{\tau} (P,R) = 0 \geq \hat{\tau}(P,Q) + \hat{\tau}(Q,R) = 0.$$
Moreover, if $P$ is a TIP by Proposition \ref{lemadesigualdad}, $Q$ and $R$ are TIPs as well.
Also, if $P,Q,R$ are all PIPs, that is, $P = I^{-}(p), Q= I^{-}(q)$ and $R = I^{-}(r)$ for some $p,q,r \in M$ then the reverse triangle inequality of $\tau$ gives us the result.
Thus, we have only two cases left, namely
\begin{itemize}
\item[(i)] $P \hat{\leq} Q \hat{\leq} R$ with $P=I^-(p)$, $Q=I^-(q)$, $R \in \hat{\partial} M$;
\item[(ii)] $P \hat{\leq} Q \hat{\leq} R$ with $P= I^{-} (p)$, $Q$, $R \in \hat{\partial} M$.
\end{itemize}

For case (i), notice that if $p\in\partial R$ then $q\not\in R$. Moreover, if
$p$, $q\in \partial R$, then $p\rightarrow q$\footnote{Recall that the \emph{horismos} relation $p\rightarrow q$ is defined as $p\leq q$ and $p\not\ll q$.} follows. Thus we only have the following possibilities:
\begin{itemize}
\item[(a)] $p \rightarrow q$ and $p, q \in \partial R$,
\item[(b)] $p \ll q$ with $q \in \partial R$,
\item[(c)] $p \ll q$ and $q \in R$,
\item[(d)] $p \rightarrow q$, $p \in R$ and $q \in \partial R$,
\item[(e)] $p \rightarrow q$ and $p,q\in R$.
\end{itemize}

In case (a) observe that $\tau(p,q)=0$, $\hat{\tau}(I^{-}(p),R)=\hat{\tau}(I^{-}(q),R)=0$ in virtue of Lemma \ref{lema:rev}. Thus the triangle inequality holds trivially.

Now, for case $(b)$  consider $R=I^{-}(\{r_{m}\})$. Since
 $p \in R$ and $\hat{\tau}(I^{-}(q),R)=0$,  take a future directed timelike sequence $\{q_{n}\}$ such that $q_{n} \in
I^{-}(q)$ and $q_{n} \rightarrow q$. Thus, $ p \ll q_{n} \ll r_{m_n}$ for some $n$ and $m_n$ natural numbers and this leads to
\[
\tau(p,r_{m_n}) \geq \tau(p,q_n) + \tau(q_n,r_{m_{n}}) \geq \tau(p,q_n)
\]
by the reverse triangle inequality of $\tau$. Since $(M,g)$ is globally hyperbolic, the time separation function $\tau$  is continuous, thus
$$\hat{\tau}(I^-({p)},R) \geq \hat{\tau}(I^{-}(p),I^{-}(q))+0= \hat{\tau}(I^{-}(p),I^{-}(q))+\hat{\tau}(I^{-}(q),R).$$

For (c) it is enough to prove
$$\hat{\tau} (I^{-} (p), R) \geq \tau (p,q) + \hat{\tau} (I^{-}(q) , R).$$
Let $\{r_m\}$ be a future-directed chain generating $R$. By definition
\begin{align*}
\hat{\tau} (I^{-} (p),R) :=& \sup \{\tau (p,r_m): R= I^{-} (\{r_m\})\} \\
\hat{\tau} (I^{-} (q),R) :=& \sup \{\tau (q,r_m): R= I^{-} (\{r_m\})\}
\end{align*}
We know that $I^{-}(p) \hat{\leq} I^{-} (q)$ if and only if $p \leq q$. Therefore $p\leq q \ll r_m$ for large $m$ and  we have
$$\tau (p,r_m) \geq \tau(p,q) + \tau (q,r_m),$$
which in turn implies
$$\hat{\tau} ( I^{-} (p), R) \geq \tau (p,q) + \hat{\tau} (I^{-}(q), R).$$

In case $(d)$ we have that
$\hat{\tau}(I^{-}(q),R)=0=\hat{\tau}(I^{-}(p),I^{-}(q))$ and since $p \in R$ we have that
$p \ll r_{m}$ for some $m \in \mathbb{N}$ which leads to $\tau(p,r_{m})>0$ and this gives $\hat{\tau}(I^{-}(p),R)>0=\hat{\tau}(I^{-}(p),I^-(q))+\hat{\tau}(I^{-}(q),R)$.

In order to prove $(e)$ note that $\hat{\tau}(I^{-}(p),I^{-}(q))=0$ and following the same argument in point $(c)$ we have that $p \leq q \ll r_{m}$ for  large $m$ and thus $\tau(p,r_{m}) \geq \tau(q,r_{m})$ by the inverse triangle inequality in $(M,g)$. Thus, we have $\hat{\tau}(I^{-}(p),R) \geq \hat{\tau}(I^{-}(q),R)$ which is the reverse triangle inequality in this case.

Similarly, for case (ii) notice that $R\not\in\hat{I}^+(Q)$, thus we have the following possibilities:

\begin{itemize}
	\item[(a)] $p \in \partial Q \cap \partial R$,
	\item[(b)] $p \in \partial Q$ and $p \in R$,
	\item[(c)] $p \in Q \cap R$.
\end{itemize}

In  case (a) observe that all the quantities involved are zero and thus the reverse triangle inequality holds trivially.

In case $(b)$ we have that
$\hat{\tau}(I^{-}(p),Q)=0$ and $\hat{\tau}(Q,R)=0$ and since $p \in R$ we have that for any future directed sequence $\{r_m\}$ that generates $R$ we have
$p \ll r_{m}$ for $m \geq m_{0}$. The latter implies that $\hat{\tau}(I^{-}(p),R) \geq \tau(p,r_{m})>0$ and thus reverse triangle inequality holds trivially again.

Finally, for (c) we have to prove that
$$\hat{\tau}(I^{-} (p) , R) \geq \hat{\tau} (I^{-}(p) , Q)$$
since $\hat{\tau} (Q,R) = 0$, by definition. In order to verify this, take $\{q_n\}$ and $\{r_m\}$ future-directed chains generating $Q$ and $R$, respectively.  Observe that $Q \subset R$ since they are causally related in $\hat{M}$, so there are $n, m_{n} \in \mathbb{N}$ such that
$p \ll q_{n} \ll r_{m_{n}}$. Therefore, we have $\tau (p,r_{m_{n}}) \geq \tau (p,q_n)$ for large $n$,
and thus
$\hat{\tau} (I^{-}(p) , R) \geq \hat{\tau} (I^{-}(p) , Q)$ holds.
\end{proof}

All that is left for $\hat{\tau}$ to be a time separation function is for it to be lower semicontinuous.

\begin{lemma}[Lower semicontinuity]\label{lowersemicontinua}
The function $\hat{\tau}$ is lower semicontinuous in $\hat{M}$, i.e. for all $P,Q \in \hat{M}$, and any $\delta > 0$, if $P_n$ and $Q_n$ are sequences that converge to $P$ and $Q$, respectively, then there exists $N \in \mathbb{N}$ such that
$$\hat{\tau} (P_n , Q_n) > \hat{\tau} (P,Q) - \delta$$
for all $n\geq N$.
\end{lemma}

\begin{proof}
If $\hat{\tau}(P,Q) = 0$, there is nothing to prove. Otherwise, if $\hat{\tau} (P,Q) > 0$, then $P \hat{\ll} Q$  and $P = I^{-} (p)$ for some $p\in M$ by Proposition \ref{lemadesigualdad} . Let $\delta > 0$ and $\{q_k\}$ be a future-directed chain generating $Q$, and for every $n$, let $\{q^n _k\}$ be a future-directed chain generating $Q_n$. Given that $P_n$ converges to $P$ and $\hat{\partial} M$ is closed in virtue of Theorem \ref{teo:clt}, there exists $N_1 \in \mathbb{N}$, such that for all $n \geq N_1$, $P_n$ is a PIP, that is, $P_n = I^{-}(p_n)$ with $p_n \in M$ and $p_n$ converges to $p$.

Since $p_n$ converges to $p$, there exists $N_2 \in \mathbb{N}$ such that for all $n \geq N_2$, $I^{-} (p_n) \in \hat{I}^{-}(Q)$, since the latter is an open neighborhood for $I^{-}(p)$. Thus, for all $n \geq N_2$ there exists $k$ such that $p_n \ll q_k \in Q$. Moreover, since $Q = \liminf (Q_n) = \limsup (Q_n)$, then every $q_k$ is in all but a finitely many of the sets $Q_n$. Therefore, there exists $q^{n}_{r}$ with $p_n \ll q_k \ll q^{n}_{r}$ and as a consequence, by the reverse triangle inequality $\tau (p_n, q^{n}_{r}) > \tau (p_n, q_k).$
On the other hand, using that $\tau$ is lower semicontinuous, we have that there exists $N_3 \in \mathbb{N}$ such that for all $n \geq N_3$
$$\tau (p_n, q_k) > \tau (p, q_k) - \frac{\delta}{2}.$$
Then, if $n \geq  N = \max \{N_1,N_2,N_3\}$
$$\tau (p_n, q^{n}_{r}) > \tau (p, q_k) - \frac{\delta}{2},$$
whence $\tau (p_n, q^{n}_{r}) + \frac{\delta}{2}$ is an upper bound for $\tau (p, q_k)$. Which implies
$$\tau (p_n, q^{n}_{r}) \geq \sup \{\tau (p, q_k): Q = I^{-}(\{q_k\}) \} - \frac{\delta}{2}.$$
By definition of least upper bound
$$\sup \{\tau (p_n, q^{n}_{r}) : Q_n = I^{-}(\{q^{n}_r\}) \} > \sup \{\tau (p, q_k): Q = I^{-}(\{q_k\}) \} - \delta,$$
which by definition of $\hat{\tau}$ means that for all $n \geq N$
$$\hat{\tau} (P_n , Q_n) > \hat{\tau} (P,Q) - \delta.$$
\end{proof}

The following theorem summarizes the results of the previous lemmas.

\begin{theorem}\label{teoremaprincipal}
Let $(M,g)$ be a globally hyperbolic space-time. Then $(\hat{M}, d_c, \hat{\ll}, \hat{\leq} , \hat{\tau})$ is a Lorentzian pre-length space.
\end{theorem}

It is important to note that the relations $\hat{\ll}$, $\hat{\leq}$ as well as the CLT topology are not affected by a conformal change on the spacetime metric $g$. However, this is not the case for the time separation $\tau$. Thus, it is expected to have different Lorentzian pre-length structures in the future causal completion within the same conformal class.

\section{Applications}\label{sec:warped}

In this section we show that a class of warped product spacetimes satisfies the condition of Thrm \ref{teo:opagree} and hence its associated chronological and CLT topologies coincide.  This result can be used to carry out explicit calculations regarding the pre-length structure of the future causal completion.\footnote{A thorough description of the future causal boundary of such spacetimes can be found in \cite{valana}.} Finally, as an illustrative example we consider the particular case of de Sitter spacetime.

Recall that a Lorentzian warped product is a manifold $V=\mathbb{R} \times M$ furnished with a metric of the form
$$g=-dt^{2} + \alpha(t) h,$$
where, $(M,h)$ is a Riemannian manifold and $\alpha$ is a smooth positive function over $\mathbb{R}$.

The chronological relation on these spacetimes can be characterized as follows (see \cite[Section 2.2]{joniluis}):
		\begin{lemma}
			Let $(V,g)$ be a Lorentzian warped spacetime. If $(t_{0},x_{0}), (t_{1},x_{1})$ points in $V$, then,
			\begin{equation*}
				(t_{0},x_{0}) \ll (t_{1},x_{1}) \Leftrightarrow d(x_{0},x_{1}) < \int_{t_{0}}^{t_{1}} \frac{ds}{\sqrt{\alpha(s)}}
			\end{equation*}
		\end{lemma}

 The future causal completion can be characterized depending on the value of $\int_{0}^{+\infty} \frac{ds}{\sqrt{\alpha(s)}}$. Indeed we have\footnote{ Here $\partial_\mathcal{B} M$ and $\partial_CM$ are the Busemann and (metric) Cauchy boundaries or $M$. Refer to \cite{joniluis,valana,isobound} for a detailed account on the structure and topology of $\hat{V}$.}

\begin{theorem}\label{XX}
Let $V=\mathbb{R} \times M$ with $g=-dt^2+\alpha(t) h$ be a warped spacetime and $(M,d)$ a locally compact metric space. Then,
\begin{enumerate} \setlength\itemsep{1em}
\item If $\int_{0}^{+\infty} \frac{ds}{\sqrt{\alpha(s)}}= \infty$ ,then, the future causal boundary $\hat{\partial} V$ is an infinite null cone with base $\partial_{\mathcal{B}}M \setminus \partial_{C} M$ with apex in $i^{+}$ and timelike lines over each point in $\partial_{C} M$ and final point in $i^{+}$. Moreover, $\hat{V}$ is homeomorphic to $V \cup (\mathbb{R} \times \partial_{C} M) \cup (\mathbb{R} \times \partial_{\mathcal{B}} M) \cup i^{+}$.

\item If $\int_{0}^{+\infty} \frac{ds}{\sqrt{\alpha(s)}} < \infty$, then, $\hat{\partial} V$ is a copy of the Cauchy completion $M^{C}$ of $(M,h)$ and timelike lines over each point in $\partial_{C} M$ that finishes in the same point at the copy at infinity of $M^{C}$. Moreover, $\hat{V}$ is homeomorphic to $V \cup (\mathbb{R} \times \partial_{C} M) \cup (\{\infty\} \times M^{C})$.
\end{enumerate}
\end{theorem}

Let $(V,g)$ be a globally hyperbolic warped product spacetime. Then the Riemannian manifold $(M,h)$ is complete as can be seen in \cite[Theorem 3.68]{beem} and hence the Cauchy boundary $\partial_{C} M = \emptyset$. Therefore, if in addition
$$\int_{0}^{+\infty} \frac{ds}{\sqrt{\alpha(s)}} < \infty ,$$
then by the previous theorem we have that the future causal completion is characterized as $\hat{V} \equiv V \cup (\{+\infty\} \times M)$ and  the causal boundary consists of a copy of $M$ at infinity. As a consequence, any TIP $P \in \hat{\partial} V$ can be identified with a set of the form $P = I^{-}(+\infty, x)$, $x \in M$.\footnote{Here, $I^{-}(+\infty,x):=\{(t_{0},x_{0})  \in V \mid \int_{t_{0}}^{+\infty} \frac{ds}{\sqrt{\alpha(s)}}< \int_{0}^{+\infty} \frac{ds}{\sqrt{\alpha(s)}} - d(x_{0},x) \}$}

\begin{proposition}\label{crono-es-hausdorff}
Let $(V,g)$ be a globally hyperbolic warped product spacetime with $\int_{0}^{+\infty} \frac{ds}{\sqrt{\alpha(s)}} < \infty$. Then the chronological topology $\hat{\mathcal{T}}_{chr}$ in $\hat{V}$ is Hausdorff.
\end{proposition}

\begin{proof}
Recall that in $i(V)$, the chronological topology agrees with the topology of the space-time $V$, which is Hausdorff. Thus there are only two cases left to consider:
\begin{enumerate}
\item $P$ is PIP and $Q$ is a TIP,
\item $P$ and $Q$ are both TIPs.
\end{enumerate}

\begin{figure}[ht!]
\centering \includegraphics[scale=0.2]{caso1hausdorff.png}
\centering{\caption{Case (1). $P$ is proper and $Q$ is terminal.}}
\end{figure}

\textbf{Case 1.} Consider $p = (t_1 , x_1)$ a point in $V$ and $P = I^{-} (p)$ its correspondent past set; and take $Q = (+\infty , x_2) \in \hat{\partial} V$. We know, by the definition of the chronological relation in $V$, that for $t_2 > t_1$, we have
$$(t_1 , x_1) \ll (t_2 , x_1) $$
because
$$ 0 = d (x_1, x_1) < \int_{t_1}^{t_2} \frac{ds}{\sqrt{\alpha (s)}}.$$
Note that $(t_2,x_2) \in Q$ since $(t_2 ,x_2) \ll (t, x_2)$ for all $t>t_2$. Hence, there exists $q_n \in Q$ with $(t_2 , x_2) \ll q_n$ and as a consequence ${I}^{-}(t_2,x_2) \hat{\ll} Q$. Therefore $Q \in \hat{I}^{+} (I^{-} (t_2, x_2)).$ On the other hand, since $(t_1, x_1) \ll (t_2,x_1)$ then $P \hat{\ll} I^{-} (t_2,x_1)$, and by definition $P \in \hat{I}^{-} ( I^{-} (t_2,x_1))$. Now we need to prove
$$\hat{I}^{+} (I^{-} (t_2, x_2)) \cap \hat{I}^{-} ( I^{-} (t_2,x_1)) = \emptyset.$$
By contradiction, suppose there exists $R = I^{-} (r,y) \in \hat{V}$, with $r \in \mathbb{R}$, such that
$$R \in \hat{I}^{+} (I^{-} (t_2, x_2)) \cap \hat{I}^{-} ( I^{-} (t_2,x_1)) = \emptyset.$$
Note that $r < \infty$ because otherwise $R$ and $Q$ would be two TIPs with $R \hat{\ll} Q$, which cannot occur since $V$ is globally hyperbolic. Thus, $R$ is proper and
$$I^{-} (t_2, x_2) \hat{\ll} I^{-} (r,y) \hat{\ll} I^{-} (t_2, x_1)$$
which happens by Proposition \ref{relacionesextendidas} if and only if
$$(t_2,x_2) \ll (r,y) \ll (t_2, x_1),$$
and by transitivity
$$(t_2, x_2) \ll (t_2, x_1),$$
that is,
$$0 \leq d(x_2,x_1) < \int_{t_2}^{t_2} \frac{ds}{\sqrt{\alpha (s)}} = 0,$$
a contradiction. This concludes case 1.

\begin{figure}[ht!]
\centering \includegraphics[scale=0.2]{caso2hausdorff.png}
\caption{Case 2. $P$ and $Q$ are both terminal.}
\end{figure}

\textbf{Case 2.} Consider $P,Q \in \hat{\partial}V$ with $P \neq Q$. We can represent these sets as
$$P = I^{-} (\infty, x_1), \qquad Q = I^{-} (\infty , x_2).$$
Since $\int_{0}^{\infty} \frac{ds}{\sqrt{\alpha (s)}} < +\infty$ we can take $t_1,t_2 \in \mathbb{R}$ such that
$$\int_{t_1}^{\infty} \frac{ds}{\sqrt{\alpha (s)} }< \frac{1}{3} d(x_1, x_2) \quad \text{and} \quad \int_{t_2}^{\infty} \frac{ds}{\sqrt{\alpha (s)} }< \frac{1}{3} d(x_1, x_2)$$
which implies $(t_1, x_1) \notin Q$ and $(t_2,x_2) \notin P$, respectively.

Without loss of generality suppose $t_1 < t_2$. We know that
$$P \in \hat{I}^{+} (I^{-} (t_2, x_1)), \qquad Q \in \hat{I}^{+} (I^{-} (t_2, x_2)),$$
so it is enough to prove that this sets are disjoint, i.e.
$$\hat{I}^{+} (I^{-} (t_2, x_1)) \cap \hat{I}^{+} (I^{-} (t_2, x_2)) = \emptyset.$$
We proceed by contradition. Take $R = I^{-} (r,y) \in \hat{V}$ such that
$$R \in \hat{I}^{+} (I^{-} (t_2, x_1)) \cap \hat{I}^{+} (I^{-} (t_2, x_2)),$$
whence
$$ I^{-} (t_2, x_1) \hat{\ll} I^{-} (r,y) \quad \text{and} \quad I^{-} (t_2, x_2) \hat{\ll} I^{-} (r,y) $$
by Proposition \ref{relacionesextendidas} if and only if
$$ (t_2,x_1) \ll (r,y) \quad \text{and} \quad (t_2,x_2) \ll (r,y)$$
which happens if and only if
$$d(x_1, y) < \int_{t_2}^{r} \frac{ds}{\sqrt{\alpha (s)}} \quad \text{and} \quad d(x_2, y) < \int_{t_2}^{r} \frac{ds}{\sqrt{\alpha (s)}}$$
where adding both inequalities and using the triangle inequality for $d$ we get
$$d (x_1, x_2) \leq d(x_1,y) + d(x_2,y) < 2 \int_{t_2}^{r} \frac{ds}{\sqrt{\alpha (s)}}$$
and so
$$\int_{t_2}^{\infty} \frac{ds}{\sqrt{\alpha (s)}}  < \int_{t_2}^{r} \frac{ds}{\sqrt{\alpha (s)}},$$
where, being $r$ finite or infinite, we get a contradiction.
\end{proof}

The next result is an immediate consequence of Proposition \ref{crono-es-hausdorff} and Theorem \ref{teo:opagree}.

\begin{coro}
Let $(V,g)$ be a globally hyperbolic warped product spacetime with $\int_{0}^{+\infty} \frac{ds}{\sqrt{\alpha(s)}} < \infty$. Then $\hat{\mathcal{T}}_{chr}=\hat{\mathcal{T}}_{c}$.
\end{coro}

An important fact to highlight  is that on these spacetimes, when $P = I^{-} (p)$ is a PIP and $Q$ is a TIP, their time separation is either zero or infinity.

\begin{proposition}
If $P,Q \in \hat{V}$ with $Q \in \hat{\partial} V$ then either $\hat{\tau} (P,Q) = 0$ or $\hat{\tau} (P,Q) = \infty$.
\end{proposition}

\begin{proof}
If $P \hat{\not\ll} Q$ or $P \in \hat{\partial} V$ then $\hat{\tau} (P,Q) = 0$ and we are done. Now let us assume that $ P = I^{-}(p) \hat{\ll} Q = I^{-} (\infty, a )$ and that $\{q_n\}$ is a future-directed chain generating $Q$. By chronology, there exists $q_N \in Q$ with $p \ll q_N \ll q_n$ for all $n > N$ and since $V$ is globally hyperbolic, $\tau (p,q_N) < \infty$. Given that $\hat{\tau}$ does not depend on the election of the chain that generates Q, therefore we can choose a chain such that
$$q_n = (t_{n}, a),$$
which is based on the same spatial constant coordinate $a \in M$ as $Q$. From here, since $Q$ is terminal, $q_n$ does not converge in $V$ and then $t_{n} \to \infty$. This implies\footnote{Note that  if $\gamma: [t_1,t_2] \to V$ is a vertical curve, that is, $\gamma (t) = ( t , a)$ with $a \in M$ fixed, then $L_g (\gamma) = \int_{t_1}^{t_2} \mid \dot{\gamma} \mid dt =  t_2 - t_1$. This indicates that the $g$-length of $\gamma$ grows to infinity as $t_2\to\infty$.}
$$\tau (p,q_n) \geq \tau (p, q_N) + \tau (q_N, q_n) > \tau (q_N,q_n) \geq t_n - t_N.$$
Thus $\hat{\tau}(P,Q) = \infty$ for any $P$ a PIP and $Q$ a TIP chronologically related.
\end{proof}

We end up this section providing an explicit example. Consider de Sitter spacetime $\mathbb{S}^4_1$, the Lorentzian spaceform of constant sectional curvature $K=1$. It can be realized as the hyperboloid
$$-v^2 + w^2 + x^2 + y^2 + z^2 = 1$$
in flat five-dimensional Minkowski space $\mathbb{R}^5_1=(\mathbb{R}^5,L)$. Recall $\mathbb{S}^4_1$ can also be described as the globally hyperbolic warped product
$$(\mathbb{R} \times \mathbb{S}^3, -dt^2 + \cosh ^2 (t) g_{\mathbb{S}^3}),$$
where $g_{\mathbb{S}^3}$ is the standard round metric on $\mathbb{S}^3$. Thus, $$\int_{0}^\infty \frac{ds}{\sqrt{\alpha(s)}} = \int_{0}^\infty \frac{ds}{\cosh (s)} = \frac{\pi}{2} < +\infty.$$
Moreover, by Theorem \ref{XX} the future causal boundary $\hat{\partial} \mathbb{S}^4_1$ is a copy at infinity of the base manifold $\mathbb{S}^3$, a picture that agrees with its standard future conformal infinity $\mathcal{I}^+$ \cite{HW}.

Let $p\leq q$ be two causally related points in $\mathbb{S}^4_1$. Since the time separation function is given by
$$\cosh (\tau (p,q)) = L(p,q),$$
then
$$\tau (p,q) = t_2 - t_1$$
for $p=(t_1,x)$, $q=(t_2,x)$ two points on the same geodesic normal to $\mathbb{S}^3$. Thus, if $t_2 \to \infty$ then $\tau (p,q) \to \infty$. If we consider a PIP $P = I^{-} (t,x)\in \hat{\mathbb{S}}^4_1$  and a TIP $Q = I^{-} (\infty , y) \in \hat{\partial} \mathbb{S}^4_1$ a TIP, then by definition
$$\hat{\tau} (P,Q) = \sup \{\tau (p,q_n): Q = I^{-}(\{ q_n \})\},$$
where $\{ q_n\}$ is a future-directed chain generating $Q$. Choose a chain such that
$$q_n = (t_{n}, y),$$
which is based on the same spatial coordinate $y \in \mathbb{S}^3$. From here, since $Q$ is terminal, $q_n$ does not converge in $\mathbb{S}^4_1$ and then $t_{n} \to \infty$

We know that if $p=(t,x) \notin Q$ or $p \in \partial Q$ we have $\hat{\tau} (I^{-} (p), Q) = 0$ by Lemma \ref{lema:rev}. On the other hand, if $p \in Q$, there exists $N \in \mathbb{N}$ with $p \ll q_N \ll q_n$ for all $n > N$. By the reverse triangle inequality for $\tau$ we have
$$\tau (p,q_n) \geq \tau (p,q_N) + \tau(q_N,q_n) > \tau(q_N,q_n) \geq t_n - t_N,$$
where the second term on the right-hand-side goes to infinity as $n \to \infty$. Therefore
$$\hat{\tau} (P,Q) = \infty.$$

In summary, the time separation function $\hat{\tau}$ for the future causal completion of de Sitter spacetime is given by
\begin{equation*}
\hat{\tau} (P,Q) =\left\lbrace\begin{array}{ccc}
\tau (p,q) & \text{if} & P = I^{-}(p), Q= I^{-}(q) \\
0 & \text{if} & P \hat{\nleq} Q \\
\infty & \text{if} & P \hat{\ll} Q \text{ and } Q \in \hat{\partial} \mathbb{S}^4_1.
\end{array}\right.
\end{equation*}

\section*{Acknowledgements} The authors would like to thank the organizing committee of SCRI21, for providing a source of great mathematical ideas for years to come. L. Ake Hau was partially supported by Conacyt under grants SNI 367994 and Estancia Posdoctoral por M\'exico 854544. S. Burgos acknowledges the support of Conacyt under the Becas Nacionales program (1089466). D. Solis was partially supported by Conacyt under grant SNI 38368.


\section{Introduction}
Trajectory inference aims at recovering the dynamics of a population of particles given  samples from its temporal marginals at various time-points. This problem arises notably in the analysis of single-cell RNA-sequencing data~\citep{schiebinger2019optimal,tong2020trajectorynet,farrell2018single}, where one has access---via a destructive measurement process---to the cell state of samples from a large population of cells that evolve in time. In this application, one would like to recover the overall dynamics of the population as well as the trajectories of individual cells so as to improve our understanding of certain biological processes, such as embryonic development or tumor progression.

Trajectory inference can be cast as a regression problem where the unknown is the law of a continuous stochastic process.
Let $\mathcal{X}$ be the ambient space containing the particles and $\Omega\coloneqq \mathcal{C}([0,1],\mathcal{X})$ the path-space, i.e. the set of all possible continuous trajectories over the time interval $[0,1]$.
We consider the problem of recovering a probability distribution over paths $P\in \mathcal{P}(\Omega)$, i.e. the law of a continuous stochastic process $(X_t)_{t\in [0,1]}$.
For the inference to be well-behaved, one should look for a distribution $P$ such that its time marginals are consistent with the observed snapshots and such that the overall dynamics it represents satisfies some notion of regularity.

The problem of trajectory inference has received significant attention over the past several years.
However, while pioneering work has focused on creating new methods, the theoretical treatment has remained limited.
One class of methods focuses on recovering a potential energy landscape that best fits the observed marginals. For example, Hashimoto \emph{et al.}~\cite{hashimoto2016learning} encode the potential with a neural network, TrajectoryNet~\cite{tong2020trajectorynet} encodes the potential as a neural ODE, and JKOnet~\cite{bunne2021jkonet} encodes the potential with a neural network architecture based on the JKO scheme~\cite{jordan1998variational} and input convex neural networks.
However, these current approaches are all nonconvex, and therefore it can be difficult to establish rigorous guarantees.

Some guarantees have been established for the equilibrium case by Weinreb \emph{et al.}~\cite{weinreb2018fundamental}, and recently consistency was established for the non-equilibrium case by Lavenant \emph{et al.}~\cite{lavenant2021towards}, who, as a regularity prior, penalize the entropy of $P$ relative to the Wiener measure on $\Omega$ (i.e. the law of the Brownian motion).
Moreover, Bunne \emph{et al.} have recently shown how to produce an improved reference process, beyond ordinary Brownian motion \cite{bunne2022recovering}.

In this work we focus on the estimator introduced by Lavevant \emph{et al.}~\cite{lavenant2021towards}, which reconstructs the process $P$ by minimizing the entropy relative to the Wiener measure over $\Omega$.
By trading off data fitting with regularization, this approach generalizes the approach of Waddington-OT~\cite{schiebinger2019optimal} to datasets with many time-points yet possibly few samples per time-point, which is required for consistency~\cite{lavenant2021towards}, yields better performance in practice~\cite{schiebinger2021recovering}, and is achievable through single-embryo profiling~\cite{MITTNENZWEIG20212825}.
While the approach of Lavevant \emph{et al.} leads to a consistent estimator, it also leads to a challenging infinite dimensional convex optimization problem. It is tackled in the original paper by discretizing the space and then applying convex optimization methods. The goal of the present paper is to show that the specific structure of this estimator makes it amenable to a newly introduced class of grid-free stochastic methods called \emph{Mean-Field Langevin} (MFL) dynamics~\citep{hu2019mean, mei2018mean}, a non-linear generalization of Langevin dynamics which enjoys quantitative global convergence guarantees~\citep{nitanda2022convex, chizat2022mean}. Instantiated in our context, this method consists in a family of point clouds (one per snapshot) coupled via entropy regularized optimal transport, a.k.a.~Schr\"odinger bridges~\cite{leonard2012schrodinger}, which evolve with noisy gradient descent.
Intuitively, these dynamics can be interpreted as a (non-linear) Langevin diffusion over the path space $\Omega$.
Our approach leads to an inference method with end-to-end theoretical guarantees that solves an interpretable model for trajectory inference. We illustrated on Figure~\ref{fig:cover} the recovered estimator and on Figure~\ref{fig:dynamics} the MFL optimization dynamics.

\paragraph{Organization of the paper} In \S~\ref{sec:estimator}, we present the estimator of~\citet{lavenant2021towards}. The heart of our theoretical contributions is in \S~\ref{sec:algorithm} where we introduce the MFL dynamics and show its convergence towards the min-entropy estimator at an exponential rate. Numerical experiments are shown in \S~\ref{sec:numerics}.

\paragraph{Notation and blanket assumptions} For two probability measures $\mu,\nu$, their relative entropy is $H(\mu|\nu)=\int \log(\d\mu/\d\nu)\d\mu$ if $\mu \ll \nu$ and $+\infty$ otherwise. For $n\in \mathbb{N}$, let $[n]\coloneqq \{1,\dots,n\}$. Throughout, the ambient space $\mathcal{X}$ is a compact convex subset of $\mathbb{R}^d$ and $(B_t)_t$ is a reflected Brownian motion on $\mathcal{X}$ with $B_0$ distributed uniformly over $\mathcal{X}$. The path space is $\mathcal{C}([0,1];\mathcal{X})$ and laws on path space are noted by capital letters, $P$ or $R$. We denote by $P_{t_1,\dots, t_T} = (e_{t_1},\dots,e_{t_T})_\# P$ the marginal of such laws under the (joint) evaluation maps $e_t:\omega \mapsto \omega(t)$. We use boldface letters for family of probability measures $\boldsymbol{\mu}\in \mathcal{P}(\mathcal{X})^T$.

\begin{figure}
    \centering
    \begin{subfigure}{0.62\linewidth}
        \includegraphics[width = \linewidth]{fig1_langevin_new.pdf}
    \end{subfigure}
    \begin{subfigure}{0.22\linewidth}
        \includegraphics[width = \linewidth]{fig1_distances.pdf}

    \end{subfigure}
    \caption{(left) Sampled time-series data and reconstructions by the proposed Mean-Field Langevin dynamics (MFL) and Global Waddington-OT~\cite{lavenant2021towards} (gWOT). Observed and reconstructed marginals are colored by the measurement time $t_i$, and illustrative sample paths are overlaid in grey. We compare a scenario with limited data at intermediate time-points ($N = 1$) to a uniformly sampled scenario ($N = 64$). (right) RMS Energy Distance to the ground truth over marginals for $N = 2^0, \ldots, 2^6$, illustrating the robustness of MFL to the low-sample regime. See \S\ref{sec:fig1_example} for details.}
    \label{fig:cover}
\end{figure}

\section{Min-Entropy Estimator in Path Space}\label{sec:estimator}

\subsection{Trajectory Inference as Stochastic Process Inference}
\paragraph{Model of population dynamics} The first step of any inference task is to determine a prior on the ground truth. For population dynamics with mass conservation, a natural prior is given by drift-diffusion processes. We thus model the population dynamics as
\begin{equation}\label{eq:SDE}
    \d X_t = -\nabla \Psi(t,X_t) \d t +\sqrt{\tau} \d B_t,\qquad \mathrm{Law}(X_0) = \mu_0,\qquad t\in [0,1],
\end{equation}
where $\tau>0$ is the temperature/diffusivity, assumed known, and $\Psi\in \mathcal{C}^2([0,1]\times \mathcal{X})$ is unknown. As discussed in prior works~\citep{hashimoto2016learning, lavenant2021towards}, it is hopeless to recover the divergence-free component of the drift: it is thus natural to assume that the drift is given by the gradient of a function $\Psi$, called the \emph{Waddington potential} or \emph{epigenetic landscape} in the context of cell development. The noise level $\tau>0$ models the inherent randomness of the ground truth process. It turns out that this noise also has a favorable effect on the algorithm we develop here. This model can be extended to allow particles to branch and die~\citep{baradat2021regularized}, this extension is discussed in Section~\ref{sec:mass-variation}.

\paragraph{Model of measurements} Let $0 \leq t_1<\dots <t_T\leq 1$ be a family of measurement times.
In contrast to the classical field of inference for stochastic processes~\citep{rao2000stochastic} where each realization $(X_t)_{t\in [0,1]}$ is observed $T$ times (with $T$ large), in trajectory inference each realization is only observed once, i.e.~we observe $(X_{t_i,j})_{i\in [T], j\in [n_i]}$ with independence for each couple $(i,j)$. The measurements are summarized by the family of snapshots $\hat \mu_{t_i} = \frac{1}{n_i}\sum_{j=1}^{n_i} \delta_{X_{t_i,j}} \in \mathcal{P}(\mathcal{X})$ for $i\in [T]$, where $n_i\geq 1$ is the number of particles observed at time $t_i$.

\paragraph{Inference problem} The general goal is to recover the law $P=\mathrm{Law}(X)\in \mathcal{P}(\Omega)$ of the stochastic process given the snapshots $(\hat \mu_{t_1},\dots, \hat \mu_{t_T})$. This object $P$ contains all the information about $X$: in particular it contains the marginals $P_t$ of the process and the transition probability between each family of marginals $P_{t_1,\dots,t_T}$. Note that, as discussed in the introduction, some works~\citep{tong2020trajectorynet,hashimoto2016learning,bunne2021jkonet} focus on directly on recovering the Waddington potential $\Psi$ within a parameterized class of functions. This is a different problem, with encouraging empirical results but weaker known theoretical guarantees.

\subsection{Min-Entropy Estimator}
Let $W^\tau \in \mathcal{P}(\Omega)$ be the law of the Brownian motion with temperature $\tau$. The min-entropy estimator\footnote{We propose this name because $R^*$ is the probability measure in path space with the smallest entropy relative to the Wiener measure among all those that fit at least equally well the observations.} introduced in~\cite{lavenant2021towards}, is defined as the unique minimizer $R^*$ of the functional $\mathcal{F}:\mathcal{P}(\Omega)\to \mathbb{R}$ defined as
\begin{align}\label{eq:estimator}
\mathcal{F}(R) \coloneqq \frac{1}{\lambda} \sum_{i=1}^T \Delta t_{i}\mathrm{Fit}_\sigma(R_{t_i}|\hat \mu_{t_i}) + \tau H(R|W^{\tau})
\end{align}
where $\lambda>0$ is the regularization strength, $\Delta t_i \coloneqq (t_{i+1}-t_{i-1})/2$ (with the convention $t_{-1}=0$ and $t_{T+1}=1$) and $\mathrm{Fit}_\sigma$ is a divergence functional that quantifies how much $R_{t_i}$ and $\hat \mu_{t_i}$ differ, parameterized by a bandwidth $\sigma>0$. For a particular choice of fitting functional (see below), the authors prove the following result:
\begin{theorem}[Consistency,~\cite{lavenant2021towards}]\label{thm:consistency}
If $(t_i)_{i\in [T]}$ becomes dense in $[0,1]$ as $T$ grows, then
\[
\lim_{\lambda, \sigma \to 0}\lim_{T\to \infty} R^* = P \qquad \text{weakly, almost surely.}
\]
\end{theorem}
In this paper, we consider the following data fitting term
\begin{align}\label{eq:fit}
\mathrm{Fit}_\sigma(R_t,\hat \mu_t)
&\coloneqq \int -\log \left[ \int \exp\Big( -\frac{\Vert x-y\Vert^2}{2\sigma^2} \Big) \d R_t(x) )\right]\d\hat \mu_t(y)\\
& = H(\hat \mu_t|R_t \ast g_{\sigma}) + H(\hat \mu_t) + C\label{eq:not-always-defined}
\end{align}
where $g_{\sigma}(x) = (2\pi\sigma^2)^{-\frac{d}{2}}e^{-\Vert x\Vert^2/(2\sigma^2)}$ is the Gaussian density and $C>0$ is a constant. Eq.~\eqref{eq:fit} defines $\mathrm{Fit}_\sigma$ and is always well-defined and finite, but Eq.~\eqref{eq:not-always-defined} only makes sense when $H(\hat \mu_t)<\infty$.
This is the negative log-likelihood under the noisy observation model $\hat X_{t_i,j}=X_{t_i,j}+\sigma Z_{i,j}$ where $\hat X_{t_i,j}$ is the observation and $Z_{i,j}\overset{iid}{\sim} \mathcal{N}(0,I)$ the noise. Intuitively, when minimizing this \emph{soft-min} objective over measures $R_t$ which are mixtures of Dirac masses, each of the observed particles from $\hat \mu_{t}$ tends to attract the particles from $R_t$ that are the closest---even if they are far away in an absolute sense---but barely influences those that are the farthest.

This data-fitting functional $\mathrm{Fit}_\sigma$ slightly differs from $H(g_\sigma \ast \hat \mu_t|R_t)$, which is the one for which Thm.~\ref{thm:consistency} has been established in~\cite{lavenant2021towards}, but it preserves its key properties, namely joint convexity in $(R_t,\hat \mu_t)$ and linearity in $\hat \mu_t$. Our choice for $\mathrm{Fit}_\sigma$ is motivated by a more favorable behavior when $R_t$ is a discrete measure and by its natural statistical interpretation. In this paper, we focus on the optimization aspects and leave the statistical analysis of the estimator $R^*$ with this specific choice of $\mathrm{Fit}_\sigma$ for future work. Let us also mention that Thm.~\ref{thm:consistency} has been established for $\mathcal{X}$ a compact manifold without boundary, while under our assumptions $\mathcal{X}$ has boundaries.

We now turn to the presentation of the algorithm. The data model of Eq.~\eqref{eq:SDE} was introduced only to motivate the min-entropy estimator via Thm.~\ref{thm:consistency}, and plays no role in the rest of the paper.

\section{Optimization Method: Mean Field Langevin Dynamics}\label{sec:algorithm}
In order to compute the min-entropy estimator $R^*\in \mathcal{P}(\Omega)$, one needs to minimize an entropy regularized functional over the space of paths $\mathcal{F}:\mathcal{P}(\Omega)\to \mathbb{R}$ which is of the form
 \begin{align}\label{eq:path-objective}
 \mathcal{F}(R) = \mathrm{Fit}(R_{t_1},\dots,R_{t_T}) +  \tau H(R|W^{\tau})
 \end{align}
where we have posed $\mathrm{Fit}(R_{t_1},\dots,R_{t_T}):=\frac{1}{\lambda} \sum_{i=1}^T \Delta t_{i}\mathrm{Fit}_\sigma(R_{t_i}|\hat \mu_{t_i})$. This is the sum of a convex functional $\mathrm{Fit}$ and the relative entropy with respect to the scaled Wiener measure. If, instead of $\mathcal{P}(\Omega)$, the optimization space was $\mathcal{P}(\mathcal{X})$, a problem with this structure could be solved by Mean-Field Langevin (MFL) dynamics. More explicitly, MFL dynamics are designed to minimize problems of the form $F(\mu) = G(\mu)+\tau H(\mu)$ where $G:\mathcal{P}(\mathcal{X})\to \mathbb{R}$ is convex and ``smooth'' and $H$ is minus the differential entropy. Considering for some $m\in \mathbb{N}^*$ the noisy gradient descent on the function $(x_1,\dots,x_m) \mapsto G(\frac1m \sum_{i=1}^m \delta_{x_i})$,  MFL dynamics are obtained in the mean-field $m\to \infty$ and vanishing step-size limit, see Section~\ref{sec:MFL} for details.

Inspired by this parallel, we now design a drift-diffusion dynamics in path space that converges to the minimizer of $\mathcal{F}$. The main idea is a reformulation of the problem as a ``reduced'' problem over $\mathcal{P}(\mathcal{X})^T$ with a structure amenable to MFL dynamics.

\subsection{Reduced Formulation}
 We first introduce the reduced functional $F$, and then state its connection to $\mathcal{F}$ in Thm.~\ref{thm:representer}.

For $\mu,\nu \in \mathcal{P}(\mathcal{X})$, let $\Pi(\mu,\nu)$ be the set of transport plans between $\mu$ and $\nu$, that is, probability measures on $\mathcal{X}\times \mathcal{X}$ with respective marginals $\mu$ and $\nu$. The entropy regularized optimal transport cost between $\mu$ and $\nu$ is defined, for some $\tau_i>0$, as
\begin{align}\label{eq:ROT}
    T_{\tau_i}(\mu,\nu) &\coloneqq \min_{\gamma \in \Pi(\mu,\nu)} \int c_{\tau_i}(x,y) \d\gamma(x,y) + \tau_i H(\gamma|\mu\otimes \nu) = \min_{\gamma \in \Pi(\mu,\nu)} \tau_i H(\gamma|p_{\tau_i}\mu\otimes \nu).
\end{align}
where $p_t(x,y)$ is the transition probability density of the Brownian motion on $\mathcal{X}$ over the time interval $[0,t]$ and $c_{\tau_i}(x,y)\coloneqq -\tau_i \log(p_{\tau_i}(x,y))$. In numerical experiments we use the approximation suggested by Varadhan's formula when $\tau_i$ is small: $c_{\tau_i}(x,y)\approx \frac12 \Vert y-x\Vert^2$ (up to an additive constant which is irrelevant when one is interested in minimizers only)~\cite{norris1997heat}.
This optimization problem is also known as the \emph{Schr\"odinger bridge problem}: its solution gives the most likely evolution of a cloud of particles following a Brownian motion, conditioned on being distributed as $\mu$ at $t=0$ and $\nu$ at $t=1$~\citep{leonard2012schrodinger}. Some background on this problem is given in Appendix~\ref{app:EOT}, including the definition of the Schr\"odinger potentials $(\varphi,\psi)$ used hereafter. Consider the function $G:\mathcal{P}(\mathcal{X})^T\to \mathbb{R}$ defined for $\boldsymbol{\mu} = (\boldsymbol{\mu}^{(1)},\dots,\boldsymbol{\mu}^{(T)})$ that represents a family of $T$ reconstructed temporal marginals, by
\begin{align}\label{eq:G-functional}
    G(\boldsymbol{\mu}) \coloneqq \mathrm{Fit}(\boldsymbol{\mu}) +  \sum_{i=1}^{T-1} \frac{1}{t_{i+1}-t_i} T_{\tau_i}(\boldsymbol{\mu}^{(i)},\boldsymbol{\mu}^{(i+1)})
\end{align}
where $\tau_i \coloneqq (t_{i+1}-t_i)\tau$.
 We now introduce the \emph{reduced objective} $F:\mathcal{P}(\mathcal{X})^T\to \mathbb{R}$, defined as
\begin{align}\label{eq:reduced-objective}
F(\boldsymbol{\mu}) \coloneqq G(\boldsymbol{\mu}) + \tau H(\boldsymbol{\mu})
\end{align}
where $H(\boldsymbol{\mu})=\sum_{i=1}^T \int \log(\boldsymbol{\mu}^{(i)}(x))\d\boldsymbol{\mu}^{(i)}(x)$ is minus the differential entropy of the family of measures $\boldsymbol{\mu}$.

The next result makes the link between minimizing $\mathcal{F}$, the objective in path space~\eqref{eq:path-objective}, and $F$, the reduced objective~\eqref{eq:reduced-objective}. It can be interpreted as a Representer Theorem for the min-entropy estimator.\footnote{We stretch a bit the term ``representer theorem'' which is usually reserved to finite-dimensional reductions.} This theorem is straightforward to deduce from~\cite[Cor.~3.5]{benamou2019entropy}, but the specific form of $F$ that we use here is central for our subsequent algorithmic developments.\footnote{The key difference with~\cite[Prop.~B.1]{lavenant2021towards} is that here $T_\tau$ involves the entropy of $\gamma$ relative to the product measure (instead of the volume measure), which makes $H(\boldsymbol{\mu})$ appear with a \emph{positive} sign in $F$ (instead of negative).}

\begin{theorem}[Representer Theorem]\label{thm:representer}
Let $\mathrm{Fit}:\mathcal{P}(\mathcal{X})^T\to \mathbb{R}$ be \emph{any} function.\begin{enumerate}[label=(\roman*)]
    \item If $\mathcal{F}$ admits a minimizer $R^*$ then $(R^*_{t_1},\dots,R^*_{t_T})$ is a minimizer for $F$.
    \item Conversely, if $F$ admits a minimizer $\boldsymbol{\mu}^* \in \mathcal{P}(\mathcal{X})^T$ then a minimizer $R^*$ for $\mathcal{F}$ is built as
\[
R^*(\cdot) = \int_{\mathcal{X}^T} W^\tau(\cdot| x_1,\dots,x_T )\d R_{t_1,\dots,t_T}(x_1,\dots,x_T)
\]
 where $W^\tau(\cdot|x_1,\dots,x_T)$ is the law of $W^\tau$ conditioned on passing through $x_1,\dots,x_T$ at times $t_1,\dots,t_T$ respectively and $R_{t_1,\dots,t_T}$ is the composition of the transport plans $\gamma_{i,i+1}$ which are optimal in the definition of $T_{\tau_i}(\boldsymbol{\mu}^{*(i)},\boldsymbol{\mu}^{*(i+1)})$, for $i=1,\dots,T-1$.
\end{enumerate}
\end{theorem}

The composition of the transport plans is obtained as
\[
R_{t_1,\dots,t_T}(\d x_1,\dots, \d x_T) = \gamma_{1,2}(\d x_1,\d x_2) \gamma_{2,3}(\d x_3|x_2) \dots \gamma_{T-1,T}(\d x_{T}|x_{T-1})
\]
where we have introduced the conditional probability (a.k.a.~``disintegrations'') characterized by $\gamma_{i,i+1}(\d x_{i},\d x_{i+1}) = \gamma_{i,i+1}(\d x_{i+1}|x_i)\mu_i(\d x_i)$.
The proof of Theorem~\ref{thm:representer} is given in Appendix~\ref{app:proof-representer}.

The importance of the reduced problem comes from the following facts:
\begin{itemize}
    \item The optimization space has been ``reduced'' from $\mathcal{P}(\Omega)$ to $\mathcal{P}(\mathcal{X})^T$. This reduction is enabled by the Markovian property of $W^\tau$. Moreover, Theorem~\ref{thm:representer} gives an explicit method to construct a minimizer for $\mathcal{F}$ from a minimizer for $F$ and the associated optimal transport plans $(\gamma_{i,i+1})_{i\in [T-1]}$. This fact was already exploited in~\citet{lavenant2021towards}.
    \item The objective function $F$ is the sum of two terms: a ``smooth'' convex function $G$ and a  differential entropy term $\tau H$. This is precisely the structure tackled by MFL dynamics. Observe how the entropy in path space $H(P|W^\tau)$ is split into two parts: the Schr\"odinger bridges $T_\tau$, included in the ``smooth'' term $G$, and minus the differential entropy $H$.
\end{itemize}

Let us now describe some useful properties of $G$ and $F$. Hereafter, the \emph{first-variation} of $G:\mathcal{P}(\mathcal{X})^T\to \mathbb{R}$ at $\boldsymbol{\mu}$ is the unique function $V[\boldsymbol{\mu}]\in \mathcal{C}(\mathcal{X})^T$ such that for all $\boldsymbol{\nu} \in \mathcal{P}(\mathcal{X})^T$,
\begin{align}\label{eq:first-variation}
\lim_{\epsilon \downarrow 0} \frac1\epsilon \big(G((1-\epsilon)\boldsymbol{\mu}+\epsilon \boldsymbol{\nu}) -G(\boldsymbol{\mu})\big) = \sum_{i=1}^T \int V^{(i)}[\boldsymbol{\mu}](x) \d(\boldsymbol{\nu}-\boldsymbol{\mu})^{(i)}(x).
\end{align}

\begin{proposition}\label{prop:first-properties}
The function $G$ is convex, weakly continuous and its \emph{first-variation} is given for $\boldsymbol{\mu}\in \mathcal{P}(\mathcal{X})^T$ and $i\in [T]$ by
\begin{align*}
    V^{(i)}[\boldsymbol{\mu}] = \frac{\delta \mathrm{Fit}}{\delta \boldsymbol{\mu}^{(i)}}[\boldsymbol{\mu}]+\frac{\varphi_{i,i+1}}{t_{i+1}-t_i} + \frac{\psi_{i,i-1}}{t_{i}-t_{i-1}},\quad \frac{\delta \mathrm{Fit}}{\delta \boldsymbol{\mu}^{(i)}}[\boldsymbol{\mu}]:x \mapsto  - \frac{\Delta t_i}{\lambda} \int \frac{g_\sigma(x-y)}{(g_{\sigma}\ast \boldsymbol{\mu}^{(i)})(y)}\d\hat \mu_{t_i}(y)
\end{align*}
where $(\varphi_{i,j},\psi_{i,j})\in \mathcal{C}^\infty(\mathcal{X})$ are the Schr\"odinger potentials for $T_{\tau_i}(\boldsymbol{\mu}^{(i)},\boldsymbol{\mu}^{(j)})$, with the convention that the corresponding term vanishes when it involves $\psi_{1,0}$ or $\varphi_{T,T+1}$.
The function $F_\tau$ admits a unique minimizer $\boldsymbol{\mu}^*$, which has an absolutely continuous density (again denoted by $\boldsymbol{\mu}^*$) characterized by
\[
(\boldsymbol{\mu}^*)^{(i)} \propto e^{-V^{(i)}[\boldsymbol{\mu}^*]/\tau}, \quad \text{for $i\in [T]$}.
\]
\end{proposition}

From now on, we thus focus on minimizing the reduced functional $F$ in Eq.~\eqref{eq:reduced-objective}, keeping in mind that this is sufficient to minimize $\mathcal{F}$, which is our main goal.

\subsection{Mean-Field Langevin Dynamics}\label{sec:MFL}
As previously mentioned, the MFL dynamics is natural when it comes to minimize functionals of the form $F=G+\tau H$. Using the first-variation $V[\boldsymbol{\mu}]$ of $G$ given in Prop.~\ref{prop:first-properties}, the MFL dynamics is defined as the solution of the following non-linear SDE of McKean-Vlasov type, for $s\geq 0$:
\begin{align}\label{eq:MFL-SDE}
\left\{
\begin{aligned}
    dX_{s}^{(i)} &= - \nabla V^{(i)}[\boldsymbol{\mu}_s](X_{s}^{(i)})\d s + \sqrt{2\tau} dB_{s}^{(i)},\quad \mathrm{Law}(X_{0}^{(i)})=\boldsymbol{\mu}_0^{(i)}\\
    \boldsymbol{\mu}_{s}^{(i)} &= \mathrm{Law}(X_{s}^{(i)}), \quad i\in [T]
    \end{aligned}
    \right.
\end{align}
with reflected boundary conditions in $\mathcal{X}$. Beware that the pseudo-time $s$ of the optimization dynamics in Eq.~\eqref{eq:MFL-SDE} should not be mistaken with the pseudo-time $t$ of the process in Eq.~\eqref{eq:SDE}, which is now represented by the discrete exponents $i\in [T]$.

The family of laws $(\boldsymbol{\mu}_{s})_{s\geq 0}$ of this stochastic process are characterized by the following system of PDEs (understood in the sense of distributions, with no-flux boundary conditions):
\begin{align}\label{eq:MFL}
    \partial_s \boldsymbol{\mu}_{s}^{(i)} = \nabla \cdot (\nabla V^{(i)}[\boldsymbol{\mu}_s]\boldsymbol{\mu}_{s}^{(i)}) + \tau \Delta \boldsymbol{\mu}_{s}^{(i)},\quad i \in [T]
\end{align}
which are coupled with the quantity $\nabla V^{(i)}[\boldsymbol{\mu}_s]$. The link between~\eqref{eq:MFL-SDE} and~\eqref{eq:MFL} follows from Ito's lemma, see e.g.~\cite{hu2019mean}. This is a multi-species PDE where each of the species $\boldsymbol{\mu}^{(i)}$ attempts to minimize $\frac{\Delta t_i}{\lambda}\mathrm{Fit}_\sigma(\cdot|\hat \mu_{t_i}) +\tau H$ via a drift-diffusion dynamics, and is connected to $\boldsymbol{\mu}^{(i-1)}$ and $\boldsymbol{\mu}^{(i+1)}$ via Schr\"odinger bridges.

\subsection{Exponential Convergence}\label{sec:theory}
Let us now state the main convergence result proved in Appendix~\ref{app:convergence}. We note that extending this result to the non-compact case raises interesting challenges and is left as an open question.

\begin{theorem}[Convergence]\label{thm:convergence}
Let $\boldsymbol{\mu}_0\in \mathcal{P}(\mathcal{X})^T$ be such that $F(\boldsymbol{\mu}_0)<\infty$. Then there exists a unique solution $(\boldsymbol{\mu}_s)_{s\geq 0}$ to the MFL dynamics~\eqref{eq:MFL} and it holds
\[
F(\boldsymbol{\mu}_s) - F(\boldsymbol{\mu}^*) \leq e^{-Cs}\big(F(\boldsymbol{\mu}_0) - F(\boldsymbol{\mu}^*)\big).
\]
where $C:= 2\tau D^{-2} e^{-(2\mathrm{osc}(c)/\Delta t + e^{D^2/(2\sigma^2)}\Delta t/\lambda)/\tau}$ with $\Delta t \coloneqq \max_{i} \{ t_{i+1}-t_i\}$ and $D$ the diameter of $\mathcal{X}$, and $\mathrm{osc}(c) \coloneqq \max_i (\sup c_{\tau_i} - \inf c_{\tau_i})$ which, converges to $D^{2}/2$ when $\max_i \tau_i\to 0$.
\end{theorem}

We can make the following comments:
\begin{itemize}
    \item Under these assumptions, the convergence of $\boldsymbol{\mu}_s$ to $\boldsymbol{\mu}^*$ in relative entropy and in Wasserstein distance, also holds, with the same rate~\cite{nitanda2022convex, chizat2022mean}.
    \item By Thm.~\ref{thm:representer}, one can map $(\boldsymbol{\mu}_s)_{s\geq 0}$ to a dynamics in $\mathcal{P}(\Omega)$ which converges to $R^*$. It is this dynamics which we refer to as ``Mean-Field Langevin in Path Space'' in the title.
    \item This result gives a convergence rate for an optimization dynamics in a non-convex landscape, so it is not even obvious to have global convergence. The convergence rate depends exponentially on $1/\tau$, a standard drawback of Langevin-like dynamics in absence of log-concavity. The very poor dependency in $D^2/\sigma^2$ could be mitigated by replacing the Gaussian smoothing kernel $g_\sigma$ by a kernel with a slower decay.
\end{itemize}

\paragraph{Summary of theoretical guarantees:}
Overall, the MFL dynamics of Eq.~\eqref{eq:MFL} converges to the unique minimizer $\boldsymbol{\mu}^*$ of the reduced problem in Eq.~\eqref{eq:reduced-objective}, from which we can build  by Thm~\ref{thm:representer} the min-entropy estimator $R^*\in \mathcal{P}(\Omega)$ which is optimal for~\eqref{eq:path-objective}. In this sense, our method leads to an inference method with end-to-end theoretical guarantees for trajectory inference.
However, this statement is about an idealized dynamics which in practice has to be discretized in time and with particles.
In the next section (\S \ref{sec:discretization}), we will see that the resulting discretization error can be controlled, leveraging the long history of work on mean-field dynamics.

\subsection{Discretization}\label{sec:discretization}
In practice, an approximation of the MFL dynamics is obtained by running noisy gradient descent on the function $G_m:(\mathcal{X}^m)^T\to \mathbb{R}$ defined as
\begin{align}\label{eq:Gm}
    G_m(\hat X) \coloneqq G(\hat \boldsymbol{\mu}_{\hat X}) &&\text{where}&& \hat \boldsymbol{\mu}_{\hat X}^{(i)} = \frac1m \sum_{j=1}^m \delta_{\hat X_j^{(i)}}
\end{align}
where $m\in \mathbb{N}$ is be number of particles used to discretized each of the time marginals $\boldsymbol{\mu}^{(i)}$. Since it can be shown that $m \nabla_{X^{(i)}_j} G_m(\hat X) = \nabla V^{(i)}[\hat \boldsymbol{\mu}_{\hat X}](\hat X^{(i)}_{j})$ (see e.g.~\cite[Prop.~2.4]{chizat2022mean}), this leads to the following update equations, for $k\geq 0$, $i\in [T]$ and $j\in [m]$:
\begin{align}\label{eq:NPGD}
\left\{
\begin{aligned}
    \hat X_{j}^{(i)}[k+1] &= \hat X_{j}^{(i)}[k] - \eta \nabla V^{(i)}[\hat \boldsymbol{\mu}[k]](\hat X^{(i)}_{j}[k]) + \sqrt{2\eta \tau} Z^{(i)}_{j,k},\quad \hat X_{j}^{(i)}[0] \overset{iid}{\sim} \boldsymbol{\mu}_0^{(i)} \\
    \hat \boldsymbol{\mu}^{(i)}[k] &= \frac{1}{m} \sum_{j=1}^m \delta_{\hat X^{(i)}_j[k]}, \quad i\in [T]
    \end{aligned}
    \right.
\end{align}
where $\eta>0$ is a step-size. The convergence of such a numerical scheme towards the MFL dynamics is standard in the mean-field literature~\cite{sznitman1991topics}, and holds here thanks to the stability of $\nabla V$ which follows from~\cite{anonymous2022stability}.

\begin{figure}
    \centering
    \begin{subfigure}{0.22\linewidth}
     \centering
        \includegraphics[width = \linewidth]{fig2_objective.pdf}

    \end{subfigure}
    \begin{subfigure}{0.77\linewidth}
        \includegraphics[width = \linewidth]{fig2_trainingdynamics.pdf}

    \end{subfigure}
    \caption{Optimization dynamics in the setting of Figure~\ref{fig:cover} with $N=64$. (left) Evolution of the objective function $F(\boldsymbol{\mu}_s)$ (Eq.~\eqref{eq:reduced-objective}) (right) Evolution of the reconstructed marginals $\hat \boldsymbol{\mu}^{(i)}[k]$ (Eq.~\eqref{eq:NPGD}) with the iteration number $k$, starting from isotropic Gaussians at $k=1$ and colored by measurement time $t_i$.}
    \label{fig:dynamics}
\end{figure}

In the related context of wide neural networks, quantitative bounds on the discretization error were studied in~\cite[Thm.~2]{mei2019mean}. It is shown that for $s=k\eta$, with high probability, $\vert G(\hat \boldsymbol{\mu}_k) - G(\boldsymbol{\mu}_{s})\vert = \tilde O(e^{Cs}(m^{-1/2}+\eta^{1/2}))$ for some $C>0$ independent of $k,\eta$ and $m$. Inspecting their proof, it can be seen that it goes through as long as $G(\mu)$ and $\nabla V[\mu]$ can be estimated at the ``fast rate'' from independent samples of $\mu$, which is our case here~\cite[Thm.~4.5]{delbarrio2022improved}. However, the exponential dependence in the pseudo-time $s$ of these bounds make them not precise enough to obtain interesting long-time guarantees. We leave for future work a direct convergence analysis of the fully discretized MFL dynamics (note that ~\cite[Cor.~2]{nitanda2022convex} gives a guarantee for the \emph{time}-discretized MFL dynamics).

Each iteration of Eq.~\eqref{eq:NPGD} requires to solve an entropic optimal transport problem in order to obtain $\nabla V$. As discussed in Appendix~\ref{app:sinkhorn}, this can be solved to precision $\epsilon$ in time $O(m^2/\epsilon)$ using Sinkhorn's algorithm.

\section{Numerical Experiments}\label{sec:numerics}

\subsection{Simulated data}\label{sec:fig1_example}
In our first experiment, we compare the behaviour of the MFL dynamics to that of Global Waddington-OT (gWOT)~\cite{lavenant2021towards} in a setting with few samples per time-point. Although both methods minimize a functional of the form \eqref{eq:estimator}, gWOT works on a fixed, discretized support: the union of all observed sample points. Thus, gWOT should perform poorly when the support set has missing regions or ``gaps''. To simulate this, we simulated a dataset with $N_i$ samples at time $t_i$, with $N_1 = N_{10} = 64$ and $N_i = N, 2 \leq i \leq 9$ for $N \in \{ 2^0, \ldots, 2^6 \}$. The samples are drawn from the marginals of a bifurcating stochastic differential equation (see Appendix~\ref{sec:simdetails} for details).

In Figure \ref{fig:cover}(left) we illustrate two extreme cases: $N = 1$ (very few samples at intermediate timepoints) and $N = 64$ (uniform sampling over time) respectively, with $\lambda = 0.05, \lambda_{\mathrm{gWOT}} = 0.0025$.
Visually, it seems that the output of MFL is robust to the few-sample regime, with relatively little qualitative difference between the reconstructed trajectories for $N = 1, 64$. On the other hand, the performance of gWOT degrades visibly once the set of support points is a poor reflection of the underlying process.

To examine this quantitatively, we applied both MFL dynamics and gWOT for various values of $N$ and computed the root-mean-square (RMS) Energy Distance \cite{szekely2013energy} over time to an approximate ground truth (see Appendix~\ref{sec:simdetails}).
As shown in Figure \ref{fig:cover}(b), for large values of $N$ we observe that both MFL and gWOT perform similarly, but with MFL outperforming gWOT for small $N$ as expected.

Finally, Figure \ref{fig:dynamics} illustrates the evolution of the reduced objective $F$ \eqref{eq:reduced-objective} and the reconstructed marginals $\hat{\boldsymbol{\mu}}^{(i)}[k]$ over the course of MFL dynamics for $N = 64, \lambda = 0.05$. To estimate the entropy term $H(\boldsymbol{\mu})$, we used the nearest neighbour estimator of \cite[Section 6]{singh2016analysis}.

\subsection{Dealing with Branching}\label{sec:mass-variation}
In the context of single cell RNA-sequencing data it is crucial to be able to deal with the birth and death of cells, a.k.a.~branching, which does not occur homogeneously in the domain $\mathcal{X}$. Assume that we dispose of a function $g\in \mathcal{C}^1([0,1]\times \mathcal{X})$ such that $g(t,x)$ is the prior knowledge about the instantaneous growth rate of the distribution of particles in $x$ at time $t$ in the model of Eq.~\eqref{eq:SDE}. That is, in time $\d t$ the probability of a particle $X_t$ dividing is $g(t, X_t) \d t$~\cite{baradat2021regularized}.
We would like to incorporate this knowledge, and also allow for additional mass variations to account for the inaccuracy of our prior $g$.

To this end, we proceed heuristically by simply replacing $T_{\tau_i}(\mu_{t_i},\mu_{t_{i+1}})$ in Eq.~\eqref{eq:ROT} with the following ``unbalanced'' Schr\"odinger bridge problem
\begin{align}\label{eq:unbalanced}

    \min_{\gamma \in \mathcal{M}_+(\mathcal{X}^2)} \int c_{\tau_i}(x,y)\d\gamma(x,y) + \rho H(\gamma_1|\tilde \mu_{t_i}) + \rho H(\gamma_2|\tilde \mu_{t_{i+1}}) + \tau_i H(\gamma | \tilde \mu_{t_i}\otimes \tilde \mu_{t_{i+1}})
\end{align}
where $\tilde \mu_{t_i} \propto \exp(-g(x)(t_{i+1}-t_i)/2)\mu_{t_i}$ and $\tilde \mu_{t_{i+1}}\propto \exp(g(x)(t_{i+1}-t_i)/2)\mu_{t_{i+1}}$ are probability measures, $\rho > 0$ is a parameter, $\mathcal{M}_+(\mathcal{X}^2)$ is the set of nonnegative measures on $\mathcal{X}$ and $\gamma_1,\gamma_2$ are the marginals of $\gamma$. This problem can be solved with a variant of Sinkhorn's algorithm~\cite{chizat2018scaling,frogner2015learning}.

The rationale behind this formula is as follows: (i) the modifications $\tilde \mu_{t_{i}}$ of $\mu_{t_i}$ are intended to approximate the growth process $\partial_t \hat \mu_t = g(t,\cdot) \hat \mu_t$ over the time interval $[t_i,t_{i+1}]$ and (ii) we relax the marginal constraints, with a parameter $\rho>0$, to account for the potential inaccuracy of our prior $g$. We illustrate on Figure~\ref{fig:growth} the practical advantage of taking into account branching in the algorithm.

\begin{figure}

    \begin{subfigure}{0.49\linewidth}
        \includegraphics[width = \linewidth]{fig3_growth_a.pdf}
    \end{subfigure}
    \begin{subfigure}{0.49\linewidth}
        \includegraphics[width = \linewidth]{fig3_growth_b.pdf}
    \end{subfigure}
    \caption{Accounting for branching rates allows for the separation of dynamics from growth. (left) Ground truth paths and sample points. (right) Reconstruction produced by MFL dynamics with and without accounting for branching. Fraction of paths terminating in the upper (resp. lower) branch is annotated in the plot.}
    \label{fig:growth}
\end{figure}

\paragraph{Simulated data experiment} Here we consider a modified version of the bistable process from \cite[Figure 18]{lavenant2021towards} (see Appendix~\ref{sec:simdetails} for details).
Figure \ref{fig:growth} shows paths sampled from the ground truth, sample points, and the inferred trajectories. Since the lower potential well
is closer to the initial condition, we observe in the ground truth that $\sim 72\%$ of paths terminate in the lower branch. This effect is opposed by the high rate of branching in the upper branch, which ultimately results in more particles being sampled in the upper branch. By accounting for branching in the MFL dynamics we are able to reconstruct dynamics that isolate the effect of the potential, as evidenced by the comparable proportion (77\%) of inferred paths ending near the lower branch. On the other hand, neglecting the presence of branching results in growth effects being confounded with the potential, and the proportion of paths ending near each branch are roughly equal.

\subsection{Reprogramming dataset}

\begin{figure}
    \centering
    \begin{subfigure}{0.745\linewidth}
        \includegraphics[width = \linewidth]{reprogramming_snapshots.pdf}
    \end{subfigure}
    \begin{subfigure}{0.245\linewidth}
        \includegraphics[width = \linewidth]{reprogramming_distances.pdf}
    \end{subfigure}
    \caption{(left) Subsampled timepoints (Subsample) from the iPSC reprogramming dataset of \cite{schiebinger2019optimal}, with 100 cells at the first and last timepoints and 10 cells at each intermediate timepoint, shown beside inferred marginals obtained using MFL dynamics and the full dataset for those timepoints. (right) Best performance as measured by Energy Distance of inferred marginals to the full dataset for MFL dynamics (Langevin) and Global Waddington-OT~\cite{lavenant2021towards} (gWOT), with the subsampled timepoints (Sample) for reference.}
    \label{fig:reprogramming}
\end{figure}

We now consider the stem cell reprogramming dataset of \cite{schiebinger2019optimal}, in which a population of differentiating cells were profiled over a time course using single-cell RNA sequencing. For the purpose of this example, we restrict our attention to days 2.5-6 inclusive making a total of 8 timepoints. As previously, we consider a regime where few samples are observed as snapshots of the time-series and we apply MFL dynamics to infer reconstructed marginals. We reason that if the underlying process is well described by a Schr\"odinger bridge (as was empirically validated in \cite{schiebinger2019optimal}), then MFL dynamics should be able to ``improve'' on the sample marginals. As a proxy for ground truth, we used the full dataset (consisting of 59154 cells over the 8 timepoints). From this, we subsampled timepoints consisting of 100 cells at days 2.5 and 6, and 10 cells at the remaining intermediate timepoints. After carrying out a series of preprocessing steps (see Appendix~\ref{sec:preprocessing} for details), we applied MFL dynamics to the subsampled data to produce reconstructed marginals. Figure \ref{fig:reprogramming}(left) shows an example of the sampled points, MFL reconstruction, and the full dataset.

A reconstruction ``error'' was then computed as the Energy Distance \eqref{eq:energydistance} of each reconstructed marginal to the corresponding snapshot in the full dataset. Figure \ref{fig:reprogramming}(right) shows the error for each timepoint for each of MFL and gWOT, for the best performing parameters found in a parameter sweep. We observe that both MFL and gWOT improve upon the subsampled snapshots, but MFL by a larger margin.

\section{Conclusion}
We introduced a grid-free numerical method for trajectory inference that computes the min-entropy estimator introduced in~\cite{lavenant2021towards}, with global and quantitative convergence guarantees in the mean-field limit. This method arises naturally when decomposing the optimization problem in a suitable way and, in practice, outperforms the fixed-grid method of~\cite{lavenant2021towards}.

Concerning limitations, our method shares those of the min-entropy estimator : since it does not incorporate fine structural prior on the structure of the Waddington potential (i.e. it is fully nonparametric), it may suffer from a limited statistical efficiency in high dimension. An interesting research direction is to quantify this statistical performance and to adapt our algorithm to learn the Waddington potential $\Psi$ in a structured parameterized class of functions, jointly with the law on paths $P$.

On a more abstract level, our main insight is that min-entropy problems in Wiener space can be solved via a multi-species diffusion dynamics coupled with Schr\"odinger bridges. This point of view raises interesting theoretical questions, e.g.~can we rigorously interpret our dynamics as a diffusion in path space?

\section*{Acknowledgements}

GS was supported by a Career Award at the Scientific Interface from the Burroughs Wellcome Fund, a NFRF Exploration Grant, a NSERC Discovery Grant, and a CIHR Project Grant.

{
\small