# Efficient Optimal Transport Algorithm by Accelerated Gradient descent

Dongsheng An  Na Lei  Xianfeng Gu

###### Abstract

Optimal transport (OT) plays an essential role in various areas like machine
learning and deep learning. However, computing discrete optimal transport plan
for large scale problems with adequate accuracy and efficiency is still highly
challenging. Recently, methods based on the Sinkhorn algorithm add an entropy
regularizer to the prime problem and get a trade off between efficiency and
accuracy. In this paper, we propose a novel algorithm to further improve the
efficiency and accuracy based on Nesterov’s smoothing technique. Basically,
the non-smooth c-transform of the Kantorovich potential is approximated by the
smooth Log-Sum-Exp function, which finally smooths the original non-smooth
Kantorovich dual functional. The smooth Kantorovich functional can be
optimized by the fast proximal gradient algorithm (FISTA) efficiently.
Theoretically, the computational complexity of the proposed method is given by
$O(n^{\frac{5}{2}}\sqrt{\log n}/\epsilon)$, which is lower than current
estimation of the Sinkhorn algorithm. Empirically, compared with the Sinkhorn
algorithm, our experimental results demonstrate that the proposed method
achieves faster convergence and better accuracy with the same parameter.

##  1 Introduction

Optimal transport (OT) is a powerful tool to compute the Wasserstein distance
between probability measures, which are widely used to model various natural
and social phenomena, including economics (Galichon, 2016), optics (Glimm &
Oliker, 2003), biology (Schiebinger et al., 2019), physics (Jordan et al.,
1998) and so on. Recently, optimal transport has been successfully applied in
the areas of machine learning and statistics, such as parameter estimation in
Bayesian nonparametric models (Nguyen, 2013), computer vision (Arjovsky et
al., 2017; Courty et al., 2017; An et al., 2020a, b), natural language
processing (Kusner et al., 2015; Yurochkin et al., 2019) and so on. In these
areas, the complex probability measures are approximated by summations of
Dirac measures supported on the samples. To obtain the Wasserstein distance
between the empirical distributions, we have to solve the discrete OT problem.

For the discrete optimal transport problem, where both the source and target
measures are discrete, the Kantorovich functional becomes a convex function
defined on a convex domain. Due to the lack of smoothness, conventional
gradient descend method can not be applied directly. Instead, it can be
optimized with the sub-differential method (Nesterov, 2005), where the
gradient is replaced by the sub-differential. In order to achieve an
approximation error less than $\varepsilon$, the sub-differential method
requires $O(1/\varepsilon^{2})$ iterations. Recently, several approximation
methods have been proposed to improve the computational efficiency. In these
methods (Cuturi, 2013; Benamou et al., 2015; Altschuler et al., 2017), a
strongly convex entropy function is added to the prime Kantorovich problem and
thus the regularized problem can be efficiently solved by the Sinkhorn
algorithm. More detailed analysis shows that the computational complexity of
the Sinkhorn algorithm is $\tilde{O}(n^{2}/\varepsilon^{2})$ (Dvurechensky et
al., 2018) by setting $\lambda=\epsilon/4\log n$. Also, a series of primal-
dual algorithms are proposed, including the APDAGD algorithm (Dvurechensky et
al., 2018) with computational complexity $\tilde{O}(n^{2.5}/\epsilon)$, the
APDAMD algorithm (Lin et al., 2019) with $\tilde{O}(n^{2}\sqrt{r}/\epsilon)$
where $r$ is a complex constant of the Bregman divergence, and the APDRCD
algorithm (Guo et al., 2020) with $\tilde{O}(n^{2.5}/\epsilon)$. But all of
the three methods need to build a matrix with space complexity $O(n^{3})$,
which makes them hard to use when $n$ is large.

Our Method In this work, instead of starting from the prime Kantorovich
problem like the Sinkhorn based methods, we directly deal with the dual
Kantorovich problem. The key idea is to approximate the original non-smooth
c-transform of the Kantorovich potential by Nesterov’s smoothing idea.
Specifically, we approximate the $\max$ function by the Log-Sum-Exp function,
which has also been used in (Schmitzer, 2019; Peyré & Cuturi, 2018), such that
the original non-smooth Kantorovich functional is converted to be an
unconstrained $(n-1)$-dimensional smooth convex energy. By using the Fast
Proximal Gradient Method (FISTA) (Beck & Teboulle, 2008), we can quickly
optimize the smoothed energy to get a precise estimate of the OT cost. In
theory, the method can achieve the approximate error $\varepsilon$ with space
complexity $O(n^{2})$ and computational complexity ${O}(n^{2.5}\sqrt{\log
n}/\varepsilon)$.

The contributions of the proposed method includes: (i) We convert the dual
Kantorovich problem to be an unconstrained convex optimization problem by
approximating the nonsmooth c-transform of the Kantorovich potential with
Nesterov’s smoothing. (ii) The smoothed Kantorovich functional can be
efficiently solved by the FISTA algorithm with computational complexity
$\tilde{O}(n^{2.5}/{\sqrt{\varepsilon}})$. At the same time, the computational
complexity of the Kantorovich functional itself is given by
$\tilde{O}(n^{2.5}/\varepsilon)$. (iii) The experiments demonstrate that
compared with the Sinkhorn algorithm, the proposed method achieves faster
convergence, better accuracy and higher stability with the same parameter
$\lambda$.

##  2 Optimal Transport Theory

In this section, we introduce some basic concepts and theorems in the
classical OT theory, focusing on Kantorovich’s approach and its generalization
to the discrete settings via c-transform. The details can be found in
Villani’s book (Villani, 2008).

######  Problem 1 (Kantorovich Problem).

Suppose $X\subset\mathbb{R}^{d}$, $Y\subset\mathbb{R}^{d}$ are two subsets of
the Euclidean space $\mathbb{R}^{d}$,
$\mu=\sum_{i=1}^{m}\mu_{i}\delta(x-x_{i}),\nu=\sum_{j=1}^{n}\nu_{j}\delta(y-y_{j})$
are two discrete probability measures defined on $X$ and $Y$ with equal total
measure, $\sum_{i=1}^{m}\mu_{i}=\sum_{j=1}^{n}\nu_{j}$. Then the Kantorovich
problem aims to find the optimal transport plan $P=(p_{ij})$ that minimizes
the total transport cost:

$\small(KP)\hskip
5.69054ptM_{c}(\mu,\nu)=\min_{P\in\pi(\mu,\nu)}\sum_{i=1}^{m}\sum_{j=1}^{n}p_{ij}c_{ij}.$
(1)

where $c_{ij}=c(x_{i},y_{j})$ is the cost that transports one unit mass from
$x_{i}$ to $y_{j}$, and
$\pi(\mu,\nu)=\\{P|\sum_{i=1}^{m}p_{ij}=\nu_{j},\sum_{j=1}^{n}p_{ij}=\mu_{i},p_{ij}\geq
0\\}$.

######  Problem 2 (Dual Kantorovich).

Given two discrete probability measures
$\mu=\sum_{i=1}^{m}\mu_{i}\delta(x-x_{i}),\nu=\sum_{j=1}^{n}\nu_{j}\delta(y-y_{j})$
supported on $X$ and $Y$, and the transport cost function $c:X\times
Y\to\mathbb{R}$, the Kantorovich problem is equivalent to maximizing the
following _Kantorovich functional_ :

$\begin{split}\small(DP1)~{}M_{c}(\mu,\nu)=\max_{\phi,\psi}\left\\{-\sum_{i=1}^{m}\phi_{i}\mu_{i}+\sum_{j=1}^{n}\psi_{j}\nu_{j}\right\\}\end{split}$
(2)

where $\phi\in L^{1}(X,\mu)$ and $\psi\in L^{1}(Y,\nu)$ are called the
_Kantorovich potentials_ , and $-\phi_{i}+\psi_{j}\leq c_{ij}$. The above
problem can be reformulated as the following minimization form with the same
constraints:

$\small(DP2)~{}M_{c}(\mu,\nu)=-\min\left\\{\sum_{i=1}^{m}\phi_{i}\mu_{i}-\sum_{j=1}^{n}\psi_{j}\nu_{j}\right\\}$
(3)

######  Definition 3 (c-transform).

Let $\phi\in L^{1}(X,\mu)$ and $\psi\in L^{1}(Y,\nu)$, we define

$\small\phi_{i}=\psi^{c}(x_{i})=\max_{j}\\{\psi_{j}-c_{ij}\\}$ (4)

With c-transform, Eqn. (3) is equivalent to solving the following
unconstrained convex optimization problem:

$\small
M_{c}(\mu,\nu)=-\min_{\psi}E(\psi)=-\min_{\psi}\\{\sum_{i=1}^{m}\mu_{i}\psi^{c}(x_{i})-\sum_{j=1}^{n}\nu_{j}\psi_{j}\\}$
(5)

Suppose $\psi^{*}$ is the solution to the problem Eqn. (5), then
$\psi^{*}+k\mathbf{1}$ is also an optimal solution for Eqn. (5). In order to
make the solution unique, we add a constraint $\psi\in H$ using the indicator
function $I_{H}$, where $H=\\{\psi|\sum_{j=1}^{n}\psi_{j}=0\\}$, and modify
the Kantorovich functional $E(\psi)$ in Eqn. (5):

$\small\tilde{E}(\psi)=E(\psi)+I_{H}(\psi),\quad
I_{H}(\psi)=\left\\{\begin{array}[]{ll}0&\psi\in H\\\ \infty&\psi\not\in
H\end{array}\right.$ (6)

Then the problem Eqn. (5) is equivalent to solving:

$\small\begin{split}M_{c}(\mu,\nu)&=-\min_{\psi}\tilde{E}(\psi)\\\
&=-\min_{\psi}\sum_{i=1}^{m}\mu_{i}\psi^{c}(x_{i})-\sum_{j=1}^{n}\nu_{j}\psi_{j}+I_{H}(\psi)\end{split}$
(7)

which is essentially an $(n-1)$-dimensional unconstrained convex problem.
According to the definition of $c$-transform in Eqn. (4), $\psi^{c}$ is non-
smooth w.r.t $\psi$.

##  3 Nesterov’s Smoothing of Kantorovich functional

In this section, we smooth the non-smooth discrete Kantorovich functional
$E(\psi)$ by approximating $\psi^{c}(x)$ with the Log-Sum-Exp function to get
the smooth Kantorovich functional $E_{\lambda}(\psi)$, following Nesterov’s
original strategy (Nesterov, 2005), which has also been applied in the OT
field (Peyré & Cuturi, 2018; Schmitzer, 2019). Then through the FISTA
algorithm (Beck & Teboulle, 2008), we can easily induce that the computation
complexity of our algorithm is $O(n^{2.5}\sqrt{\log n}/\varepsilon)$, with
$\tilde{E}(\psi^{t})-\tilde{E}(\psi^{*})\leq\epsilon$. By abuse of notation,
in the following we call both $E(\psi)$ and $\tilde{E}(\psi)$ the Kantorovich
functional, both $E_{\lambda}(\psi)$ and $\tilde{E}_{\lambda}(\psi)$ the
smooth Kantorovich functional.

######  Definition 4 ($(\alpha,\beta)$-smoothable).

A convex function $f$ is called $(\alpha,\beta)$-smoothable if for any
$\lambda>0$, $\exists$ a convex function $f_{\lambda}$ s.t.

$\small\begin{split}f_{\lambda}(x)&\leq f(x)\leq
f_{\lambda}(x)+\beta\lambda\\\ f_{\lambda}(y)&\leq
f_{\lambda}(x)+\langle\nabla
f_{\lambda}(x),y-x\rangle+\frac{\alpha}{2\lambda}(y-x)^{T}\nabla^{2}f_{\lambda}(x)(y-x)\end{split}$

Here $f_{\lambda}$ is called a $\frac{1}{\lambda}$-smooth approximation of $f$
with parameters $(\alpha,\beta)$.

In the above definition, the parameter $\lambda$ defines a tradeoff between
the approximation accuracy and the smoothness, the smaller the $\lambda$, the
better approximation and the less smoothness.

######  Lemma 5 (Nesterov’s Smoothing).

Given $f:\mathbb{R}^{n}\to\mathbb{R}$, $f(x)=\max\\{x_{j}:j=1,\dots,n\\}$, for
any $\lambda>0$, we have its $\frac{1}{\lambda}$-smooth approximation with
parameters $(1,\log n)$

$\small
f_{\lambda}(x)=\lambda\log\left(\sum_{j=1}^{n}e^{{x_{j}}/{\lambda}}\right)-\lambda\log
n,$ (8)

Recall the definition of c-transform of the Kantorovich potential in Eqn. (4),
we obtain the Nesterov’s smoothing of $\psi^{c}$ by applying Eqn. (8)

$\small\psi^{c}_{\lambda}(x_{i})=\lambda\log\left(\sum_{j=1}^{n}e^{(\psi_{j}-c_{ij})/\lambda}\right)-\lambda\log
n.$ (9)

We use $\psi_{\lambda}^{c}$ to replace $\psi^{c}$ in Eqn. (7) to approximate
the Kantorovich functional. Then the Nesterov’s smoothing of the Kantorovich
functional becomes

$\small
E_{\lambda}(\psi)=\lambda\sum_{i=1}^{m}\mu_{i}\log\left(\sum_{j=1}^{n}e^{(\psi_{j}-c_{ij})/\lambda}\right)-\sum_{j=1}^{n}\nu_{j}\psi_{j}-\lambda\log
n$ (10)

and its gradient is given by

$\small\frac{\partial
E_{\lambda}(\psi)}{\partial\psi_{j}}=\sum_{i=1}^{m}\mu_{i}\frac{e^{(\psi_{j}-c_{ij})/\lambda}}{\sum_{k=1}^{n}e^{(\psi_{k}-c_{ik})/\lambda}}-\nu_{j},~{}~{}\forall~{}j\in[n]$
(11)

Furthermore, we can directly compute the Hessian matrix of
$E_{\lambda}(\psi)$. Let $K_{ij}=e^{-c_{ij}/\lambda}$ and
$v_{j}=e^{\psi_{j}/\lambda}$, and set
$E_{\lambda}^{i}:=\lambda\log\sum_{j=1}^{n}K_{ij}v_{j},~{}\forall~{}i\in[m]$.
Direct computation gives the following Hessian matrix:

$\small\begin{split}\nabla^{2}E_{\lambda}&=\frac{1}{\lambda}\sum_{i=1}^{m}\mu_{i}\left(\frac{1}{\mathbf{1}^{T}V_{i}}\Lambda_{i}-\frac{1}{(\mathbf{1}^{T}V_{i})^{2}}V_{i}V_{i}^{T}\right)\end{split}$
(12)

where $V_{i}=(K_{i1}v_{1},K_{i2}v_{2},\dots,K_{in}v_{n})^{T}$, and
$\Lambda_{i}=\text{diag}(K_{i1}v_{1},K_{i2}v_{2},\dots,K_{in}v_{n})$. By the
Hessian matrix, we can show that $E_{\lambda}$ is a smooth approximation of
$E$.

######  Lemma 6.

$E_{\lambda}(\psi)$ is a $\frac{1}{\lambda}$-smooth approximation of $E(\psi)$
with parameters $(1,\log n)$.

######  Lemma 7.

Suppose $E_{\lambda}(\psi)$ is the $\frac{1}{\lambda}$-smooth approximation of
$E(\psi)$ with parameters $(1,\log n)$, $\psi_{\lambda}^{*}$ is the optimizer
of $E_{\lambda}(\psi)$, the approximate OT plan is unique and given by

$\small\begin{split}(P_{\lambda}^{*})_{ij}=\mu_{i}\frac{e^{((\psi_{\lambda}^{*})_{j}-c_{ij})/\lambda}}{\sum_{k=1}^{n}e^{((\psi_{\lambda}^{*})_{k}-c_{ik})/\lambda}}=\frac{\mu_{i}K_{ij}v_{j}^{*}}{K_{i}v^{*}}\end{split}$
(13)

where $K_{i}$ is the $ith$ row of $K$ and
$v^{*}=e^{\psi_{\lambda}^{*}/\lambda}$.

Similar to the discrete Kantorovich functional in Eqn. (5), the optimizer of
the smooth Kantorovich functional in Eqn. (10) is also not unique: given an
optimizer $\psi_{\lambda}^{*}$, then $\psi_{\lambda}^{*}+k\mathbf{1}$,
$k\in\mathbb{R}$ is also an optimizer. We can eliminate the ambiguity by
adding the indicator function as Eqn. (6),
$\tilde{E}_{\lambda}(\psi)=E_{\lambda}(\psi)+I_{H}(\psi)$,

$\small\begin{split}\tilde{E}_{\lambda}(\psi)&=\lambda\sum_{i=1}^{m}\mu_{i}\log\left(\sum_{j=1}^{n}e^{(\psi_{j}-c_{ij})/\lambda}\right)-\sum_{j=1}^{n}\nu_{j}\psi_{j}-\lambda\log
n+I_{H}(\psi)\end{split}$ (14)

This energy can be optimized effectively through the accelerated proximal
algorithms (Nesterov, 2005), such as the Fast Proximal Gradient Method (FISTA)
(Beck & Teboulle, 2008). The FISTA iterations are as follows:

$\small\begin{split}z^{t+1}&=\Pi_{\eta_{t}I_{H}}\left(\psi^{t}-\eta_{t}\nabla
E_{\lambda}(\psi^{t})\right)\\\
\psi^{t+1}&=z^{t+1}+\frac{\theta_{t}-1}{\theta_{t+1}}(z^{t+1}-z^{t})\end{split}$
(15)

with initial conditions $\psi^{0}=v^{0}=\mathbf{0}$, $\theta_{0}=1$,
$\eta_{t}=\lambda$ and
$\theta_{t+1}=\frac{1}{2}\left(1+\sqrt{1+4\theta_{t}^{2}}\right)$. Here
$\Pi_{\eta_{t}I_{H}}(z)=z-\frac{1}{n}\sum_{j=1}^{n}z_{j}$ is the projection of
$z$ to $H$ (the proximal function of $I_{H}(x)$) (Parikh & Boyd, 2014).
Similar to the Sinkhorn’s algorithm, this algorithm can be parallelized, since
all the operations are row based.

######  Theorem 8.

Given the cost matrix $C=(c_{ij})$, the source measure $\mu$ and target
measure $\nu$ with $\sum_{i=1}^{m}\mu_{i}=\sum_{j=1}^{n}\nu_{j}=1$, $\psi^{*}$
is the optimizer of the discrete dual Kantorovich functional
$\tilde{E}(\psi)$, and $\psi_{\lambda}^{*}$ is the optimizer of the smooth
Kantorovich functional $\tilde{E}_{\lambda}(\psi)$. Then the approximation
error is

$\small|\tilde{E}(\psi^{*})-\tilde{E}_{\lambda}(\psi_{\lambda}^{*})|\leq
2\lambda\log n$

######  Corollary 9.

Suppose $\lambda$ is fixed and $\psi^{0}=\mathbf{0}$, then for any
$t\geq\sqrt{\frac{{2\|\psi_{\lambda}^{*}\|^{2}}}{\lambda\varepsilon}}$, we
have

$\small\tilde{E}_{\lambda}(\psi^{t})-\tilde{E}_{\lambda}(\psi_{\lambda}^{*})\leq\varepsilon,$

where $\psi_{\lambda}^{*}$ is the optimizer of $\tilde{E}_{\lambda}(\psi)$.

######  Theorem 10.

If $\lambda=\frac{\varepsilon}{2\log n}$, then for any
$t\geq\frac{\sqrt{8\|\bar{C}\|^{2}n\log n}}{\epsilon}$ with
$\bar{C}=C_{\max}-\lambda\log\nu_{min}$, we have

$\small\tilde{E}(\psi^{t})-\tilde{E}(\psi^{*})<\varepsilon$

where $\psi^{t}$ is the solver of $\tilde{E}_{\lambda}(\psi)$ after $t$ steps
in the iterations in Eqn. (15), and $\psi^{*}$ is the optimizer of
$\tilde{E}(\psi)$. Then the total computational complexity is
$O(\frac{n^{2.5}\sqrt{\log n}}{\varepsilon})$.

Since $\tilde{E}(\psi^{t})>\tilde{E}(\psi_{\lambda}^{*})$, thus we have

$\small 0<\tilde{E}(\psi_{\lambda}^{*})-\tilde{E}(\psi^{*})<2\lambda\log n$
(16)

This shows that $\tilde{E}(\psi_{\lambda}^{*})$ at least linearly converges to
$\tilde{E}(\psi^{*})$ with respect to $\lambda$.

|  |  |
---|---|---|---
Figure 1: Comparison with the Sinkhorn algorithm (Cuturi, 2013) under
different costs.

##  4 Experiments

In this section, we investigate the performance of the proposed algorithm by
comparing with the Sinkhorn algorithm (Cuturi, 2013) with computational
complexity $\tilde{O}(n^{2}/\epsilon^{2})$ (Dvurechensky et al., 2018). All of
the codes are written using MATLAB with GPU acceleration. The experiments are
conducted on a Windows laptop with Intel Core i7-7700HQ CPU, 16 GB memory and
NVIDIA GTX 1060Ti GPU.

Cost Matrix We test the performance of the algorithm with different cost
matrix. Specifically, we set $\mu=\sum_{i=1}^{m}\mu_{i}\delta(x-x_{i})$,
$\nu=\sum_{j=1}^{n}\nu_{j}\delta(y-y_{j})$. After the settings of $\mu_{i}$s
and $\nu_{j}$s, they are normalized by
$\mu_{i}=\frac{\mu_{i}}{\sum_{k=1}^{m}\mu_{k}}$,
$\nu_{j}=\frac{\nu_{j}}{\sum_{k=1}^{n}\nu_{k}}$. To build the cost matrix, we
use the squared Euclidean distance and the spherical distance.

  * •

For the squared Euclidean distance (SED) experiment, like (Altschuler et al.,
2017), we randomly choose one pair of images from the MNIST dataset (LeCun &
Cortes, 2010), and then add negligible noise $0.01$ to each background pixel
with intensity 0. $\mu_{i}$ and $x_{i}$ ($\nu_{j}$ and $y_{j}$) are set to be
the value and the coordinate of each pixel in the source (target) image. Then
the squared Euclidean distance between $x_{i}$ and $y_{j}$ are given by
$c(x_{i},y_{j})=\|x_{i}-y_{j}\|^{2}$.

  * •

For the spherical distance (SD) experiment, we set $m=n=500$. Both $\mu_{i}$s
and $\nu_{j}$s are randomly generated from the uniform distribution
$Uni([0,1])$. $x_{i}$’s are randomly sampled from the Gaussian distribution
$\mathbb{N}(3\mathbf{1_{d}},I_{d})$ and $y_{j}$’s are randomly sampled from
the Uniform distribution $Uni([0,1]^{d})$. Then we normalize $x_{i}$ and
$y_{j}$ by $x_{i}=\frac{x_{i}}{\|x_{i}\|_{2}}$ and
$y_{j}=\frac{y_{j}}{\|y_{j}\|_{2}}$. As a result, both $x_{i}$’s and $y_{j}$’s
are located on the sphere. The spherical distance is given by
$c(x_{i},y_{j})=\arccos{(\langle x_{i},y_{j}\rangle)}$.

Evaluation Metrics W mainly use two metrics to evaluate the proposed method:
the first is the total transport cost, which is defined by Eqn. (5) and given
by $-E(\psi)$; and the second is the $L_{1}$ distance from the computed
transport plan $P_{\lambda}$ to the admissible distribution space
$\pi(\mu,\nu)$ defined as
$D(P_{\lambda})=\|P_{\lambda}\mathbf{1}-\mu\|_{1}+\|P_{\lambda}^{T}\mathbf{1}-\nu\|_{1}$.

Then we compare with the Sinkhorn algorithm (Cuturi, 2013) w.r.t both the
convergence rate and the approximate accuracy. We manually set
$\lambda=\frac{\max(C)-\min(C)}{T}$ with $T=500$ to get a good estimate of the
OT cost, and adjust the step length $\eta$ to get the best convergence of the
proposed method. For fair comparison, we use the same $\lambda$ for the
Sinkhorn algorithm. We summarize the results in Fig. 1, where the green curves
represent the groundtruth computed by linear programming, the blue curves are
for the Sinkhorn algorithm, and the red curves give the results of our method.

From the figure, we can observe that in both experiments, our method achieves
faster convergence than the Sinkhorn algorithm. Also, $-E(\psi_{\lambda}^{*})$
gives a better approximate of the OT cost than $\langle P_{\lambda},C\rangle$
used in the Sinkhorn algorithm with the same small $\lambda$, especially for
the squared euclidean distance cost. Basically, to achieve $\epsilon$
precision, $\langle P_{\lambda}^{*},C\rangle$ needs to set
$\lambda=\frac{\epsilon}{4\log n}$ (Dvurechensky et al., 2018), which is
smaller than our requirement for $\lambda=\frac{\epsilon}{2\log n}$.

We also report the iterations and the corresponding running time with $T=700$
in Tab. 1. The stop condition is set to be
$|E(\psi^{t+1})-E(\psi^{t})|/|E(\psi^{t})|<10^{-3}$. In each iteration,
Sinkhorn requires $2n^{2}$ times of operation, and the proposed method
requires $3n^{2}$ times of operation. From Tab. 1 we can see that the proposed
method needs less iterations and time for convergence.

Table 1: The running iterations and time comparison between the proposed method and Sinkhorn (Cuturi, 2013) with $T=700$. | Sinkhorn | Ours
---|---|---
Iterations | Time(s) | Iterations | Time(s)
SED | 209 | 0.0431 | 29 | 0.0197
SD | 297 | 0.0564 | 22 | 0.0142

##  5 Conclusion

In this paper, we propose a novel algorithm to improve the accuracy for
solving the discrete OT problem based on Nesterov’s smoothing technique. The
c-transform of the Kantorovich potential is approximated by the smooth Log-
Sum-Exp function, and the smoothed Kantorovich functional can be solved by the
fast proximal gradient algorithm (FISTA) efficiently. Theoretically, the
computational complexity of the proposed method is given by
$O(n^{2.5}\sqrt{\log n}/\varepsilon)$, which is lower than current estimation
of the Sinkhorn method. Empirically, our experimental results demonstrate that
the proposed method achieves faster convergence and better accuracy than the
Sinkhorn algorithm.

## References

  * Altschuler et al. (2017) Altschuler, J., Niles-Weed, J., and Rigollet, P.  Near-linear time approximation algorithms for optimal transport via sinkhorn iteration.  In _Advances in Neural Information Processing Systems 30_. 2017.
  * An et al. (2020a) An, D., Guo, Y., Lei, N., Luo, Z., Yau, S.-T., and Gu, X.  Ae-ot: A new generative model based on extended semi-discrete optimal transport.  In _International Conference on Learning Representations_ , 2020a.
  * An et al. (2020b) An, D., Guo, Y., Zhang, M., Qi, X., Lei, N., and Gu, X.  AE-OT-GAN: Training GANs from data specific latent distribution.  In _European Conference on Computer Vision (ECCV)_ , pp. 548–564, 2020b.
  * Arjovsky et al. (2017) Arjovsky, M., Chintala, S., and Bottou, L.  Wasserstein generative adversarial networks.  In _ICML_ , pp. 214–223, 2017.
  * Beck & Teboulle (2008) Beck, A. and Teboulle, M.  A fast iterative shrinkage-thresholding algorithm for linear inverse problems.  _SIAM J. IMAGINGSCIENCES_ , 2008.
  * Benamou et al. (2015) Benamou, J.-D., Carlier, G., Cuturi, M., Nenna, L., and Peyré, G.  Iterative bregman projections for regularized transportation problems.  _SIAM Journal on Scientific Computing_ , 2015.
  * Courty et al. (2017) Courty, N., Flamary, R., Tuia, D., and Rakotomamonjy, A.  Optimal transport for domain adaptation.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 39(9):1853–1865, 2017.
  * Cuturi (2013) Cuturi, M.  Sinkhorn distances: Lightspeed computation of optimal transportation distances.  In _International Conference on Neural Information Processing Systems_ , 2013.
  * Dvurechensky et al. (2018) Dvurechensky, P., Gasnikov, A., and Kroshnin, A.  Computational optimal transport: Complexity by accelerated gradient descent is better than by sinkhorn’s algorithm.  In _Proceedings of the 35th International Conference on Machine Learning_. PMLR, 2018.
  * Galichon (2016) Galichon, A.  _Optimal Transport Methods in Economics_.  Princeton University Press, 2016.
  * Glimm & Oliker (2003) Glimm, T. and Oliker, V.  Optical design of single reflector systems and the monge–kantorovich mass transfer problem.  _Journal of Mathematical Sciences_ , 117(3):4096–4108, Sep 2003.
  * Guo et al. (2020) Guo, W., Ho, N., and Jordan, M.  Fast algorithms for computational optimal transport and wasserstein barycenter.  In _International Conference on Artificial Intelligence and Statistics (AISTATS)_ , 2020.
  * Horn & Johnson (1991) Horn, T. A. and Johnson, C. R.  _Topics in Matrix Analysis_.  Cambridge, 1991.
  * Jordan et al. (1998) Jordan, R., Kinderlehrer, D., and Otto, F.  The variational formulation of the fokker–planck equation.  _SIAM Journal on Mathematical Analysis_ , 29(1):1–17, 1998.
  * Kusner et al. (2015) Kusner, M., Sun, Y., Kolkin, N., and Weinberger, K.  From word embeddings to document distances.  In _Proceedings of the 32nd International Conference on Machine Learning_ , pp. 957–966, 2015.
  * LeCun & Cortes (2010) LeCun, Y. and Cortes, C.  MNIST handwritten digit database.  2010\.  URL http://yann.lecun.com/exdb/mnist/.
  * Lin et al. (2019) Lin, T., Ho, N., and Jordan, M.  On efficient optimal transport: An analysis of greedy and accelerated mirror descent algorithms.  In _International Conference on Machine Learning_ , 2019.
  * Nesterov (2005) Nesterov, Y.  Smooth minimization of non-smooth functions.  _Mathematical programming_ , 2005.
  * Nguyen (2013) Nguyen, X.  Convergence of latent mixing measures in finite and infinite mixture models.  _Ann. Statist_ , 41, 2013.
  * Parikh & Boyd (2014) Parikh, N. and Boyd, S.  _Proximal Algorithms_.  Foundations and Trends in Optimization, 2014.
  * Peyré & Cuturi (2018) Peyré, G. and Cuturi, M.  _Computational Optimal Transport_.  https://arxiv.org/abs/1803.00567, 2018.
  * Schiebinger et al. (2019) Schiebinger, G., Shu, J., Tabaka, M., Cleary, B., Subramanian, V., Solomon, A., Gould, J., Liu, S., Lin, S., Berube, P., Lee, L., Chen, J., Brumbaugh, J., Rigollet, P., Hochedlinger, K., Jaenisch, R., Regev, A., and Lander, E.  Optimal-transport analysis of single-cell gene expression identifies developmental trajectories in reprogramming.  _Cell_ , 2019.
  * Schmitzer (2019) Schmitzer, B.  Stabilized sparse scaling algorithms for entropy regularized transport problems.  _SIAM Journal on Scientific Computing_ , 2019.
  * Villani (2008) Villani, C.  _Optimal transport: old and new_ , volume 338.  Springer Science & Business Media, 2008.
  * Yurochkin et al. (2019) Yurochkin, M., Claici, S., Chien, E., Mirzazadeh, F., and Solomon, J. M.  Hierarchical optimal transport for document representation.  In _Advances in Neural Information Processing Systems 32_. 2019.

##  Appendix A The Proofs

In this section, we give the proofs of the lemmas, theorems and propositions
in the main paper.

###### Proof of Lemma 5.

We have $\forall x\in\mathbb{R}^{n}$,

$\begin{split}f_{\lambda}(x)&\leq\lambda\log(n\max_{j}e^{x_{j}/\lambda})-\lambda\log(n)=f(x)\\\
f(x)&=\lambda\log\max_{j}e^{x_{j}/\lambda}<\lambda\log\sum_{j=1}^{n}e^{x_{j}/\lambda}=f_{\lambda}(x)+\lambda\log
n\end{split}$

Furthermore, it is easy to prove that $f_{\lambda}(x)$ is
$\frac{1}{\lambda}$-smooth. Therefore, $f_{\lambda}(x)$ is an approximation of
$f(x)$ with parameters $(1,\log n)$. ∎

###### Proof of Lemma 6.

From Eqn. (12), we see
$\nabla^{2}E_{\lambda}^{i}=\frac{1}{\lambda}\left(\frac{1}{\mathbf{1}^{T}V_{i}}\Lambda_{i}-\frac{1}{(\mathbf{1}^{T}V_{i})^{2}}V_{i}V_{i}^{T}\right)$
has $O=\\{k\mathbf{1}:k\in\mathbb{R}\\}$ as its null space. In the orthogonal
complementary space of $O$, $\nabla^{2}E_{\lambda}^{i}$ is diagonal dominant,
therefore strictly positive definite.

Weyl’s inequality (Horn & Johnson, 1991) states that the eigen value of
$A=B+C$ is no greater than the maximal eigenvalue of $B$ minus the minimal
eigenvalue of $C$, where $B$ is an exact matrix and $C$ is a perturbation
matrix. Hence the maximal eigenvalue of $\nabla^{2}E_{\lambda}^{i}$, denoted
as $\lambda_{i}$, has an upper bound,

$0\leq\lambda_{i}\leq\frac{1}{\lambda}\frac{1}{\mathbf{1}^{T}V_{i}}\max_{j}\\{K_{ij}v_{j}\\}\leq\frac{1}{\lambda}.$

Thus the maximal eigenvalue of $\nabla^{2}E_{\lambda}(\psi)$ is no greater
than $\sum_{i=1}^{n}\mu_{i}\lambda_{i}=\frac{1}{\lambda}$. It is easy to find
that $E_{\lambda}(\psi)\leq E(\psi)\leq E_{\lambda}(\psi)+\lambda\log n$.
Thus,$E_{\lambda}(\psi)$ is a $\frac{1}{\lambda}$-smooth approximation of
$E(\psi)$ with parameters $(1,\log n)$. ∎

###### Proof of Lemma 7.

By the gradient formula Eqn. (11) and the optimizer $\psi_{\lambda}^{*}$, we
have

$\frac{\partial
E_{\lambda}(\psi_{\lambda}^{*})}{\partial\psi_{j}}=\sum_{i=1}^{m}(P_{\lambda}^{*})_{ij}-\nu_{j}=0\quad\forall~{}j=1,\cdots,n.$

On the other hand, by the definition of $P_{\lambda}^{*}$, we have

$\sum_{j=1}^{n}(P_{\lambda}^{*})_{ij}=\mu_{i},\quad\forall~{}i=1,\dots,m,$

compare the above two equations, we obtain that
$P_{\lambda}^{*}\in\pi(\mu,\nu)$ is the approximate OT plan.

The optimal solvers of $E_{\lambda}(\psi)$ has the form
$\psi_{\lambda}^{*}+\mathbf{1}$, all of them induce the same approximate
transport plan given by Eqn. (13). Thus, $P_{\lambda}^{*}$ is unique.

∎

###### Proof of Theorem 8.

Assume $\psi^{*}$ and $\psi_{\lambda}^{*}$ are the minimizers of $E(\psi)$ and
$E_{\lambda}(\psi)$ respectively. Then by the inequality in Eqn. (A)

$\begin{array}[]{llll}E_{\lambda}(\psi^{*})&\leq E(\psi^{*})&\leq
E(\psi_{\lambda}^{*})&\leq E_{\lambda}(\psi_{\lambda}^{*})+\lambda\log n\\\
E_{\lambda}(\psi_{\lambda}^{*})&\leq E_{\lambda}(\psi^{*})&\leq
E(\psi^{*})&\leq E_{\lambda}(\psi^{*})+\lambda\log n\end{array}$

This shows
$|E_{\lambda}(\psi^{*})-E_{\lambda}(\psi_{\lambda}^{*})|\leq\lambda\log n$.
Removing the indicator functions, we can get

$\begin{split}|\tilde{E}(\psi^{*})-\tilde{E}_{\lambda}(\psi_{\lambda}^{*})|&=|E(\psi^{*})-E_{\lambda}(\psi_{\lambda}^{*})|\\\
&\leq|E(\psi^{*})-E_{\lambda}(\psi^{*})|+|E_{\lambda}(\psi^{*})-E_{\lambda}(\psi_{\lambda}^{*})|\\\
&\leq 2\lambda\log n.\end{split}$

∎

###### Proof of Corollary 9.

Under the settings of the smoothed Kantorovich problem Eqn. (10),
$E_{\lambda}(\psi)$ is convex and differentiable with
$\nabla^{2}E_{\lambda}(\psi)\preceq\frac{1}{\lambda}I$, $I_{H}(\psi)$ is
convex and its proximal function is given by $\Pi_{H}(v)$. Thus, directly
applying Thm 4.4 of (Beck & Teboulle, 2008) and set $L=\frac{1}{\lambda}$, we
can get
$\tilde{E}_{\lambda}(\psi^{t})-\tilde{E}_{\lambda}(\psi_{\lambda}^{*})\leq\frac{2}{\lambda(t+1)^{2}}\|\psi_{\lambda}^{*}-\psi^{0}\|^{2}$.
Set $\psi^{0}=\textbf{0}$ and
$\frac{2}{\lambda(t+1)^{2}}\|\psi_{\lambda}^{*}\|^{2}\leq\varepsilon$, then we
get that when
$t\geq\sqrt{\frac{2\|\psi_{\lambda}^{*}\|^{2}}{\lambda\varepsilon}}$, we have
$\tilde{E}_{\lambda}(\psi^{t})-\tilde{E}_{\lambda}(\psi_{\lambda}^{*})\leq\varepsilon$.
∎

###### Proof of Theorem 10.

We set the initial condition $\psi^{0}=\mathbf{0}$. For any given
$\varepsilon>0$, we choose iteration step $t$, such that
$\frac{2}{\lambda(t+1)^{2}}\|\psi_{\lambda}^{*}\|^{2}\leq\frac{\varepsilon}{2}$,
$t\geq\frac{\sqrt{8\|\psi_{\lambda}^{*}\|^{2}\log n}}{\epsilon}$, where
$\psi_{\lambda}^{*}$ is the optimizer of $\tilde{E}_{\lambda}(\psi)$. By Thm
4.4 of (Beck & Teboulle, 2008), we have

$\begin{split}\tilde{E}_{\lambda}(\psi^{t})-\tilde{E}(\psi^{*})&\leq\tilde{E}_{\lambda}(\psi^{t})-\tilde{E}_{\lambda}(\psi^{*})\\\
&\leq\tilde{E}_{\lambda}(\psi^{t})-\tilde{E}_{\lambda}(\psi_{\lambda}^{*})\\\
&\leq\frac{2}{\lambda(t+1)^{2}}\|\psi_{\lambda}^{*}\|^{2}\\\
&\leq\frac{\varepsilon}{2}\end{split}$

By Eqn. (14), we have

$\begin{split}\tilde{E}(\psi^{t})-\tilde{E}(\psi^{*})&=E(\psi^{t})+I_{H}(\psi^{t})-E(\psi^{*})-I_{H}(\psi^{*})\\\
&=(E(\psi^{t})-E_{\lambda}(\psi^{t}))+(E_{\lambda}(\psi^{t})+I_{H}(\psi^{t}))\\\
&~{}~{}~{}~{}-(E(\psi^{*})+I_{H}(\psi^{*}))\\\
&\leq|E(\psi^{t})-E_{\lambda}(\psi^{t})|+(\tilde{E}_{\lambda}(\psi^{t})-\tilde{E}(\psi^{*}))\\\
&\leq\lambda\log
n+\frac{\varepsilon}{2}=\frac{\varepsilon}{2}+\frac{\varepsilon}{2}\\\
&=\varepsilon\end{split}$

Next we show that $\|\psi_{\lambda}^{*}\|^{2}\leq n\|\bar{C}\|^{2}$ by proving
$|(\psi_{\lambda}^{*})_{j}|\leq\bar{C}~{}\forall j\in[n]$. According to Eq.
(13),

$\begin{split}\nu_{j}&=\sum_{i=1}^{m}\mu_{i}\frac{e^{((\psi_{\lambda}^{*})_{j}-c_{ij})/\lambda}}{\sum_{k=1}^{n}e^{((\psi_{\lambda}^{*})_{k}-c_{ik})/\lambda}}\\\
&=\sum_{i=1}^{m}\mu_{i}\frac{e^{(\psi_{\lambda}^{*})_{j}/\lambda}}{\sum_{k=1}^{n}e^{(\psi_{\lambda}^{*})_{k}/\lambda}e^{(c_{ij}-c_{ik})/\lambda}}\end{split}$
(17)

Assume $(\psi_{\lambda}^{*})_{\max}$ is the maximal element of
$\psi_{\lambda}^{*}$, we have
$\sum_{k=1}^{n}e^{(\psi_{\lambda}^{*})_{k}/\lambda}e^{(c_{ij}-c_{ik})/\lambda}\geq
e^{(\psi_{\lambda}^{*})_{\max}/\lambda}e^{-C_{\max}/\lambda}$, where
$C_{\max}$ is the maximal element of the matrix of $C$. Thus,

$\begin{split}\nu_{j}&\leq\sum_{i=1}^{m}\mu_{i}\frac{e^{(\psi_{\lambda}^{*})_{j}/\lambda}}{e^{(\psi_{\lambda}^{*})_{\max}/\lambda}e^{-C_{\max}/\lambda}}\\\
&=~{}~{}~{}~{}~{}~{}~{}~{}\frac{e^{(\psi_{\lambda}^{*})_{j}/\lambda}}{e^{(\psi_{\lambda}^{*})_{\max}/\lambda}e^{-C_{\max}/\lambda}}\end{split}$

Then,
$(\psi_{\lambda}^{*})_{\max}\leq(\psi_{\lambda}^{*})_{j}+C_{\max}-\lambda\log\nu_{j}$
and

$\begin{split}(\psi_{\lambda}^{*})_{\max}&\leq\frac{1}{n}\sum_{j=1}^{n}\\{(\psi_{\lambda}^{*})_{j}+C_{\max}-\lambda\log\nu_{j}\\}\\\
&\leq C_{\max}-\lambda\log\nu_{min}\end{split}$ (18)

According to the inequality of arithmetic and geometric means , we have
$\sum_{k=1}^{n}e^{(\psi_{\lambda}^{*})_{k}/\lambda}\geq
ne^{\frac{1}{n}(\sum_{k=1}^{n}(\psi_{\lambda}^{*})_{k}/\lambda)}=n$. Thus,
$\nu_{j}\leq\frac{e^{(\psi_{\lambda}^{*})_{j}/\lambda}}{ne^{-C_{\max}/\lambda}}$.

$\begin{split}(\psi_{\lambda}^{*})_{j}&\geq\lambda\log
n-C_{\max}+\lambda\log\nu_{j}\\\
&\geq\lambda\log\nu_{min}-C_{\max}\end{split}$ (19)

Combine Eqn. (18) and (19) together, we have $|(\psi_{\lambda}^{*})_{j}|\leq
C_{\max}-\lambda\log\nu_{min}=\bar{C}$. Hence, we obtain when
$t\geq\frac{\sqrt{8\|\bar{C}\|^{2}n\log n}}{\epsilon}$,
$\tilde{E}(\psi^{t})-\tilde{E}(\psi^{*})<\varepsilon$.

For each iteration in Eqn. (15), we need $O(n^{2})$ times of operations, thus
total the computational complexity of the proposed method is
$O(\frac{n^{2.5}\sqrt{\log n}}{\varepsilon})$. ∎


###### Abstract

Giant impacts can generate transient hydrogen-rich atmospheres, reducing
atmospheric carbon. The reduced carbon will form hazes that rain out onto the
surface and can become incorporated into the crust. Once heated, a large
fraction of the carbon would be converted into graphite. The result is that
local regions of the Hadean crust were plausibly saturated with graphite. We
explore the consequences of such a crust for a prebiotic surface hydrothermal
vent scenario. We model a surface vent fed by nitrogen-rich volcanic gas from
high-temperature magmas passing through graphite-saturated crust. We consider
this occurring at pressures of $1$–$1000\;{\rm bar}$ and temperatures of
$1500$–$1700\,{\rm{}^{\circ}C}$. The equilibrium with graphite purifies the
left-over gas, resulting in substantial quantities of nitriles ($0.1\%$ HCN
and $1\,{\rm ppm}$ HC_3N) and isonitriles ($0.01\%$ HNC) relevant for
prebiotic chemistry. We use these results to predict gas-phase concentrations
of methyl isonitrile of $\sim 1\,{\rm ppm}$. Methyl isocyanide can participate
in the non-enzymatic activation and ligation of the monomeric building blocks
of life, and surface, or shallow, hydrothermal environments provide its only
known equilibrium geochemical source.

###### keywords:

origin of life; volcanism on the early earth; hydrothermal vents; graphite

1 1 0 A Surface Hydrothermal Source of Nitriles and Isonitriles Paul B. Rimmer
1, Oliver Shorttle 2,3 Paul B. Rimmer, Oliver Shorttle

##  1 Introduction

Because the synthetic chemistry at life’s origins is a many-step process, a
requirement for a prebiotic environment is clean productive chemistry. If the
prebiotic environment is too diverse and complex (if the number of different
reacting species is too large), then the chemical parameter space inhabited by
a geochemical environment becomes large, and the desired products and
intermediates are lost in a morass of many thousands of other molecules. This
is what we mean by messy chemistry: chemistry of such diversity and
complexity, that desired chemical products and behaviors are harder to realize
Islam and Powner (2017); Sutherland (2017); Rimmer et al. (2018).

The requirement for clean chemistry is related to the arithmetic demon: if a
step-wise reaction does not provide products with high and selective yields,
and does not have a way of purifying and preserving desired products, then as
the fraction of useful product becomes the reactant for the next step
subsequent yields become exponentially diluted. Chemical reaction yields of
10% over a sequence of twenty steps will take a starting solution with
reactants at high concentration ($1\;{\rm M}$) to a final product with a
concentration of less than a molecule per cubic centimeter of solution. At
some point, this ceases to be chemistry and becomes homeopathy.

Prebiotic chemistry is more likely to be successful when it is clean and
productive. This places constraints on the environments in which successful
prebiotic chemistry can occur. A prebiotic environment that hosts clean
chemistry is one that facilitates selective, high-yield chemical reactions to
occur. A productive environment for prebiotic chemistry is one that
facilitates the synthesis of complex organic molecules. These two
environmental conditions are in tension with each other. The only conditions
that have been experimentally demonstrated to be productive are reducing
environments, and reducing environments tend to be messy. We will keep this
tension in mind as we consider one of the particularly promising chemical
starting points for prebiotic synthesis: nitriles.

Nitriles feature prominently in the majority of prebiotic systems chemistry
Ruiz-Mirazo et al. (2014); Benner et al. (2019). They carry the same redox
state of nitrogen found in biomolecules Sasselov et al. (2020), and have
remained a constant in the highly fruitful iterative discovery of
geochemically plausible UV-driven prebiotic synthesis of life’s monomeric
building blocks Sasselov et al. (2020); Green et al. (2021). Isonitrile
chemistry, specifically chemistry involving methyl isocyanide, has been
discovered to have astounding chemical properties of activating and ligating
monomeric building blocks Mariani et al. (2018), including nucleotides Liu et
al. (2020), phospholipids Bonfio et al. (2020) and amino acids Liu et al.
(2020); Wu et al. (2021).

Paradoxes often bear fruit in scientific exploration Benner (2014). In
addition to clean, productive chemistry, ideal conditions for the prebiotic
chemistry that forms nitriles admit several other paradoxes.

  * •

Access to ultraviolet light at wavelengths between $200$–$400\,{\rm nm}$ Patel
et al. (2015); and,
shielding of ultraviolet light at wavelengths between $200$–$400\,{\rm nm}$
Mariani et al. (2018).

  * •

Near-freezing temperatures Sanchez et al. (1966); Miyakawa et al. (2002);
Rimmer et al. (2018, 2021); and,
near-boiling temperatures Mansy and Szostak (2008).

  * •

Low pH Mariani et al. (2018); and,
neutral to high pH Mariani et al. (2018).

  * •

Water activity $\sim 1$, the conditions under which most nitrile-based
prebiotic chemistry takes place; and,
water activity $\ll 1$, required for phosphorylation of nucleosides, and
helpful for other condensation reactions Hulshof and Ponnamperuma (1976).

It is not possible for a single static environment to fulfill all of these
conditions. These can be satisfied in principle by a dynamic and heterogeneous
environment. We show that surface hydrothermal vents fed by gas from high-
temperature magma on early Earth can qualify.

Surface hydrothermal vents are exposed to ultraviolet light where they are in
contact with the atmosphere, and shielded from ultraviolet light at depth or
in crevices. Hydrothermal systems like the glaciovolcanic hydrothermal vents
on Iceland today have a wide range of temperatures, pH and chemistries, and
can provide good analogues for anoxic systems in the past Cousins et al.
(2013); Ilanko et al. (2019). For a single hydrothermal vent, temperatures can
be near freezing at their surfaces and above boiling at greater depth. Fluid
flow through small channels of rock with natural mineral buffers towards a
surface in contact with CO2, or phosphate-rich alkaline lake water Toner and
Catling (2020), can provide a steep gradient from high to low pH. Most of the
vent has water activity near unity, but the surface edges of the vent can dry
or freeze, lowering the water activity.

Hydrothermal vents, whether shallow surface vents or underwater vents, have
redox gradients generated by serpentinization and radiolysis Lin et al.
(2005); Shibuya et al. (2015). It is unlikely that these processes would have
generated sufficiently reducing conditions for the generation of nitriles or
isonitriles Lin et al. (2005); Klein et al. (2013); Rimmer and Shorttle
(2019). Primordial abiotic kerogen could provide sufficient reducing power,
but it is debatable whether significant amounts of kerogen was present in the
crust or upper mantle before life Gold and Soter (1980); Glasby (2006).

Even if reducing conditions are accessible, it is unlikely that they would
result in selective chemistry Rimmer and Shorttle (2019). Often the choice is
reducing conditions or chemical selection. Clean selective chemistry tends to
be oxidizing, or at least neutral. Reducing chemistry often results in tar Kim
et al. (2016).

In this paper, we show that clean productive chemistry rich in nitriles and
isonitriles can be found together in surface hydrothermal vents because of
graphitization. We present the surface vent scenario in Section 2. We discuss
the model used to predict the surface vent chemistry in Section 3 and show our
results in Section 4. Section 5 contains discussion and conclusions.

##  2 The Scenario

We present a scenario that we predict to result in clean and productive
prebiotic chemistry. It is a scenario that is both prebiotically plausible and
well supported by observation, experiment and models. A schematic of this
scenario is given in Figure 1. We will discuss some of the simplifying
assumptions of this scenario and other ways the same chemistry could emerge in
Section 5.

The Hadean Eon spans 500 million years of Earth’s history, from planet
formation to 4 Ga. After the moon-forming impact at $\sim 4.5\,{\rm Ga}$,
Earth likely had an atmosphere dominated by CO2 and N2, with some CO and H2O
Zahnle et al. (2010), around $1-5\%$ H2 Zahnle et al. (2018), and
comparatively trace amounts of sulfur-bearing compounds, SO2 and H2S Ranjan et
al. (2018), with low-to-mid, stable concentrations of sulfites in most natural
waters Ranjan et al. (2023).

At a time around $\sim 4.3\,{\rm Ga}$, Earth was likely hit by a roughly moon-
sized object Genda et al. (2017). The iron in this giant impactor would have
reacted with ocean water, producing large amounts of hydrogen at high
temperatures. Such a hydrogen-dominated atmosphere equilibrates with the
surface magma generated to give $\sim$1 bar partial pressure of H2 Itcovitz et
al. (2022). In the high temperatures of this post-impact atmosphere, hydrogen
would have reacted with carbon dioxide and nitrogen to produce methane and
ammonia Zahnle et al. (2020); Wogan et al. (2023). These giant impacts thereby
initiated transient and global highly-reducing conditions in the atmosphere
and on the crust of the Hadean Earth. Many of the reducing molecules,
hydrogen, methane and ammonia, are greenhouse gases, and the surface of Earth
at this time would have been hot, likely above the boiling point of liquid
water at 1 bar pressure Zahnle et al. (2020).

In these conditions, the hydrogen, methane, and either nitrogen or ammonia in
the atmosphere would have been photodissociated, with their products combining
to produce the nitriles HCN and HC3N Wogan et al. (2023). Many complex
organics would also have formed during this epoch, these would condense out of
the atmosphere, forming a tholin-like haze Trainer et al. (2004); Arney et al.
; Zahnle et al. (2020). This nitrogen-rich haze would have rained out onto the
hot surface as a thick tar Benner et al. (2020); Ritson et al. (2020); Zahnle
et al. (2020). It is likely some of this tar was incorporated into the crust,
either by tectonic, magmatic, or impact churning of the surface. Over the
period of a million years or more, the atmosphere would have returned to a
neutral chemical state through the conversion of methane and ammonia back to
carbon dioxide, nitrogen and hydrogen, and the escape of hydrogen into space
Zahnle et al. (2020); Wogan et al. (2023). The surface of Earth would then
cool to near freezing Kadoya et al. (2020).

Photochemically-produced tar mixed into the crust would have experienced
episodic heating to $>1500\,{\rm{}^{\circ}C}$ by the early high-temperature
(komatiitic) magmas known to have been an important constituent of the early
Earth magmatism Nisbet et al. (1993). We will show that this heating likely
broke apart the tar, transforming most of it into graphite, molecular hydrogen
and molecular nitrogen (though some of the hydrogen and nitrogen may have been
complexed with the graphite at this stage). Magmatic gas, interacting with the
graphitized crustal material, was likely transformed into HCN, HC3N and
isonitriles, along with sulfide, carbon monoxide, and little else. The result
is clean productive chemistry, which could degas through shallow and surface
hydrothermal vents on ancient volcanic islands.

Figure 1: A schematic representation of the scenario we propose here for clean
high-yield production of prebiotic feedstock. Events move around from top left
clockwise: first Earth has a neutral atmosphere; this is reduced following
giant impact at 4.3 Ga by oxidation of the impactor’s metal core to produce a
massive H2 atmosphere with significant methane and ammonia; This atmosphere
quickly cools (in $<1$ kyr) with photochemistry producing a tholin-rich haze
that deposits complex nitrogen-rich organics; these organics become
progressively buried and graphitized by interaction with magma, the atmosphere
clears as H2 is lost to space and becomes neutral again; finally, magmatic
gases interact with the graphite, and are scrubbed to produce high yields of
clean HCN, HC3N and isonitriles.

##  3 The Model

We will now model the latter part of this scenario: the interaction of gas
initially in equilibrium with $f_{O2}\approx{\rm QFM-1}$ magma (i.e., one log
unit below the quartz-magnetite-fayalite buffer in terms of its oxygen
fugacity) as it flows through graphite. We start with a fiducial model at
$1700\,{\rm{}^{\circ}C}$, 100 bar gas at QFM$-1$ plus a fixed nitrogen content
of $\ce{N}=5.7\%$. Graphite is added to the system. QFM$-1$ conditions are
incorporated by setting initial conditions to the values from Table 1.

We are not here invoking nitrogen-rich magma. The nitrogen is included in the
magmatic gas only for convenience; it makes no difference if it is included
initially, or included with the graphite, because the quantity is not varied.
The source of the nitrogen is expected to be the same as the source of the
graphite. This nitrogen concentration represents the predicted nitrogen-rich
nature of the post-impact organics Wogan et al. (2023), that were then
incorporated into the crust.

To predict the change in concentrations of chemical species as a function of
added graphite, we solve the following equation,

$\dfrac{d[\ce{X}]}{dt}=P_{\rm X}-L_{\rm X}\ce{[X]},$ (1)

where $[\ce{X}]\,{\rm(cm^{-3})}$ is the concentration of species X, $P_{\rm
X}\,{\rm(cm^{-3}\,s^{-1})}$ is the rate of production for that species, and
$L_{\rm X}\,{\rm(s^{-1})}$ is the rate of loss for that species. The terms
$P_{\rm X}$ and $L_{\rm X}$ are made up of concentrations of other species and
rate constants, $k$, which themselves have units depending on the reaction
order, and, when reversible, are set to reproduce chemical equilibrium by
assuring that, for the generic reaction:

$\ce{A}+\ce{B}+\ce{C}+...\rightleftarrows\ce{X}+\ce{Y}+\ce{Z}+...,$ (2)

with forward rate constant $k_{+}$ and reverse rate constant $k_{-}$, the
reverse rate constant is set such that:

$\dfrac{k_{+}}{k_{-}}=K_{\rm
eq}=e^{-\Delta_{f}G/RT}=\dfrac{\\{\ce{X}\\}\\{\ce{Y}\\}\\{\ce{Z}\\}...}{\\{\ce{A}\\}\\{\ce{B}\\}\\{\ce{C}\\}...},$
(3)

where $\\{\ce{A}\\}$, etc. are the activities of the different chemical
species, and $\Delta_{f}G/RT$ is determined using NASA coefficients, mostly
from Burcat & Ruscic Burcat and Ruscic (2005), and subsequent updates to the
database. More details can be found in prior presentations of the underlying
model Rimmer and Helling (2016); Rimmer et al. (2021).

For these calculations, we use an updated gas-phase chemical network based on
STAND-2020 Rimmer et al. (2021). This model includes H/C/N/O/S chemistry, a
very limited P network, and some reactions involving various heavier elements
such as Fe, Mg, Ti. It includes 6279 reactions involving 511 chemical species
The full list of species and network are available at
https://doi.org/10.7910/DVN/FKKYY3. The main addition is the two reactions:

$\ce{C(g)}\rightleftarrows\ce{C(s)},$ (4)

where C(g) is gas-phase carbon and C(s) is solid carbon as graphite and the
rate constants are set to reproduce equilibrium. The model is run for one day
model time.

The kinetics model and the FastChem model use different vapor pressures for
graphite. For the kinetics model we use Brewer et al. (1948); Joseph et al.
(2002):

$\log_{10}\,p_{\rm vap}=6.455-\dfrac{2.7709\times
10^{4}}{T-3.549}-\dfrac{10^{7}}{T^{2}},$ (5)

where $p\,({\rm bar})$ is the pressure and $T\,({\rm K})$ is the temperature.

We want to compare our kinetics results to equilibrium. We predict the
equilibrium of the magmatic gas using FastChem Stock et al. (2022), to which
we have added the thermochemical data for cyanoacetylene (HC3N; Knight et al.,
1985), given as the constant of mass action

$\ln{K(T)}=\dfrac{2.96\times 10^{5}}{T}-4.83\ln{T}-28.47+\big{(}1.91\times
10^{-3}\big{)}\,T-\big{(}1.10\times 10^{-7}\big{)}\,T^{2},$ (6)

where $T$ is the temperature in Kelvin. For FastChem the vapor pressure works
out to be approximately:

$\log_{10}\,p_{\rm vap,FC}=4.855-\dfrac{2.5709\times
10^{4}}{T}-\dfrac{10^{7}}{T^{2}}.$ (7)

We run our model for a wide range of conditions to determine the sensitivity
of our results to elemental composition, temperature and pressure. We run the
model from $1300\,{\rm{}^{\circ}C}$ to $1800\,{\rm{}^{\circ}C}$ and from
$1\,{\rm bar}$ to $1000\,{\rm bar}$, varying H, C, N and O elemental
abundances.

Our network does not include methyl isocyanide (CH3NC), the prebiotically-
relevant compound, and the gas-phase kinetics of this species would require
more investigation before it could be reliably included in this model. We show
HNC as an indication of the overall concentration of isonitriles, and we use
the BURCAT thermochemistry data to estimate the concentration of CH3NC.

We can use the Gibbs free energy $(\Delta_{r}G,\,{\rm kJ/mol})$ for the
following reaction:

$\ce{HCN}+\ce{CH4}\rightarrow\ce{CH3NC}+\ce{H2},$ (8)

to predict CH3NC concentration. This can be expressed as:

$[\ce{CH3NC}]=e^{-\Delta_{r}G/RT}\dfrac{[\ce{HCN}][\ce{CH4}]}{\ce{[H2]}}.$ (9)

where [HCN] and [CH4] are the concentrations of HCN and CH4, respectively;
$R=8.3145\;{\rm J/(mol\,K}$ is the gas constant, $T\,({\rm K})$ is the
temperature.

Table 1: Initial conditions for the kinetics model. We set temperature equal
to $1700\,{\rm{}^{\circ}C}$ and pressure is $100\,{\rm bar}$. Mixing ratios of
species are given below. Abundances have been rounded. Exact abundances used
can be found at https://doi.org/10.7910/DVN/FKKYY3.

Species |  C(g) |  CO |  H2 |  N2 |  O2 |  CH4 |  CO2 |  H2O
---|---|---|---|---|---|---|---|---
Mixing Ratio | * | $0.35$ | $0.1$ | $0.057$ | $7.6\times 10^{-10}$ | $1.4\times 10^{-6}$ | 0.23 | 0.26

* Graphite content is varied from zero to saturation (elemental abundance $\approx$ 80%)

##  4 Results

At $1700\,{\rm{}^{\circ}C}$ and 100 bar, the model achieves equilibrium after
$10^{3}$ – $10^{4}\,{\rm s}$ when the graphite concentration is $\lesssim
0.1$. The kinetics and equilibrium models then diverge. The equilibrium model
does not change with graphite added beyond a concentration of 0.8, which is
the saturation limit of graphite. The kinetics model continues to deviate from
equilibrium, but much more slowly. The results below saturation are given in
Figure 2 for HCN, HNC and HC3N, and Figure 3 for major species. A datafile of
the full results is available at https://doi.org/10.7910/DVN/FKKYY3. Results
are given as a function of carbon and oxygen fractions.

The main reason for the deviation between the two models is the choice of
graphite vapor pressure, where the kinetics model uses Eq. (5) and the
equilibrium model uses Eq. (6). If the kinetic model is run with the graphite
vapor pressure equal to Eq. (6), the results converge, as can be seen in
Figure 4.

We also show how the concentrations of HCN, HNC and HC3N at graphite
saturation depend on the magma temperature, see Figure 5. The results of the
sensitivity analysis are presented in the Appendix. This analysis was only run
for the chemical kinetics model.

Figure 2: Mixing ratios of HCN, HNC and HC3N (left y-axis) as a function of
added graphite (x-axis) for a QFM-1 buffered magmatic gas with $5.7\%$
elemental nitrogen content held at $1700\,{\rm{}^{\circ}C}$ and $100\,{\rm
bar}$. The right y-axis shows the deviation from QFM in log units. Solid lines
indicate the kinetics (this paper) results and dashed lines the equilibrium
(FastChem) results, with the range between the two shaded in. Figure 3: Mixing
ratios of major species (y-axis) for a QFM-buffered magmatic gas with $5.7\%$
elemental nitrogen content held at $1700\,{\rm{}^{\circ}C}$ and $100\,{\rm
bar}$. Solid lines indicate the kinetics (this paper) results and dashed lines
the equilibrium (FastChem) results, with the range between the two shaded in.
Figure 4: Mixing ratios of HCN, HNC and HC3N (left y-axis) as a function of
added graphite (x-axis) for a QFM-buffered magmatic gas with $5.7\%$ elemental
nitrogen content held at $1700\,{\rm{}^{\circ}C}$ and $100\,{\rm bar}$. Here
the kinetics model and the equilibrium model both use graphite vapor pressure
equal to Eq. (6). The right y-axis shows the deviation from QFM in log units.
Solid lines indicate the kinetics (this paper) results and dashed lines the
equilibrium (FastChem) results, with the range between the two shaded in.
Figure 5: Mixing ratios of HCN, HNC and HC3N (y-axis) at graphite saturation
as a function of temperature (∘C, x-axis) for a magmatic gas stating at
QFM$-1$ with $5.7\%$ elemental nitrogen content held at $100\,{\rm bar}$. Only
the kinetics results are shown. Indicative eruptive temperatures of
terrestrial magmas are given above the plot: ‘MORB’ = mid-ocean ridge basalts
(Coogan et al., 2014), ‘LIP’ = Large igneous province (Coogan et al., 2014),
and komatiite temperatures are given for Archean-age examples (Nisbet et al.,
1993).

The most important result is that graphite formation cleans up the chemistry
considerably. While Rimmer & Shorttle Rimmer and Shorttle (2019) predict that
the majority carbon-containing species is diacetylene, our study finds that
the majority carbon-containing species is graphite, and diacetylene has
decreased from $42\%$ to $<0.1\%$. This is significant because diacetylene at
lower temperatures will polymerize and will plausibly generate a mass of large
inert hydrocarbons, effective “tarrification” Kim et al. (2016). Graphite
formation resolves much of this problem. Since graphite is removed from the
gas phase, the remaining therochemically stable gas is enriched, and so much
higher concentrations of cyanide (HCN) and cyanoacetylene (HC_3N) are
achieved, to a maximum of $\sim 1\%$ and $\sim 1\,{\rm ppm}$ of the gas-phase
respectively.

An unexpected prediction is the formation of higher-than-expected
concentrations of hydrogen isocyanide (HNC). At graphite saturation,
$[\ce{HNC}]\sim 0.01\%$.

We use the data from BURCAT to calculate the Gibbs free energy of Reaction
(8), and show this free energy in Figure 6. We apply the value of the Gibbs
free energy at $1700\,{\rm{}^{\circ}C}$ of $114\,{\rm kJ/mol}$, and
$[\ce{H2}]=0.2$, $[\ce{CH4}]=3-7\times 10^{-3}$ and
$[\ce{HCN}]=10^{-3}-10^{-2}$ to Equation (9) to predict the equilibrium
concentration: $[\ce{CH3NC}]$ is between $7\times 10^{-7}$ and $2\times
10^{-5}$.

Figure 6: The Gibbs free energy of the reaction:
$\ce{HCN}+\ce{CH4}\rightarrow\ce{CH3NC}+\ce{H2}$ (kJ/mol), as a function of
temperature (K). We can see that for all temperatures, this formation reaction
for CH3NC is endergonic, but not strongly so.

At lower pressure, if the magma is graphite saturated, then HCN, HNC and HC3N
remain effectively unchanged. In systems which are not graphite saturated, HCN
and, in the right range, HNC become even more favored, with HNC achieving
mixing ratios of $\sim 0.1\%$ of the gas phase. Cyanoacetylene drops
significantly, by roughly two orders of magnitude per order of magnitude of
decreasing pressure, to $\lesssim 100\,{\rm ppt}$ levels at 1 bar. In any
event, experiments to test the predictions of this model at ambient pressure,
when the system is graphite saturated, may be tractable in the near future.

At lower temperature, graphite reigns. Rimmer & Shorttle Rimmer and Shorttle
(2019) predict that, at $1200\,{\rm{}^{\circ}C}$, significant amounts of HCN
and HC3N could still be formed at $>10\,{\rm bar}$ pressures. With our present
model, the carbon is largely removed into graphite, and species like methane,
molecular hydrogen and carbon monoxide dominate. For model results over a full
range of oxygen, carbon and hydrogen fractions, see Appendix A.

##  5 Discussion and Conclusions

In this paper, we present a new calculation of magmatic gas-phase chemistry
for nitrogen-rich, otherwise of standard oxidation state (QFM$-1$), magmas at
$1700{\rm{}^{\circ}C}$ temperatures and 100 bar pressure, with graphite added
to the point of saturation. This is the specific sequence we have modeled, but
it may well turn out that the QFM$-1$ magma is nitrogen-poor and instead that
at temperatures of $\sim 1700{\rm{}^{\circ}C}$ nitrogen is made available to
the gas from the its complex with graphite. This model differs from Rimmer &
Shorttle 2019 Rimmer and Shorttle (2019) in that we now include phase
equilibrium with graphite. We find that graphite helps to clean up the
chemistry, and increases the concentrations of prebiotically relevant
compounds, nitriles and isonitriles, if the magma temperature is $\sim
1700\,{\rm{}^{\circ}C}$, consistent with early experimental studies of
nitrogen and graphite at high temperatures Palmer (1967). If the magma
temperature is lower, the concentration of nitriles and isonitriles drops
significantly, see Figure 5.

Including graphite is not predicted to result in a decrease of nitriles so
long as the temperature is sufficiently high ($\gtrsim 1500{\rm{}^{\circ}C}$).
The results at lower temperatures indeed show a significant suppression of
nitriles, in line with previous model expectations for hydrogen solubility,
such as those of Wogan et al. Wogan et al. (2020). It is important to note
that these environments are locally at least an order of magnitude more
reducing than the minimum $f_{\ce{O2}}$ considered by Wogan et al. Even under
these conditions, the roughly estimated kg gas/kg magma for CH4 at 100 bar and
1700 K would be $\lesssim 10^{-4}$, in reasonable agreement with Wogan et al.
If the magma was still present at the point of gas-graphite interaction,
solubility experiments would need to be performed under these unusual
conditions in order to determine the true fate of gas-phase nitriles and other
hydrogen- and nitrogen-bearing species.

The required conditions for clean chemistry with nitriles and isonitriles is
graphite-saturated magma at $\gtrsim 1500{\rm{}^{\circ}C}$, else the carbon
and nitrogen that would form nitriles is locked up in N2 and graphite. See the
appendix for plots showing the dependence of the nitrile concentration on
temperature and concentration of carbon (Figures 15 \- 17). Cyanoacetylene is
favored over a wide range of parameter space at higher pressure ($\gtrsim
100\,{\rm bar}$; see the Appendix, Figures 11 \- 14), but can we predict that
it can be produced in abundance at pressures as low as $1\,{\rm bar}$ so long
as the system is saturated with graphite. See also Table 2 for a summary of
these comparisons.

The required high temperatures effectively limit the source of nitriles and
isonitriles to more niche Hadean environments, particularly those that run at
higher temperatures, such as Komatiite magmas Nisbet et al. (1993); Herzberg
et al. (2010), for which temperatures can surpass $1600\,^{\circ}C$ Takahashi
and Scarfe (1985). These high-temperature magmas are thought to be more
prevalent (though not ubiquitous) on the Hadean and Archaean Earth Campbell et
al. (1989); Nebel et al. (2014). This development of the theory admits to
several potential scenarios beyond the one we presented in Section 2.

It may simply be that certain regions of the early upper mantle were highly
reduced, possessing abiotically-generated kerogen-like material, in a manner
hypothesized by Thomas Gold, among others Gold and Soter (1980). The natural
high temperatures of early volcanic systems could heat and reprocess this
material.

Cosmic dust is kerogen-like, and would have been much more ubiquitous on early
Earth, with as much as 50% of the dust being cosmic during the Hadean Walton
et al. (2023). The chemistry predicted here could arise in environments where
this dust is concentrated and then heated volcanically, or by a subsequent
large impact.

Radiolysis and serpentinization can generate large redox gradients Lin et al.
(2005); Rempfert et al. (2023), with certain regions of the crust becoming
much more reduced while others become more oxidized. If the reduced regions
intersected a magmatic flow at $>10\,{\rm bar}$, the heat would convert much
of the chemistry to nitriles and isonitriles. Radiolysis especially would have
been much more intense during the Hadean Pastorek et al. (2020).

Though we favor the scenario presented in Section 2 for Early Earth, based on
current geological and experimental evidence, we are encouraged that similar
prebiotic chemistry can emerge in a wide range of conditions expected on early
Earth, Mars, and on exoplanets. It is also worth noting that all of these
scenarios are compatible with shallow hydrothermal systems, which admit most
of the physical and chemical advantages of underwater hydrothermal events
Barge and Price (2022). It is also worthy of note that some of these systems,
such as the surface hydrothermal fields near Erebus volcano, Antarctica,
exhibit these vast temperature shifts, with high-temperature magmatic systems
intersecting with vents that release gas into ice Ilanko et al. (2019), a
system that would favor concentrating some of these prebiotic feed-stocks into
eutectic phases. In addition, these systems show fascinating possibly-abiotic
redox behavior, and can incorporate nitrogen from the atmosphere into their
magmatic systems Ilanko et al. (2019).

The provision of chemical feedstocks from high-temperature magmas may pave the
way for order amidst geochemical chaos, with clean chemical equilibrium
mixtures at high-concentration, segregated by pool or stream on the basis of
the magmatic source: magmas with higher carbon content provide isonitriles,
such as methyl isocyanide. All of this is mediated and regulated by the
formation of graphite, removing from the gas excess carbon and the
combinatorial mess that comes with it, and high temperatures, that favor
relatively simple gas-phase mixtures of the starting material required for
productive prebiotic chemistry. Natural prebiotic environments need not
produce messy chemistry; the environment can constrain the chemistry to be
clean and productive. Even if the environments favouring more chemically
ordered and promising prebiotic chemistry turn out to be a very small fraction
of the total environment, the selection pressure for productive synthesis
could outweigh their relative rarity.

Whether this chemical solution provides a “buffet lunch” for prebiotic
chemistry, or an unappetizing and unusable mess depends on the kinetic
stability and solubility of the gas-phase mixture once it is quenched and
enters into the surface waters. This question can only be resolved with future
experiments.

This research was funded by the IPLU grant on “Quantifying Aliveness”
(Cambridge Planetary Science and Life in the Universe Research Grants Scheme
2021-22).

All data required to reproduce the results of this paper can be found at
https://doi.org/10.7910/DVN/FKKYY3.

###### Acknowledgements.

We thank Craig Walton and David Catling for helpful conversations about this
work. We thank the anonymous referees and Shang-Min Tsai for constructive
comments on the draft that helped to improve the work. The author declares no
conflict of interest. no

##  Appendix A Sensitivity Analysis

Here we show the results of varying the gas-phase chemistry of the magma, by
adjusting the elemental abundances of H, C and O with a fixed N at 5%. These
results are shown for a pressure of 100 bar and a temperature of
$1600\rm{\,{}^{\circ}C}$: Figures 7 to 10, as well as for a pressure of 10 bar
and a temperature of $1600\rm{\,{}^{\circ}C}$: 11 to 14, and for a pressure of
100 bar and a temperature of $1300\rm{\,{}^{\circ}C}$: 15 to 17. For the last
case, HC3N is not shown because its abundance is below 1 ppm over the entire
parameter space. Sensitivities based on these results are summarized in Table
2. Full results, including all molecules predicted by the model and network,
are given at https://doi.org/10.7910/DVN/FKKYY3.

We find that in all cases that nitrile concentrations change linearly when
varying the abundance of N. The results when varying N are therefore not
plotted.

Table 2: Sensitivity of species concentrations on temperature and pressure, expressed as ratios. $[\ce{X}_{a}/\ce{X}_{b}]$ is the ratio of $[\ce{X}]$ under the conditions specified by $a$ over $[\ce{X}]$ under the conditions specified by $b$. $a$ and $b$ are given in the top row. Unless otherwise stated, temperature is $1700\,{\rm{}^{\circ}C}$ and pressure is $100\,{\rm bar}$.  Species | Temperature ($1300\,{\rm{}^{\circ}C}_{a}$/$1700\,{\rm{}^{\circ}C}_{b}$) | Pressure $\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;({\rm 10\,bar_{a}/100\,bar_{b})}$
---|---|---
$[\ce{HCN}_{a}/\ce{HCN}_{b}]$ | $10^{-3}$ | $10^{-1}$
$[\ce{HC3N}_{a}/\ce{HC3N}_{b}]$ | $10^{-5}$ | $10^{-3}$
$[\ce{HNC}_{a}/\ce{HNC}_{b}]$ | $10^{-3}$ | $10^{-1}$
Figure 7: Fraction of graphite in the mixture as a function of carbon (x-axis)
and oxygen (y-axis) fractions for a magmatic gas held at
$1600\,{\rm{}^{\circ}C}$ and $100\,{\rm bar}$. The nitrogen fraction is set to
$5\%$. The hydrogen fraction makes up the difference, if any. Figure 8:
Hydrogen cyanide concentrations, as $\log_{10}([\ce{HCN}])$, as a function of
carbon (x-axis) and oxygen (y-axis) fractions for a magmatic gas held at
$1600\,{\rm{}^{\circ}C}$ and $100\,{\rm bar}$. The nitrogen fraction is set to
$5\%$. The hydrogen fraction makes up the difference, if any.  Figure 9:
Cyanoacetylene concentrations, as $\log_{10}([\ce{HC_{3}N}])$, as a function
of carbon (x-axis) and oxygen (y-axis) fractions for a magmatic gas held at
$1600\,{\rm{}^{\circ}C}$ and $100\,{\rm bar}$. The nitrogen fraction is set to
$5\%$. The hydrogen fraction makes up the difference, if any.  Figure 10:
Hydrogen isocyanide concentrations, as $\log_{10}([\ce{HNC}])$, as a function
of carbon (x-axis) and oxygen (y-axis) fractions for a magmatic gas held at
$1600\,{\rm{}^{\circ}C}$ and $100\,{\rm bar}$. The nitrogen fraction is set to
$5\%$. The hydrogen fraction makes up the difference, if any.  Figure 11:
Fraction of graphite in the mixture as a function of carbon (x-axis) and
oxygen (y-axis) fractions for a magmatic gas held at $1600\,{\rm{}^{\circ}C}$
and $10\,{\rm bar}$. The nitrogen fraction is set to $5\%$. The hydrogen
fraction makes up the difference, if any. Figure 12: Hydrogen cyanide
concentrations, as $\log_{10}([\ce{HCN}])$, as a function of carbon (x-axis)
and oxygen (y-axis) fractions for a magmatic gas held at
$1600\,{\rm{}^{\circ}C}$ and $10\,{\rm bar}$. The nitrogen fraction is set to
$5\%$. The hydrogen fraction makes up the difference, if any.  Figure 13:
Cyanoacetylene concentrations, as $\log_{10}([\ce{HC_{3}N}])$, as a function
of carbon (x-axis) and oxygen (y-axis) fractions for a magmatic gas held at
$1600\,{\rm{}^{\circ}C}$ and $10\,{\rm bar}$. The nitrogen fraction is set to
$5\%$. The hydrogen fraction makes up the difference, if any.  Figure 14:
Hydrogen isocyanide concentrations, as $\log_{10}([\ce{HNC}])$, as a function
of carbon (x-axis) and oxygen (y-axis) fractions for a magmatic gas held at
$1600\,{\rm{}^{\circ}C}$ and $10\,{\rm bar}$. The nitrogen fraction is set to
$5\%$. The hydrogen fraction makes up the difference, if any.  Figure 15:
Fraction of graphite in the mixture as a function of carbon (x-axis) and
oxygen (y-axis) fractions for a magmatic gas held at $1300\,{\rm{}^{\circ}C}$
and $100\,{\rm bar}$. The nitrogen fraction is set to $5\%$. The hydrogen
fraction makes up the difference, if any. Figure 16: Hydrogen cyanide
concentrations, as $\log_{10}([\ce{HCN}])$, as a function of carbon (x-axis)
and oxygen (y-axis) fractions for a magmatic gas held at
$1300\,{\rm{}^{\circ}C}$ and $100\,{\rm bar}$. The nitrogen fraction is set to
$5\%$. The hydrogen fraction makes up the difference, if any.  Figure 17:
Hydrogen isocyanide concentrations, as $\log_{10}([\ce{HNC}])$, as a function
of carbon (x-axis) and oxygen (y-axis) fractions for a magmatic gas held at
$1300\,{\rm{}^{\circ}C}$ and $100\,{\rm bar}$. The nitrogen fraction is set to
$5\%$. The hydrogen fraction makes up the difference, if any.

References

## References

  * Islam and Powner (2017) Islam, S.; Powner, M.W.  Prebiotic Systems Chemistry: Complexity Overcoming Clutter.  Chem 2017, 2, 470–501.  https://doi.org/10.1016/j.chempr.2017.03.001.
  * Sutherland (2017) Sutherland, J.D.  Opinion: Studies on the origin of life—the end of the beginning.  Nature Reviews Chemistry 2017, 1, 0012.
  * Rimmer et al. (2018) Rimmer, P.B.; Xu, J.; Thompson, S.J.; Gillen, E.; Sutherland, J.D.; Queloz, D.  The origin of RNA precursors on exoplanets.  Science Advances 2018, 4, eaar3302, [arXiv:astro-ph.EP/1808.02718].  https://doi.org/10.1126/sciadv.aar3302.
  * Ruiz-Mirazo et al. (2014) Ruiz-Mirazo, K.; Briones, C.; de la Escosura, A.  Prebiotic systems chemistry: new perspectives for the origins of life.  Chemical reviews 2014, 114, 285–366.
  * Benner et al. (2019) Benner, S.A.; Kim, H.J.; Biondi, E.  Prebiotic Chemistry that Could Not Not Have Happened.  Life 2019, 9, 84.  https://doi.org/10.3390/life9040084.
  * Sasselov et al. (2020) Sasselov, D.D.; Grotzinger, J.P.; Sutherland, J.D.  The origin of life as a planetary phenomenon.  Science Advances 2020, 6, eaax3419.  https://doi.org/10.1126/sciadv.aax3419.
  * Green et al. (2021) Green, N.J.; Xu, J.; Sutherland, J.D.  Illuminating life’s origins: UV photochemistry in abiotic synthesis of biomolecules.  Journal of the American Chemical Society 2021, 143, 7219–7236.
  * Mariani et al. (2018) Mariani, A.; Russell, D.A.; Javelle, T.; Sutherland, J.D.  A light-releasable potentially prebiotic nucleotide activating agent.  Journal of the American Chemical Society 2018, 140, 8657–8661.
  * Liu et al. (2020) Liu, Z.; Wu, L.F.; Xu, J.; Bonfio, C.; Russell, D.A.; Sutherland, J.D.  Harnessing chemical energy for the activation and joining of prebiotic building blocks.  Nature Chemistry 2020, 12, 1–6.  https://doi.org/10.1038/s41557-020-00564-3.
  * Bonfio et al. (2020) Bonfio, C.; Russell, D.A.; Green, N.J.; Mariani, A.; Sutherland, J.D.  Activation chemistry drives the emergence of functionalised protocells.  Chemical Science 2020, 11, 10688–10697.
  * Wu et al. (2021) Wu, L.F.; Liu, Z.; Sutherland, J.D.  pH-Dependent peptide bond formation by the selective coupling of $\alpha$-amino acids in water.  Chemical Communications 2021, 57, 73–76.
  * Benner (2014) Benner, S.A.  Paradoxes in the Origin of Life.  Origins of Life and Evolution of the Biosphere 2014, 44, 339–343.  https://doi.org/10.1007/s11084-014-9379-0.
  * Patel et al. (2015) Patel, B.H.; Percivalle, C.; Ritson, D.J.; Duffy, C.D.; Sutherland, J.D.  Common origins of RNA, protein and lipid precursors in a cyanosulfidic protometabolism.  Nature Chemistry 2015, 7, 301–307.  https://doi.org/10.1038/nchem.2202.
  * Sanchez et al. (1966) Sanchez, R.; Ferris, J.; Orgel, L.E.  Conditions for Purine Synthesis: Did Prebiotic Synthesis Occur at Low Temperatures?  Science 1966, 153, 72–73.  https://doi.org/10.1126/science.153.3731.72.
  * Miyakawa et al. (2002) Miyakawa, S.; James Cleaves, H.; Miller, S.L.  The Cold Origin of Life: A. Implications Based On The Hydrolytic Stabilities Of Hydrogen Cyanide And Formamide.  Origins of Life and Evolution of the Biosphere 2002, 32, 195–208.  https://doi.org/10.1023/A:1016514305984.
  * Rimmer et al. (2021) Rimmer, P.B.; Thompson, S.J.; Xu, J.; Russell, D.A.; Green, N.J.; Ritson, D.J.; Sutherland, J.D.; Queloz, D.P.  Timescales for Prebiotic Photochemistry Under Realistic Surface Ultraviolet Conditions.  Astrobiology 2021, 21, 1099–1120.  https://doi.org/10.1089/ast.2020.2335.
  * Mansy and Szostak (2008) Mansy, S.S.; Szostak, J.W.  Thermostability of model protocell membranes.  Proceedings of the National Academy of Science 2008, 105, 13351–13355.  https://doi.org/10.1073/pnas.0805086105.
  * Mariani et al. (2018) Mariani, A.; Bonfio, C.; Johnson, C.M.; Sutherland, J.D.  pH-Driven RNA strand separation under prebiotically plausible conditions.  Biochemistry 2018, 57, 6382–6386.
  * Hulshof and Ponnamperuma (1976) Hulshof, J.; Ponnamperuma, C.  Prebiotic condensation reactions in an aqueous medium: A review of condensing agents.  Origins of Life 1976, 7, 197–224.  https://doi.org/10.1007/BF00926938.
  * Cousins et al. (2013) Cousins, C.R.; Crawford, I.A.; Carrivick, J.L.; Gunn, M.; Harris, J.; Kee, T.P.; Karlsson, M.; Carmody, L.; Cockell, C.; Herschy, B.; et al.  Glaciovolcanic hydrothermal environments in Iceland and implications for their detection on Mars.  Journal of Volcanology and Geothermal Research 2013, 256, 61–77.  https://doi.org/10.1016/j.jvolgeores.2013.02.009.
  * Ilanko et al. (2019) Ilanko, T.; Fischer, T.P.; Kyle, P.; Curtis, A.; Lee, H.; Sano, Y.  Modification of fumarolic gases by the ice-covered edifice of Erebus volcano, Antarctica.  Journal of Volcanology and Geothermal Research 2019, 381, 119–139.  https://doi.org/10.1016/j.jvolgeores.2019.05.01710.31223/osf.io/jsc25.
  * Toner and Catling (2020) Toner, J.D.; Catling, D.C.  A carbonate-rich lake solution to the phosphate problem of the origin of life.  Proceedings of the National Academy of Science 2020, 117, 883–888.  https://doi.org/10.1073/pnas.1916109117.
  * Lin et al. (2005) Lin, L.H.; Hall, J.; Lippmann-Pipke, J.; Ward, J.A.; Sherwood Lollar, B.; Deflaun, M.; Rothmel, R.; Moser, D.; Gihring, T.M.; Mislowack, B.; et al.  Radiolytic H2 in continental crust: Nuclear power for deep subsurface microbial communities.  Geochemistry, Geophysics, Geosystems 2005, 6, Q07003.  https://doi.org/10.1029/2004GC000907.
  * Shibuya et al. (2015) Shibuya, T.; Yoshizaki, M.; Sato, M.; Shimizu, K.; Nakamura, K.; Omori, S.; Suzuki, K.; Takai, K.; Tsunakawa, H.; Maruyama, S.  Hydrogen-rich hydrothermal environments in the Hadean ocean inferred from serpentinization of komatiites at 300 C and 500 bar.  Progress in Earth and Planetary Science 2015, 2, 1–11.
  * Klein et al. (2013) Klein, F.; Bach, W.; McCollom, T.M.  Compositional controls on hydrogen generation during serpentinization of ultramafic rocks.  Lithos 2013, 178, 55–69.  https://doi.org/10.1016/j.lithos.2013.03.008.
  * Rimmer and Shorttle (2019) Rimmer, P.; Shorttle, O.  Origin of Life’s Building Blocks in Carbon- and Nitrogen-Rich Surface Hydrothermal Vents.  Life 2019, 9, 12, [arXiv:astro-ph.EP/1901.08542].  https://doi.org/10.3390/life9010012.
  * Gold and Soter (1980) Gold, T.; Soter, S.  The Deep-Earth-Gas Hypothesis.  Scientific American 1980, 242, 154–161.  https://doi.org/10.1038/scientificamerican0680-154.
  * Glasby (2006) Glasby, G.P.  Abiogenic origin of hydrocarbons: An historical overview.  Resource Geology 2006, 56, 83–96.
  * Kim et al. (2016) Kim, H.J.; Furukawa, Y.; Kakegawa, T.; Bita, A.; Scorei, R.; Benner, S.A.  Evaporite Borate-Containing Mineral Ensembles Make Phosphate Available and Regiospecifically Phosphorylate Ribonucleosides: Borate as a Multifaceted Problem Solver in Prebiotic Chemistry.  Angewandte Chemie 2016, 128, 16048–16052.  https://doi.org/10.1002/ange.201608001.
  * Zahnle et al. (2010) Zahnle, K.; Schaefer, L.; Fegley, B.  Earth’s earliest atmospheres.  Cold Spring Harbor perspectives in biology 2010, 2, a004895.
  * Zahnle et al. (2018) Zahnle, K.J.; Gacesa, M.; Catling, D.C.  Strange messenger: A new history of hydrogen on Earth as told by xenon.  In Proceedings of the AGU Fall Meeting Abstracts, 2018, Vol. 2018, pp. P44B–01.
  * Ranjan et al. (2018) Ranjan, S.; Todd, Z.R.; Sutherland, J.D.; Sasselov, D.D.  Sulfidic Anion Concentrations on Early Earth for Surficial Origins-of-Life Chemistry.  Astrobiology 2018, 18, 1023–1040, [arXiv:astro-ph.EP/1801.07725].  https://doi.org/10.1089/ast.2017.1770.
  * Ranjan et al. (2023) Ranjan, S.; Abdelazim, K.; Lozano, G.G.; Mandal, S.; Zhou, C.Y.; Kufner, C.L.; Todd, Z.R.; Sahai, N.; Sasselov, D.D.  Geochemical and photochemical constraints on S [IV] concentrations in natural waters on prebiotic Earth.  AGU Advances 2023, 4, e2023AV000926.
  * Genda et al. (2017) Genda, H.; Brasser, R.; Mojzsis, S.J.  The terrestrial late veneer from core disruption of a lunar-sized impactor.  Earth and Planetary Science Letters 2017, 480, 25–32, [arXiv:astro-ph.EP/1709.07554].  https://doi.org/10.1016/j.epsl.2017.09.041.
  * Itcovitz et al. (2022) Itcovitz, J.P.; Rae, A.S.P.; Citron, R.I.; Stewart, S.T.; Sinclair, C.A.; Rimmer, P.B.; Shorttle, O.  Reduced Atmospheres of Post-impact Worlds: The Early Earth.  Planetary Science Journal 2022, 3, 115, [arXiv:astro-ph.EP/2204.09946].  https://doi.org/10.3847/PSJ/ac67a9.
  * Zahnle et al. (2020) Zahnle, K.J.; Lupu, R.; Catling, D.C.; Wogan, N.  Creation and Evolution of Impact-generated Reduced Atmospheres of Early Earth.  Planetary Science Journal 2020, 1, 11, [arXiv:astro-ph.EP/2001.00095].  https://doi.org/10.3847/PSJ/ab7e2c.
  * Wogan et al. (2023) Wogan, N.F.; Catling, D.C.; Zahnle, K.J.; Lupu, R.  Origin-of-life Molecules in the Atmosphere after Big Impacts on the Early Earth.  Planetary Science Journal 2023, 4, 169, [arXiv:astro-ph.EP/2307.09761].  https://doi.org/10.3847/PSJ/aced83.
  * Trainer et al. (2004) Trainer, M.G.; Pavlov, A.A.; Curtis, D.B.; McKay, C.P.; Worsnop, D.R.; Delia, A.E.; Toohey, D.W.; Toon, O.B.; Tolbert, M.A.  Haze Aerosols in the Atmosphere of Early Earth: Manna from Heaven.  Astrobiology 2004, 4, 409–419.  https://doi.org/10.1089/ast.2004.4.409.
  * (39) Arney, G.N.; Meadows, V.S.; Domagal-Goldman, S.D.; Deming, D.; Robinson, T.D.; Tovar, G.; Wolf, E.T.; Schwieterman, E.  Pale Orange Dots: The Impact of Organic Haze on the Habitability and Detectability of Earthlike Exoplanets.  Astrophysical Journal, 836, 49, [arXiv:astro-ph.EP/1702.02994].  https://doi.org/10.3847/1538-4357/836/1/49.
  * Benner et al. (2020) Benner, S.A.; Bell, E.A.; Biondi, E.; Brasser, R.; Carell, T.; Kim, H.J.; Mojzsis, S.J.; Omran, A.; Pasek, M.A.; Trail, D.  When did life likely emerge on Earth in an RNA-first process?  ChemSystemsChem 2020, 2, e1900035.
  * Ritson et al. (2020) Ritson, D.J.; Mojzsis, S.J.; Sutherland, J.D.  Supply of phosphate to early Earth by photogeochemistry after meteoritic weathering.  Nature Geoscience 2020, 13, 344–348.  https://doi.org/10.1038/s41561-020-0556-7.
  * Kadoya et al. (2020) Kadoya, S.; Krissansen-Totton, J.; Catling, D.C.  Probable Cold and Alkaline Surface Environment of the Hadean Earth Caused by Impact Ejecta Weathering.  Geochemistry, Geophysics, Geosystems 2020, 21, e2019GC008734.  https://doi.org/10.1029/2019GC008734.
  * Nisbet et al. (1993) Nisbet, E.G.; Cheadle, M.J.; Arndt, N.T.; Bickle, M.J.  Constraining the potential temperature of the Archaean mantle: A review of the evidence from komatiites.  Lithos 1993, 30, 291–307.  https://doi.org/10.1016/0024-4937(93)90042-B.
  * Burcat and Ruscic (2005) Burcat, A.; Ruscic, B.  Third millenium ideal gas and condensed phase thermochemical database for combustion (with update from active thermochemical tables).  Technical report, Argonne National Lab.(ANL), Argonne, IL (United States), 2005.
  * Rimmer and Helling (2016) Rimmer, P.B.; Helling, C.  A Chemical Kinetics Network for Lightning and Life in Planetary Atmospheres.  Astrophysical Journal Supplemental Series 2016, 224, 9, [arXiv:astro-ph.EP/1510.07052].  https://doi.org/10.3847/0067-0049/224/1/9.
  * Rimmer et al. (2021) Rimmer, P.B.; Jordan, S.; Constantinou, T.; Woitke, P.; Shorttle, O.; Hobbs, R.; Paschodimas, A.  Hydroxide Salts in the Clouds of Venus: Their Effect on the Sulfur Cycle and Cloud Droplet pH.  Planetary Science Journal 2021, 2, 133, [arXiv:astro-ph.EP/2101.08582].  https://doi.org/10.3847/PSJ/ac0156.
  * Brewer et al. (1948) Brewer, L.; Gilles, P.W.; Jenkins, F.A.  The Vapor Pressure and Heat of Sublimation of Graphite.  Journal of Chemical Physics 1948, 16, 797–807.  https://doi.org/10.1063/1.1746999.
  * Joseph et al. (2002) Joseph, M.; Sivakumar, N.; Manoravi, P.  High temperature vapour pressure studies on graphite using laser pulse heating.  Carbon 2002, 40, 2031–2034.  https://doi.org/10.1016/S0008-6223(02)00158-6.
  * Stock et al. (2022) Stock, J.W.; Kitzmann, D.; Patzer, A.B.C.  FASTCHEM 2 : an improved computer program to determine the gas-phase chemical equilibrium composition for arbitrary element distributions.  Monthly Notices of the Royal Astronomical Society 2022, 517, 4070–4080, [arXiv:astro-ph.EP/2206.08247].  https://doi.org/10.1093/mnras/stac2623.
  * Knight et al. (1985) Knight, J.; Freeman, C.; McEwan, M.; Adams, N.; Smith, D.  Selected-ion flow tube studies of HC3N.  International journal of mass spectrometry and ion processes 1985, 67, 317–330.
  * Coogan et al. (2014) Coogan, L.; Saunders, A.; Wilson, R.  Aluminum-in-olivine thermometry of primitive basalts: Evidence of an anomalously hot mantle source for large igneous provinces.  Chemical Geology 2014, 368, 1–10.
  * Palmer (1967) Palmer, D.J.  Reaction of Hydrogen with Graphite.  Nature 1967, 215, 388–389.  https://doi.org/10.1038/215388a0.
  * Wogan et al. (2020) Wogan, N.; Krissansen-Totton, J.; Catling, D.C.  Abundant Atmospheric Methane from Volcanism on Terrestrial Planets Is Unlikely and Strengthens the Case for Methane as a Biosignature.  Planetary Science Journal 2020, 1, 58, [arXiv:astro-ph.EP/2009.07761].  https://doi.org/10.3847/PSJ/abb99e.
  * Herzberg et al. (2010) Herzberg, C.; Condie, K.; Korenaga, J.  Thermal history of the Earth and its petrological expression.  Earth and Planetary Science Letters 2010, 292, 79–88.
  * Takahashi and Scarfe (1985) Takahashi, E.; Scarfe, C.M.  Melting of peridotite to 14 GPa and the genesis of komatiite.  Nature 1985, 315, 566–568.  https://doi.org/10.1038/315566a0.
  * Campbell et al. (1989) Campbell, I.H.; Griffiths, R.W.; Hill, R.I.  Melting in an Archaean mantle plume: heads it’s basalts, tails it’s komatiites.  Nature 1989, 339, 697–699.  https://doi.org/10.1038/339697a0.
  * Nebel et al. (2014) Nebel, O.; Campbell, I.H.; Sossi, P.A.; Van Kranendonk, M.J.  Hafnium and iron isotopes in early Archean komatiites record a plume-driven convection cycle in the Hadean Earth.  Earth and Planetary Science Letters 2014, 397, 111–120.  https://doi.org/10.1016/j.epsl.2014.04.028.
  * Walton et al. (2023) Walton, C.R.; Rigley, J.K.; Lipp, A.; Law, R.; Suttle, M.D.; Schönbächler, M.; Wyatt, M.; Shorttle, O.  Cosmic dust fertilization of prebiotic chemistry on early Earth.  Nature Astronomy 2023, p. Accepted.
  * Rempfert et al. (2023) Rempfert, K.R.; Nothaft, D.B.; Kraus, E.A.; Asamoto, C.K.; Evans, R.D.; Spear, J.R.; Matter, J.M.; Kopf, S.H.; Templeton, A.S.  Subsurface biogeochemical cycling of nitrogen in the actively serpentinizing Samail Ophiolite, Oman.  Frontiers in Microbiology 2023, 14, 1139633.
  * Pastorek et al. (2020) Pastorek, A.; Ferus, M.; Čuba, V.; Šrámek, O.; Ivanek, O.; Civiš, S.  Primordial radioactivity and prebiotic chemical evolution: Effect of $\gamma$ radiation on formamide-based synthesis.  The Journal of Physical Chemistry B 2020, 124, 8951–8959.
  * Barge and Price (2022) Barge, L.M.; Price, R.E.  Diverse geochemical conditions for prebiotic chemistry in shallow-sea alkaline hydrothermal vents.  Nature Geoscience 2022, 15, 976–981.  https://doi.org/10.1038/s41561-022-01067-1.