# Soundness and completeness of the cirquent calculus system CL6 for
computability logic∗11footnotetext: This work was supported by the NNSF
(60974082) of China.

Wenyan Xu and Sanyang Liu
Department of Mathematics, Xidian University, Xi’an, 710071, PR China

###### Abstract

Computability logic is a formal theory of computability. The earlier article
“Introduction to cirquent calculus and abstract resource semantics” by
Japaridze proved soundness and completeness for the basic fragment CL5 of
computability logic. The present article extends that result to the more
expressive cirquent calculus system CL6, which is a conservative extension of
both CL5 and classical propositional logic.

Keywords: Cirquent calculus; Computability logic.

##  1 Introduction

Computability logic(CoL), introduced by G. Japaridze [1]-[3], is a semantical
and mathematical platform for redeveloping logic as a formal theory of
computability. Formulas in CoL represent interactive computational problems,
understood as games between a machine and its environment (symbolically named
as $\top$ and $\bot$, respectively); logical operators stand for operations on
such problems; “truth” of a problem/game means existence of an algorithmic
solution, i.e. $\top$’s effective winning strategy; and validity of a logical
formula is understood as such truth under every particular interpretation of
atoms. The approach induces a rich collection of (old or new) logical
operators. Among those, relevant to this paper are $\neg$ (negation), $\vee$
(parallel disjunction) and $\wedge$ (parallel conjunction). Intuitively,
$\neg$ is a role switch operator: $\neg A$ is the game $A$ with the roles of
$\top$ and $\bot$ interchanged ($\top$’s legal moves and wins become those of
$\bot$, and vice versa). Both $A\wedge B$ and $A\vee B$ are games playing
which means playing the two components $A$ and $B$ simultaneously (in
parallel). In $A\wedge B$, $\top$ is the winner if it wins in both components,
while in $A\vee B$ winning in just one component is sufficient. The symbols
$\top$ and $\bot$, together with denoting the two players, are also used to
denote two special (the simplest) sorts of games. Namely, $\top$ is a moveless
(“elementary”) game automatically won by the player $\top$, and $\bot$ is a
moveless game automatically won by $\bot$.

Cirquent calculus is a refinement of sequent calculus. Unlike the more
traditional proof theories that manipulate tree-like objects (formulas,
sequents, hypersequents, etc.), cirquent calculus deals with graph-style
structures termed cirquents, with its main characteristic feature thus being
allowing to explicitly account for sharing subcomponents between different
subcomponents. The approach was introduced by Japaridze [4] as a new deductive
tool for CoL and was developed later in [5]-[7]. The paper [4] constructed a
cirquent calculus system CL5 for the basic $(\neg,\wedge,\vee)$-fragment of
CoL, and proved its soundness and completeness with respect to the semantics
of CoL.

The atoms of CL5 represent computational problems in general, and are said to
be general atoms. The so called elementary atoms, representing computational
problems of zero degree of interactivity (such as the earlier-mentioned games
$\top$ and $\bot$) and studied in other pieces of literature on CoL, are not
among them. Thus, CL5 only describes valid computability principles for
general problems. This is a significant limitation of expressive power. For
example, the problem $A\rightarrow A{\wedge}A$ is not valid in CoL when $A$ is
a general atom, but becomes valid (as any classical tautology for that matter)
when $A$ is elementary. So the language of CL5 naturally calls for an
extension.

Japaridze [4] claimed without a proof that the soundness and completeness
result for CL5 could be extended to the more expressive cirquent calculus
system CL6 (reproduced later), which is a conservative extension of both CL5
and classical propositional logic. This article is devoted to a soundness and
completeness proof for system CL6, thus contributing to the task of extending
the cirquent-calculus approach so as to accommodate incrementally expressive
fragments of CoL.

##  2 Preliminaries

This paper primarily targets readers already familiar with Japaridze [4], and
can essentially be treated as a technical appendix to the latter. However, in
order to make it reasonably self-contained, in this section we reproduce the
basic concepts from [4] on which the later parts of the paper will rely. An
interested reader may consult [4] for additional explanations, illustrations
and examples.

The language of CL6 is more expressive than that of CL5 in that, along with
the old atoms of CL5 called general, it has an additional sort of atoms called
elementary, including non-logical elementary atoms and logical atoms $\top$
and $\bot$. On the other hand, all general atoms are non-logical. We use the
uppercase letters $P,Q,R,S$ as metavariables for general atoms, and the
lowercase $p,q,r,s$ as metavariables for non-logical elementary atoms. A
CL6-formula is built from atoms in the standard way using the connectives
$\neg$,$\vee$,$\wedge$, with $F\rightarrow G$ understood as an abbreviation
for $\neg F\vee G$ and $\neg$ limited only to non-logical atoms, where
$\neg\neg F$ is understood as $F$, $\neg(F\wedge G)$ as $\neg F\vee\neg G$,
$\neg(F\vee G)$ as $\neg F\wedge\neg G$, $\neg\top$ as $\bot$, and $\neg\bot$
as $\top$. An atom $P$ (resp. $p$) and its negation $\neg P$ (resp. $\neg p$)
is called a literal, and the two literals are said to be opposite. A
CL6-formula is said to be elementary iff it does not contain general atoms.
Throughout the rest of this paper, unless otherwise specified, by an “atom” or
a “formula” we mean one of the language of CL6.

Where $k\geq 0$, a $k-$ary pool is a sequence $\langle
F_{1},F_{2},\ldots,F_{k}\rangle$ of $k$ formulas. Since we may have
$F_{i}=F_{j}$ for some $i\neq j$ in such a sequence, we use the term oformula
to refer to a formula together with a particular occurrence of it in the pool.
For example, the pool $\langle E,F,G,E\rangle$ has three formulas but four
oformulas. Similarly, the terms “oliteral”,“oatom”, etc. will be used in this
paper to refer to the corresponding entities together with particular
occurrences. A $k-$ary structure is a finite sequence
St$=\langle\Gamma_{1},\ldots,\Gamma_{m}\rangle$, where $m\geq 0$ and each
$\Gamma_{i}$, said to be a group of St, is a subset of $\\{1,\ldots,k\\}$.
Again, to differentiate between a group as such and a particular occurrence of
a group in the structure, we use the term ogroup for the latter. For example,
the structure $\langle\\{2,3\\},\\{2,3\\},\\{1,4\\},\emptyset\rangle$ has
three groups but four ogroups.

A $k$-ary ($k\geq 0$) cirquent is a pair $C=({\bf St}^{C},{\bf Pl}^{C})$,
where ${\bf St}^{C}$, called the structure of $C$, is a $k$-ary structure, and
${\bf Pl}^{C}$, called the pool of $C$, is a $k$-ary pool. An ogroup of such a
$C$ will mean an ogroup of ${St}^{C}$, and an oformula of $C$ will mean an
oformula of ${\bf Pl}^{C}$. Usually, we understand the groups of a cirquent as
sets of its oformulas rather than sets of the corresponding ordinal numbers.
Thus, if ${\bf Pl}^{C}=\langle E,F,G,E\rangle$ and $\Gamma=\\{2,4\\}$, we
would think of $\Gamma$ simply as the set $\\{F,E\\}$, and say that $\Gamma$
contains $F$ and $E$. When both the pool and the structure of a cirquent $C$
are empty, i.e. $C=(\langle\rangle,\langle\rangle)$, we call it the empty
cirquent.

Rather than writing cirquents as ordered tuples in the above-described style,
we prefer to represent them through (and identify them with) diagrams. Below
is such a representation for the cirquent whose pool is $\langle
E,F,G,H\rangle$ and whose structure is
$\langle\\{1,2\\},\\{2\\},\\{3,4\\}\rangle$.

$E\ \ \ \ \ \ F\ \ \ \ \ G\ \ \ \ \ \ H$$\bullet$$\bullet$$\bullet$

The top level of a diagram thus indicates the oformulas of the cirquent, and
the bottom level gives its ogroups. An ogroup $\Gamma$ is represented by a
$\bullet$, and the lines connecting $\Gamma$ with oformulas, called arcs, are
pointing to the oformulas that $\Gamma$ contains. Finally, we put a horizontal
line at the top of the diagram to indicate that this is one cirquent rather
than two or more cirquents put together.

A model is a function $M$ that assigns a truth value — true (1) or false (0) —
to each atom, with $\top$ being always assigned true and $\bot$ false, and
extends to compound formulas in the standard classical way. Let $M$ be a
model, and $C$ a cirquent. We say that a group $\Gamma$ of $C$ is true in $M$
iff at least one of its oformulas is so. And $C$ is true in $M$ if every group
of $C$ is so. Otherwise, $C$ is false. Finally, $C$ or a group $\Gamma$ of it
is a tautology iff it is true in every model.

A substitution is a function $\sigma$ that sends every general atom $P$ to
some formula $\sigma(P)$, and sends every elementary atom to itself. If, (for
every general atom $P$), such a $\sigma(P)$ is an atom, then $\sigma$ is said
to be an atomic-level substitution.

Let $A$ and $B$ be cirquents. We say that $B$ is an instance of $A$ iff
$B=\sigma(A)$ for some substitution $\sigma$, where $\sigma(A)$ is the result
of replacing in all oformulas of $A$ every (general or elementary) atom
$\alpha$ by $\sigma(\alpha)$; and $B$ is an atomic-level instance of $A$ iff
$B=\sigma(A)$ for some atomic-level substitution $\sigma$.

A cirquent is said to be binary iff no general atom has more than two
occurrences in it. A binary cirquent is said to be normal iff, whenever it has
two occurrences of a general atom, one occurrence is negative and the other is
positive. A binary tautology (resp. normal binary tautology) is a binary
(resp. normal binary) cirquent that is a tautology.

The set of rules of CL6 is obtained from that of CL5 by adding to it $\top$ as
an additional axiom, plus the rule of contraction limited only to elementary
formulas. Below we reproduce those rules from [4], followed by illustrations.

Axioms (A): Axioms are “rules” with no premises. There are three sorts of
axioms in CL6. The first one is the empty cirquent. The second one is any
cirquent that has exactly two oformulas $F$ and $\neg F$, for some arbitrary
formula $F$, and an ogroup that contains $F$ and $\neg F$. In other words,
this is the cirquent $(\langle\\{1,2\\}\rangle,\langle F,\neg F\rangle)$. The
third one is a cirquent that has exactly one oformula $\top$ and one ogroup
that contains $\top$, i.e. the cirquent
$(\langle\\{1\\}\rangle,\langle\top\rangle)$.

Mix (M): According to this rule, the conclusion can be obtained by simply
putting any two cirquents (premises) together, thus creating one cirquent out
of two.

Exchange (E): This rule comes in two versions: oformula exchange and ogroup
exchange. The conclusion of oformula exchange is obtained by interchanging in
the premise two adjacent oformulas $E$ and $F$, and redirecting to $E$ (resp.
$F$) all arcs that were originally pointing to $E$ (resp. $F$). Ogroup
exchange is the same, with the only difference that the objects interchanged
are ogroups.

Weakening (W): This rule also comes in two versions: ogroup weakening and pool
weakening. A conclusion of ogroup weakening is obtained by adding in the
premise a new arc between an existing ogroup and an existing oformula. As for
pool weakening, a conclusion is obtained through inserting a new oformula
anywhere in the pool of the premise.

Duplication (D): A conclusion of this rule is obtained by replacing in the
premise some ogroup $\Gamma$ by two adjacent ogroups that, as groups, are
identical with $\Gamma$.

Contraction (C): According to this rule, if a cirquent (a premise) has two
adjacent elementary oformulas $F$ (the first), $F$ (the second) that are
identical, then a conclusion can be obtained by merging $F$,$F$ into $F$ and
redirecting to the latter all arcs that were originally pointing to the first
or the second $F$.

$\vee-$introduction ($\vee$): For the convenience of description, we explain
this rule in the bottom-up view. According to this rule, if a cirquent (the
conclusion) has an oformula $E\vee F$ that is contained by at least one
ogroup, then the premise can be obtained by splitting the original $E\vee F$
into two adjacent oformulas $E$ and $F$, and redirecting to both $E$ and $F$
all arcs that were originally pointing to $A\vee B$.

$\wedge-$introduction ($\wedge$): This rule, again, is more conveniently
described in the bottom-up view. According to this rule, if a cirquent (the
conclusion) has an oformula $E\wedge F$ that is contained by at least one
ogroup, then the premise can be obtained by splitting the original $E\wedge F$
into two adjacent oformulas $E$ and $F$, and splitting every ogroup $\Gamma$
that originally contained $E\wedge F$ into two adjacent ogroups $\Gamma^{E}$
and $\Gamma^{F}$, where $\Gamma^{E}$ contains $E$ (but not $F$), and
$\Gamma^{F}$ contains $F$ (but not $E$), with all other ($\neq E\wedge F$)
oformulas of $\Gamma$ contained by both $\Gamma^{E}$ and $\Gamma^{F}$.

Below we provide illustrations for all rules, in each case an abbreviated name
of the rule standing next to the horizontal line separating the premises from
the conclusions. Our illustrations for the axioms (the “A” labeled rules) are
specific cirquents or schemate of such; our illustrations for all other rules
are merely examples chosen arbitrarily. Unfortunately, no systematic ways for
schematically representing cirquent calculus rules have been elaborated so
far. This explains why we appeal to examples instead.

A$F\ \ \ \ \ \ \ \neg F$A$\bullet$$\top$A$\bullet$

$E$$\bullet$M$F\ \ \ \ \ G$$\bullet$$E$$\bullet$$F\ \ \ \ \
G$$\bullet$oformula exchange$E\ \ \ \ \ F\ \ \ \
G$$\bullet$$\bullet$$\bullet$E$F\ \ \ \ \ E\ \ \ \
G$$\bullet$$\bullet$$\bullet$ogroup exchange$E\ \ \ \ \ F\ \ \ \
G$$\bullet$$\bullet$$\bullet$E$E\ \ \ \ \ F\ \ \ \
G$$\bullet$$\bullet$$\bullet$

ogroup weakening$F\ \ \ \ G\ \ \ \ \ H$$\bullet$$\bullet$$\bullet$W$F\ \ \ \
G\ \ \ \ \ H$$\bullet$$\bullet$$\bullet$pool weakening$F\ \ \ \ \ \ \ \ \ \
H$$\bullet$$\bullet$W$F\ \ \ \ G\ \ \ \ H$$\bullet$$\bullet$$E\ \ \ \ \ F\ \ \
\ G$$\bullet$$\bullet$D$E\ \ \ \ \ F\ \ \ \ G$$\bullet$$\bullet$$\bullet$

$F$required to beelementary$E\ \ \ \ \ F\ \ \ \ F$$\bullet$$\bullet$C$E\ \ \ \
\ \ \ \ \ \ F$$\bullet$$\bullet$$H\ \ E\ \ F\ \ H$$\bullet$$\bullet$$\vee$$H\
\ E\vee F\ \ H$$\bullet$$\bullet$$F\ \ \ E\ \ \ F\ \
G$$\bullet$$\bullet$$\bullet$$\bullet$$\bullet$$\wedge$$F\ \ E\wedge F\ \
G$$\bullet$$\bullet$$\bullet$

The above are all eight rules of CL6. As a warm-up exercise, the reader may
try to verify that CL6 proves $p\rightarrow p\wedge p$ but does not prove
$P\rightarrow P\wedge P$.

As an aside, the earlier mentioned system CL5 differs from CL6 in that the
$\top-$axiom and the contraction rules are absent there. Also, as noted, the
language of CL5 does not allow elementary atoms. In next section we will see
that our proofs are carried out purely syntactically, based on the soundness
and completeness of system CL2 (introduced in Japaridze [8]) with respect to
the semantics of CoL. That is to say we do not directly use the semantics of
CoL. So, below we only explain what the language of CL2 and its rules are,
without providing any formal definitions (on top of the brief informal
explanations given in Section 1) of the underlying CoL semantics. If
necessary, such definitions can be found in [3].

The language of CL2 is more expressive than the one in which formulas of CL6
are written because, on top of $\neg$,$\vee$,$\wedge$, it has the binary
connectives $\sqcap$ and $\sqcup$, called choice operators. The CL2-formulas
are built from atoms (including general atoms and elementary atoms) in the
standard way using the connectives $\neg$,$\vee$,$\wedge$,$\sqcap$,$\sqcup$.
As in the case of CL6-formulas, the operator $\neg$ is only allowed to be
applied to non-logical atoms. A CL2-formula is said to be elementary iff it
contains neither general atoms nor $\sqcap$,$\sqcup$. A positive occurrence
(resp. negative occurrence) of an atom is one that is not (resp. is) in the
scope of $\neg$. A surface occurrence of a subformula of a CL2-formula is an
occurrence that is not in the scope of $\sqcap$,$\sqcup$. A general literal is
$P$ or $\neg P$, where $P$ is a general atom. The elementarization of a
CL2-formula $A$ is the result of replacing in $A$ every positive surface
occurrence of each general literal by $\bot$, every surface occurrence of each
$\sqcup-$subformula by $\bot$, and every surface occurrence of each
$\sqcap-$subformula by $\top$. A CL2-formula is said to be stable iff its
elementarization is a tautology of classical logic.

CL2 has the following three inference rules.

Rule (a): $\overrightarrow{H}\mapsto F$, where $F$ is stable and
$\overrightarrow{H}$ is the smallest set of formulas such that, whenever $F$
has a surface occurrence of a subformula $G_{1}\sqcap G_{2}$, for both
$i$$\in\\{$1,2$\\}$, $\overrightarrow{H}$ contains the result of replacing
that occurrence in $F$ by $G_{i}$.

Rule (b): $H\mapsto F$, where $H$ is the result of replacing in $F$ a surface
occurrence of a subformula $G_{1}\sqcup G_{2}$ by $G_{1}$ or $G_{2}$.

Rule (c): $H\mapsto F$, where $H$ is the result of replacing in $F$ two — one
positive and one negative — surface occurrences of some general atom by a non-
logical elementary atom that does not occur in $F$.

The set $\overrightarrow{H}$ of the premises of Rule (a) may be empty, in
which case the rule (its conclusion, that is) acts like an axiom. Otherwise,
the system has no (other) axioms.

##  3 Soundness and completeness of CL6

In what follows, we may use names such as (AME) to refer to the subsystem of
CL6 consisting only of the rules whose names are listed between the
parentheses. So, (AME) refers to the system that only has axioms, exchange and
mix. The same notation can be used next to the horizontal line separating two
cirquents to indicate that the lower cirquent (“conclusion”) can be obtained
from the upper cirquent (“premise”) by whatever number of applications of the
corresponding rules. The following Lemmas 1, 2, 3, 4 are precisely Lemmas 4,
5, 10 and 11 of [4], so we state them without proofs (such proofs are given in
[4]).

###### Lemma 1

All of the rules of CL6 preserve truth in the top-down direction. Taking no
premises, (the conclusion of) axioms are thus tautologies.

###### Lemma 2

The rules of mix, exchange, duplication, contraction, $\vee$-introduction and
$\wedge$-introduction preserve truth in the bottom-up direction as well.

###### Lemma 3

The rules of mix, exchange, duplication, $\vee$-introduction and
$\wedge$-introduction preserve binarity and normal binarity in both top-down
and bottom-up directions.

###### Lemma 4

Weakening preserves binarity and normal binarity in the bottom-up direction.

###### Lemma 5

If CL6 proves a cirquent $C$, then it also proves every instance of $C$.

Proof. Let $T$ be a proof tree of an arbitrary cirquent $C$, $C^{\prime}$ be
an arbitrary instance of $C$, and $\sigma$ be a substitution with
$\sigma(C)=C^{\prime}$. Replace every oformula $F$ of every cirquent of $T$ by
$\sigma(F)$. It is not hard to see that the resulting tree $T^{\prime}$, which
uses exactly the same rules as $T$ does, is a proof of $C^{\prime}$.

###### Lemma 6

Contraction preserves binarity and normal binarity in both top-down and
bottom-up directions.

Proof. This is so because contraction limited to elementary formulas can never
affect what general atoms occur in a cirquent and how many times they occur.

###### Lemma 7

A cirquent is provable in CL6 iff it is an instance of a binary tautology.

Proof. $(\Rightarrow)$ Consider an arbitrary cirquent $A$ provable in CL6. By
induction on the height of its proof tree, we want to show that $A$ is an
instance of a binary tautology.

The above is obvious when $A$ is an axiom.

Suppose now $A$ is derived by exchange from $B$. Let us just consider oformula
exchange, with ogroup exchange being similar. By the induction hypothesis, $B$
is an instance of a binary tautology $B^{\prime}$. Let $A^{\prime}$ be the
result of applying exchange to $B^{\prime}$ “at the same place” as it was
applied to $B$ when deriving $A$ from it, as illustrated in the following
example:

$P\ \ \ s\ \ \ \neg P\ \ \ \ P\vee r\ \ \ \neg P\wedge\neg r$${\bf E}$$P\ \
\neg P\ \ \ s\ \ \ \ \ P\vee r\ \ \ \neg P\wedge\neg r$$B:$$A:$$E\ \ \ s\ \ \
\neg E\ \ \ R\vee r\ \ \ \ \neg R\wedge\neg r$${\bf E}$$E\ \ \neg E\ \ \ s\ \
\ \ R\vee r\ \ \ \ \neg R\wedge\neg r$$B^{\prime}:$$A^{\prime}:$

Obviously $A$ will be an instance of $A^{\prime}$. It remains to note that, by
Lemmas 1 and 3, $A^{\prime}$ is a binary tautology.

The rules of duplication, $\vee$-introduction and $\wedge$-introduction can be
handled in a similar way.

Next, suppose $A$ is derived from $B$ and $C$ by mix. By the induction
hypothesis, $B$ and $C$ are instances of some binary tautologies $B^{\prime}$
and $C^{\prime}$, respectively. We may assume that no general atom $P$ occurs
in both $B^{\prime}$ and $C^{\prime}$, for otherwise, in one of the cirquents,
rename $P$ into another general atom $Q$ different from everything else. Let
$A^{\prime}$ be the result of applying mix to $B^{\prime}$ and $C^{\prime}$.
By Lemmas 1 and 3, $A^{\prime}$ is a binary tautology. And, as in the cases of
the other rules, it is evident that $A$ is an instance of $A^{\prime}$.

Suppose $A$ is derived from $B$ by weakening. If this is ogroup weakening, $A$
is an instance of a binary tautology for the same reason as in the case of
exchange, duplication, $\vee$-introduction or $\wedge$-introduction. Assume
now we are dealing with pool weakening, so that $A$ is the result of inserting
a new oformula $F$ into $B$. By the induction hypothesis, $B$ is an instance
of a binary tautology $B^{\prime}$. Let $P$ be a general atom not occurring in
$B^{\prime}$. And let $A^{\prime}$ be the result of applying weakening to
$B^{\prime}$ that inserts $P$ “at the same place” into $B^{\prime}$ as the
above application of weakening inserted $F$ into $B$ when deriving $A$.
Obviously $A^{\prime}$ inherits binarity from $B^{\prime}$; by Lemma 1, it
inherits from $B^{\prime}$ tautologicity as well. And, for the same reason as
in all previous cases, $A$ is an instance of $A^{\prime}$.

Finally, suppose $A$ is derived from $B$ by contraction. Then the contracted
formula $F$ should be elementary. By the induction hypothesis, $B$ is an
instance of a binary tautology $B^{\prime}$. Let $\sigma$ be a substitution
such that $B=\sigma(B^{\prime})$. And let $F^{\prime}_{1}$, $F^{\prime}_{2}$
be two oformulas in $B^{\prime}$ “at the same place” as $F$, $F$ are in $B$,
with $\sigma(F^{\prime}_{1})=F$ and $\sigma(F^{\prime}_{2})=F$. Let $\delta$
be the substitution such that, for any general atom $P$, $\delta(P)=\sigma(P)$
if $P$ occurs in $F^{\prime}_{1}$ or $F^{\prime}_{2}$, and $\delta(P)=P$
otherwise. Thus, $\delta(F^{\prime}_{1})=\delta(F^{\prime}_{2})=F$. And let
$B^{\prime\prime}=\delta(B^{\prime})$. Obviously — for the same reasons as in
classical logic — substitution does not destroy tautologicity, so
$B^{\prime\prime}$ is a tautology because $B^{\prime}$ is so. Further, the
substitution $\delta$ does not introduce any new occurrences of general atoms,
so it does not destroy the binarity of $B^{\prime}$, either. To summarize,
$B^{\prime\prime}$ is a binary tautology. Also, of course, $B$ is an instance
of $B^{\prime\prime}$. Notice that $B^{\prime\prime}$ has $F$ and $F$ where
$B$ has the contracted oformulas $F$ and $F$. So, let $A^{\prime}$ be the
result of applying contraction to $B^{\prime\prime}$ “at the same place” as it
was applied to $B$ when deriving $A$ from it, as illustrated in the following
example:

$P\ \ r\wedge s\ \ \ r\wedge s\ \ \ \ \ \neg P\ \ \ \ P\vee q\ \ \ \neg
P\wedge\neg q$${\bf C}$$P\ \ \ \ \ \ \ \ \ \ r\wedge s\ \ \ \ \ \ \neg P\ \ \
\ P\vee q\ \ \ \neg P\wedge\neg q$$B:$$A:$$R\ \ r\wedge s\ \ \ r\wedge s\ \ \
\ \neg R\ \ \ \ \ Q\vee q\ \ \ \neg Q\wedge\neg q$${\bf C}$$R\ \ \ \ \ \ \ \ \
\ r\wedge s\ \ \ \ \ \ \neg R\ \ \ \ Q\vee q\ \ \ \neg Q\wedge\neg
q$$B^{\prime\prime}:$$A^{\prime}:$

Obviously $A$ will be an instance of $A^{\prime}$. And, by Lemma 1 and Lemma
6, $A^{\prime}$ is a binary tautology.

$(\Leftarrow)$ Consider an arbitrary cirquent $A$ that is an instance of a
binary tautology $A^{\prime}$. In view of Lemma 5, it would suffice to show
that CL6 proves $A^{\prime}$. We construct a proof of $A^{\prime}$, in the
bottom-up fashion, as follows. Starting from $A^{\prime}$, we keep applying
$\vee$-introduction and $\wedge$-introduction until we hit an essentially
literal cirquent111An essentially literal cirquent, defined in [4], is one
every oformula of whose pool either is an oliteral or is homeless. $B$. As in
the proof of Theorem 6 of [4], such a cirquent $B$ is guaranteed to be a
tautology, and $A^{\prime}$ follows from it in $(\vee\wedge)$. Furthermore, in
view of Lemma 3, $B$ is in fact a binary tautology. The tautologicity of $B$
means that every ogroup of it contains either a $\top$, or at least one pair
of opposite (general or elementary) non-logical oliterals. For each ogroup of
$B$ that contains a $\top$, pick one occurrence of $\top$ and apply to $B$ a
series of weakenings to first delete all arcs but the arc pointing to the
chosen occurrence, and next delete all homeless oformulas if any such
oformulas are present. For each ogroup of the resulting cirquent that contains
a pair of opposite non-logical oliterals, pick one such pair, and continue
applying a series of weakenings, as in the proof of Theorem 6 of [4], until a
tautological cirquent $C$ is hit with no homeless oformulas, where every
ogroup only has either a $\top$ or a pair of opposite non-logical oliterals.
By Lemma 4, $C$ remains binary. Our target cirquent $A^{\prime}$ is thus
derivable from $C$ in (W$\vee\wedge$). Apply a series of contractions to $C$
to separate all shared $\top$ and all shared elementary non-logical oliterals
$p$ or $\neg p$, as illustrated below; as a result, we get a cirquent $D$
which is still a binary tautology, but whose ogroups no longer share any
elementary oformulas.

$P\ \ \ \ \neg P\ \ \ r\ \ \ \ \neg r\ \ \neg r\ \ \ \ \ r\ \ \ \ r\ \ \ \neg
Q\ \ \neg r\ \ \neg s\ \ \neg s\ \ \ Q\ \ \ s\ \ \ \ s\ \ \top\ \ \ \top\ \ \
\top$(C)$P\ \ \ \ \neg P\ \ \ r\ \ \ \ \neg r\ \ \ \ \ r\ \ \ \neg Q\ \ \neg
r\ \ \ \neg s\ \ \ \ Q\ \ \ s\ \ \ \top\ \ \ \ \ \ \top$$D:$$C:$

It is easy to see that the binarity of $D$ implies that there are no shared
general oliterals $P$ or $\neg P$ in it except the cases when they are shared
by identical-content ogroups. Applying to $D$ a series of duplications, as
illustrated below, yields a cirquent $E$ that no longer has identical-content
ogroups and hence no longer has any shared oformulas.

$P\ \ \ \ \neg P\ \ \ r\ \ \ \ \neg r\ \ \neg r\ \ \ \ \ r\ \ \ \ r\ \ \ \neg
Q\ \ \neg r\ \ \neg s\ \ \neg s\ \ \ Q\ \ \ s\ \ \ \ s\ \ \top\ \ \ \top\ \ \
\top$(D)$E:$$P\ \ \ \ \neg P\ \ \ r\ \ \ \ \neg r\ \ \neg r\ \ \ \ \ r\ \ \ \
r\ \ \ \neg Q\ \ \neg r\ \ \neg s\ \ \neg s\ \ \ Q\ \ \ s\ \ \ \ s\ \ \top\ \
\ \top\ \ \ \top$$D:$

$A^{\prime}$ is thus derivable from $E$ in (DCW$\vee\wedge$). In turn, $E$ is
obviously provable in (AME). So, CL6 proves $A^{\prime}$.

###### Lemma 8

A cirquent is an instance of a binary tautology iff it is an atomic-level
instance of some normal binary tautology.

Proof. Our proof here almost literally follows the proof of Lemma 9 of [4].

The “if” part is trivial. For the “only if” part, assume $A$ is an instance of
a binary tautology $B$. Let $P_{1},\ldots,P_{n}$ be all of the general atoms
of $B$ that have two positive or two negative occurrences in $B$. Let
$Q_{1},\ldots,Q_{n}$ be any pairwise distinct general atoms not occurring in
$B$. Let $C$ be the result of replacing in $B$ one of the two occurrences of
$P_{i}$ by $Q_{i}$, for each $i=1,\ldots,n$. Then obviously $C$ is a normal
binary cirquent, and $B$ an instance of it. By transitivity, $A$ (as an
instance of $B$) is also an instance of $C$.

We want to see that $C$ is a tautology. Deny this. Then there is a classical
model $M$ in which $C$ is false. Let $M^{\prime}$ be the model such that:

  * •

$M^{\prime}$ agrees with $M$ on all atoms that are not among
$P_{1},\ldots,P_{n},Q_{1},\ldots,Q_{n}$;

  * •

for each $i\in\\{1,\ldots,n\\}$, $M^{\prime}(P_{i})=M^{\prime}(Q_{i})=$false
if $P_{i}$ and $Q_{i}$ are positive in $C$; and
$M^{\prime}(P_{i})=M^{\prime}(Q_{i})=$true if $P_{i}$ and $Q_{i}$ are negative
in $C$.

By induction on complexity, it can be easily seen that, for every subformula
$F$ of a formula of $C$, whenever $F$ is false in $M$, so is it in
$M^{\prime}$. This extends from (sub)formulas to groups of $C$ and hence $C$
itself. Thus $C$ is false in $M^{\prime}$ because it is false in $M$. But
$M^{\prime}$ does not distinguish between $P_{i}$ and $Q_{i}$ (any $1\leq
i\leq n$). This clearly implies that $C$ and $B$ have the same truth value in
$M^{\prime}$. That is, $B$ is false in $M^{\prime}$, which is however
impossible because $B$ is a tautology. From this contradiction we conclude
that $C$ is a (normal binary) tautology.

Let $\sigma$ be a substitution such that $A=\sigma(C)$. Let $\sigma^{\prime}$
be a substitution such that, for each general atom $P$ of $C$,
$\sigma^{\prime}(P)$ is the result of replacing in $\sigma(P)$ each occurrence
of each general atom by a new general atom in such a way that: no general atom
occurs more than once in $\sigma^{\prime}(P)$, and whenever $P\neq Q$, no
general atom occurs in both $\sigma^{\prime}(P)$ and $\sigma^{\prime}(Q)$.
Since $C$ is a binary tautology and is its own instance, by Lemma 7, CL6
proves $C$. Then, by Lemma 5, CL6 proves $\sigma^{\prime}(C)$ (an instance of
$C$). In view of Lemma 1, we immediately get that $\sigma^{\prime}(C)$ is a
tautology. $\sigma^{\prime}(C)$ can also be easily seen to be a normal binary
cirquent, because $C$ is so. Finally, with a little thought, $A$ can be seen
to be an atomic-level instance of $\sigma^{\prime}(C)$.

###### Lemma 9

A CL6-formula $F$ is provable in CL2 iff it is an instance of a binary
tautology.

Proof. Again, it should be acknowledged that the present proof very closely
follows the proof of Lemma 27 of [4], even though there are certain
differences.

$(\Rightarrow)$ Consider an arbitrary CL6-formula $F$ provable in CL2. Fix a
CL2-proof of $F$ in the form of a sequence $\langle
F_{n},F_{n-1},\ldots,F_{1}\rangle$ of formulas, with $F_{1}=F$. We may assume
that this sequence has no repetitions or other redundancies. We claim that,
for each $i$ with $1\leq i\leq n$, the following conditions are satisfied:

Condition 1: $F_{i}$ does not contain $\sqcap,\sqcup$.

Condition 2: Whenever $F_{i}$ contains an elementary atom not occurring in
$F$, that atom is non-logical, and has exactly two — one positive and one
negative — occurrences in $F_{i}$.

Condition 3: If $i<n$, then $F_{i}$ is derived from $F_{i+1}$ by Rule (c).

Condition 4: $F_{n}$ is derived (from the empty set of premises) by Rule (a).

Condition 4 is obvious, because it is only Rule (a) that may take no premises.
That Conditions 1-3 are also satisfied can be verified by induction on $i$.
For the basis case of $i=1$, Conditions 1 and 2 are immediate. $F_{1}$ can not
be derived by Rule (b) because, by Condition 1, $F_{1}$ does not contain any
$\sqcup$. Nor can it be derived by Rule (a) unless $n=1$, for otherwise either
$F_{1}$ would have to contain a $\sqcap$ (which is not the case according to
Condition $1$), or the proof of $F$ would have redundancies as $F_{1}$ would
not really need any premises. Thus, if $1<n$, the only possibility for $F_{1}$
is to be derived from $F_{2}$ by Rule (c). For the induction step, assume
$i<n$ and the above conditions are satisfied for $F_{i}$. According to
Condition 3, $F_{i}$ is derived by Rule (c) from $F_{i+1}$ . This obviously
implies that $F_{i+1}$ inherits Conditions 1 and 2 from $F_{i}$. And that
Condition 3 also holds for $F_{i+1}$ can be shown in the same way as we did
for $F_{1}$.

As the conclusion of Rule (a) (Condition 4), $F_{n}$ is stable. Let $G$ be the
elementarization of $F_{n}$. The stability of $F_{n}$ means that $G$ is a
tautology. Let $H$ be the result of replacing in $G$ every occurrence of
$\top$ and $\bot$ (except those inherited from $F$) by a general atom, in such
a way that different occurrences of $\top$, $\bot$ are replaced by different
atoms. In view of Condition 2 (applied to $F_{n}$), we see that, on top of
these new general atoms and the elementary atoms inherited from $F$, the only
additional atoms that $H$ contains are elementary atoms with exactly two — one
positive and one negative — occurrences. Let $H^{\prime}$ be the result of
replacing in $H$ every occurrence of every such elementary atom by a general
atom not occurring in $H$, in such a way that different elementary atoms are
replaced by different general atoms. Then it is not hard to see that
$H^{\prime}$ is binary and $F_{n}$ is an instance of $H^{\prime}$. With
Condition 3 in mind, by induction, one can further see that the formulas
$F_{n-1}$, $F_{n-2}$, $\ldots$ are also instances of $H^{\prime}$. Thus, $F$
is an instance of $H^{\prime}$. It remains to show that $H^{\prime}$ is a
tautology. But this is indeed so because $H^{\prime}$ results from the
tautological $G$ by replacing positive occurrences of $\bot$ and replacing two
— one positive and one negative — occurrences of elementary atoms by general
atoms. It is known from classical logic that such replacements do not destroy
truth and hence tautologicity of formulas.

$(\Leftarrow)$ Assume $F$ is a CL6-formula which is an instance of a binary
tautology $T$. In view of Lemma 8, we may assume that $T$ is normal and $F$ is
an atomic-level instance of it. Let us call the general atoms that only have
one occurrence in $T$ single, and the general atoms that have two occurrences
married. Let $\sigma$ be the substitution with $\sigma(T)=F$. Let $G$ be the
formula resulting from $T$ by the following steps: substituting each single
atom $P$ by $\sigma(P)$; substituting each married atom $Q$ by $\sigma(Q)$ if
$\sigma(Q)$ is elementary; substituting each married atom $R$ by a non-logical
elementary atom $r$ not occurring in $F$ if $\sigma(R)$ is general. It is
clear that then $F$ can be derived from $G$ by a series of applications of
Rule (c), with each such application replacing two — a positive and a negative
— occurrences of some non-logical elementary atom $r$ by $\sigma(R)$. So, in
order to show that CL2 proves $F$, it would suffice to verify that $G$ is
stable and hence it can be derived from the empty set of premises by Rule (a).
But $G$ is indeed stable. To see this, consider the elementarization
$G^{\prime}$ of $G$. It results from $T$ by replacing the only occurrence of
each single general atom by some elementary atom, and doing the same with both
occurrences of each married general atom. In other words, $G^{\prime}$ is an
instance of $T$. Hence, as $T$ is a tautology, so is $G^{\prime}$, meaning
that $G$ is stable.

###### Theorem 10

A formula is provable in CL6 iff it is valid in computability logic.

Proof. This theorem is an immediate corollary of Lemma 9, Lemma 7 and the
known fact (proven in [8]) that CL2 is sound and complete with respect to the
semantics of computability logic.

## References

  * [1] G.Japaridze, Introduction to computability logic, Annals of Pure and Applied Logic 123(2003)1-99.
  * [2] G.Japaridze, Computability logic: a formal theory of interaction, In: D.Goldin, S.Smolka, P.Wegner(Eds.), Interactive Computation: The New Paradigm, Springer Verlag, Berlin, 2006, pp. 183-223.
  * [3] G.Japaridze, In the beginning was game semantics, In: O.Majer, A.-V.Pietarinen, T. Tulenheimo(Eds.), Games: Unifying Logic, Language, and Philosophy, Springer Netherlands, Dordrecht, 2009, pp. 249-350.
  * [4] G.Japaridze, Introduction to cirquent calculus and abstract resource semantics, Journal of Logic and Computation 16(4)(2006)489-532.
  * [5] G.Japaridze, Cirquent calculus deepened, Journal of Logic and Computation 18(6)(2008)983-1028.
  * [6] G.Japaridze, From formulas to cirquents in computability logic, Logical Methods in Computer Science 7(2)(2011)1-55.
  * [7] Wen-yan Xu, San-yang Liu, Deduction theorem for symmetric cirquent calculus, Quantitative Logic and Soft Computing 2010, Advances in Intelligent and Soft Computing, 82(2010)121-126.
  * [8] G.Japaridze, Propositional computability logic II, ACM Transactions on Computational Logic 7(2)(2006)331-362.


# H2 formation on interstellar grains and the fate of reaction energy

Stefano Pantaleone present address: Dipartimento di Chimica and Nanostructured
Interfaces and Surfaces (NIS) Centre, Università degli Studi di Torino, via P.
Giuria 7, 10125, Torino, Italy.  and  Dipartimento di Chimica, Biologia e
Biotecnologie, Università di Perugia, Via Elce di Sotto 8, 06123 Perugia,
Italy  Univ. Grenoble Alpes, CNRS, Institut de Planétologie et d’Astrophysique
de Grenoble (IPAG), 38000 Grenoble, France  Joan Enrique-Romero Univ. Grenoble
Alpes, CNRS, Institut de Planétologie et d’Astrophysique de Grenoble (IPAG),
38000 Grenoble, France  Departament de Química, Universitat Autònoma de
Barcelona, Bellaterra, 08193, Catalonia, Spain  Cecilia Ceccarelli Univ.
Grenoble Alpes, CNRS, Institut de Planétologie et d’Astrophysique de Grenoble
(IPAG), 38000 Grenoble, France  Stefano Ferrero Departament de Química,
Universitat Autònoma de Barcelona, Bellaterra, 08193, Catalonia, Spain  Nadia
Balucani Univ. Grenoble Alpes, CNRS, Institut de Planétologie et
d’Astrophysique de Grenoble (IPAG), 38000 Grenoble, France  Dipartimento di
Chimica, Biologia e Biotecnologie, Università di Perugia, Via Elce di Sotto 8,
06123 Perugia, Italy  Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5,
50125 Firenze, Italy  Albert Rimola Departament de Química, Universitat
Autònoma de Barcelona, Bellaterra, 08193, Catalonia, Spain  Piero Ugliengo
Dipartimento di Chimica and Nanostructured Interfaces and Surfaces (NIS)
Centre, Università degli Studi di Torino, via P. Giuria 7, 10125, Torino,
Italy.

###### Abstract

Molecular hydrogen is the most abundant molecular species in the Universe.
While no doubts exist that it is mainly formed on the interstellar dust grain
surfaces, many details of this process remain poorly known. In this work, we
focus on the fate of the energy released by the H2 formation on the dust icy
mantles, how it is partitioned between the substrate and the newly formed H2,
a process that has a profound impact on the interstellar medium. We carried
out state-of-art ab-initio molecular dynamics simulations of H2 formation on
periodic crystalline and amorphous ice surface models. Our calculations show
that up to two thirds of the energy liberated in the reaction ($\sim$300
kJ/mol $\sim$3.1 eV) is absorbed by the ice in less than 1 ps. The remaining
energy ($\sim$140 kJ/mol $\sim$1.5 eV) is kept by the newly born H2. Since it
is ten times larger than the H2 binding energy on the ice, the new H2 molecule
will eventually be released into the gas-phase. The ice water molecules within
$\sim$4 Å from the reaction site acquire enough energy, between 3 and 14
kJ/mol (360–1560 K), to potentially liberate other frozen H2 and, perhaps,
frozen CO molecules. If confirmed, the latter process would solve the long
standing conundrum of the presence of gaseous CO in molecular clouds. Finally,
the vibrational state of the newly formed H2 drops from highly excited states
($\nu=6$) to low ($\nu\leq 2$) vibrational levels in a timescale of the order
of ps.

H2 formation, ice mantles, chemical desorption, energy dissipation, molecular
dynamics, DFT

††journal: ApJ

##  1 Introduction

Molecular hydrogen is the most abundant molecule in the Universe. Its
formation is also the first step of the interstellar chemistry and, therefore,
a fundamental reaction. In molecular clouds, H2 is mainly formed via the H + H
association reaction on the interstellar dust grain surfaces, which act as a
third body capable to partially absorb the energy ($\sim$440 kJ/mol $\sim$4.56
eV) released by the chemical reaction and, consequently, stabilise the newly
formed H2 molecule (e.g. van de Hulst, 1946; Hollenbach & Salpeter, 1970,
1971; Smoluchowski, 1983). Although no doubts exist on the occurrence of this
process, many specific details remain poorly known (see e.g. Vidali, 2013;
Wakelam et al., 2017).

Here, we focus on the fate of the energy released by the reaction, which has
been source of debate for decades. Specifically: how much of the reaction
energy is absorbed by the dust grain and in what timescale? What fraction of
the chemical energy does go into the kinetic energy of the newly formed H2? Is
this energy large enough as H2 breaks the interaction with the surface and
leaves into the gas-phase? How much energy does the H2 molecule possess when
it leaves the grain surface and in what ro-vibrational state? And, finally, is
the energy transmitted to the dust grain enough for locally warming it up and
make adjacent frozen molecules sublimate? These points have a profound impact
on several aspects of astrochemistry and, more generally, the physics and
chemistry of the interstellar medium (ISM) (Duley & Williams, 1993).

Answers to the above questions largely depend on the nature of the substrate,
namely the specific dust grain surface or, in practice, the interstellar
environment where H2 forms. Here we focus on cold ($\sim 10$ K) molecular
clouds. In these environments, the grain refractory cores are coated by water-
dominated icy mantles either in polycrystalline (PCI) or amorphous (ASW for
amorphous solid water) phases, the latter dominating over the former (Smith et
al., 1988; Whittet, 1993; Boogert et al., 2015).

Unfortunately, laboratory experiments cannot entirely reproduce the
interstellar conditions so that, while they can suggest routes and processes,
they cannot provide definitive answers to the above-exposed questions (see
e.g. Vidali, 2013). For example, obtaining experimentally how the H2 formation
nascent energy is partitioned is very difficult (Roser et al., 2003; Hornekær
et al., 2003; Hama et al., 2012; Watanabe et al., 2010). Computational
chemistry methods can be regarded as a complementary tool and, sometimes, as a
unique alternative to laboratory experiments. So far, a limited number of
studies devoted to the energy dissipation of the H2 formation reaction on icy
grains has appeared in the literature (Takahashi et al., 1999a; Takahashi &
Uehara, 2001; Herbst & Cuppen, 2006).

Here, we present a new theoretical study on the dissipation of the energy
released by the $H+H\longrightarrow H_{2}$ reaction on water ice using ab-
initio molecular dynamics simulations (AIMDs). We simulated the reaction
adopting a Langmuir-Hinshelwood (LH) mechanism (i.e. both reactants are
adsorbed one the surface in neighboring sites) both on amorphous and
crystalline periodic models of interstellar icy mantles.

##  2 Computational details

###  2.1 Methods

All our calculations were carried out with the CP2K package (Hutter et al.,
2014). Core electrons were described with the Goedecker-Teter-Hutter
pseudopotentials (Goedecker et al., 1996) and valence ones with a mixed
Gaussian and Plane Wave (GPW) approach (Lippert et al., 1997). The DFT
(Density Functional Theory) PBE (Perdew-Burke-Ernzerhof) method was adopted
(Perdew et al., 1996), combined with a triple-$\zeta$ basis set for valence
electrons with a single polarization function (TZVP), and the a posteriori D3
Grimme correction to account for dispersion forces (Grimme et al., 2010,
2011). The plane-wave cutoff was set to 600 Ry. Ice surfaces were thermalized
(see sec. Ice models), and during the geometry optimization only the reacting
H atoms (those forming H2) were allowed to move, keeping the water molecules
fixed at their thermalized positions. All calculations were carried out within
the unrestricted formalism as we deal with open-shell systems (see SI for more
details). The binding energy (BE) of H2 was calculated as:

$BE_{H_{2}}=E_{CPLX}-(E_{Ice}+E_{H_{2}})$ (1)

where $E_{CPLX}$ is the energy of the H2/Ice system, $E_{Ice}$ that of the
bare ice surface, and $E_{H_{2}}$ the energy of the H2 alone.

AIMDs were carried out within the NVE (N = number of particles, V = volume, E
= energy) ensemble, where the total energy VTOT (i.e. potential + kinetic) is
conserved. The evolution of the system was followed for 5 ps for the
crystalline model and for 1 ps for the amorphous model, using time-steps of
0.2 fs.

###  2.2 Ice models

For the crystalline ice model, a periodic slab cut from the hexagonal ice bulk
structure was used. The periodic cell parameters defining the system are:
$\textit{a}=26.318$ Å and $\textit{b}=28.330$ Å, while the slab thickness is
$\sim 21$ Å (corresponding to seven layers), for a total of 576 water
molecules. The c parameter of the simulation box, i.e. the non-periodic one,
was set to 50 Å to avoid interactions among the fictitious slab replicas. The
size of the slab model was chosen to avoid non-physical temperature increase
due to the extremely large exothermic reaction (Pantaleone et al., 2020).

The amorphous model was obtained by performing a classical molecular dynamics
(MD) simulation on the crystalline structure. The MD was carried out with the
TIP3P force-field (Jorgensen et al., 1983) for 200 ps (with a timestep of 0.5
fs) at 300 K within the NVT (N = number of particles, V = volume, T =
Temperature) ensemble. A second NVT simulation was performed at 10 K, to cool
down the system at the temperature of cold molecular clouds. Finally, a
geometry optimization and another NVT-MD simulation at 10 K were run using
PBE, to recover the potential energy surface at the same theory level as
described in the previous section.

###  2.3 H2 vibrational state

To calculate the H2 vibrational state during the AIMDs the anharmonic
oscillator model was employed. As a first step, the PES (potential energy
surface) of the H2 isolated molecule was explored by performing a rigid scan
of the H–H distance, ranging from 0.3 to 3.0 Å with a step of 0.01 Å. Computed
data were fitted with the Morse equation to obtain the force constant of the
oscillator, its dissociation energy and, hence, the vibrational levels of H2.
The H2 turning points calculated with our model were compared with those of
the AIMD simulations and the vibrational level was assigned at each H2
oscillation during the AIMDs.

##  3 Results

###  3.1 Reactants plus product positions in the ice models

Figure 1: Top view of the crystalline (a), and top (b) and side (c and d) views of the ASW ice models. The yellow circles represent the centers of mass of the two H atoms starting positions. In the case of ASW, the numbers mark the positions (Pos1, Pos2 and Pos3) discussed in the text. The initial H to H distance is about 4 Å (see Supplementary Fig. LABEL:h2_ene). O atoms are in red and H atoms of the ice are in white. Figure 2: PBE-D3/TZVP optimized geometries of reactants (left panels) and product (right panels) of the H2 reaction formation on the crystalline (top) and ASW Pos1 (bottom) ice models. The numbers in parenthesis correspond to the relative energy in kJ/mol with respect to the reactants. The H atoms involved in the H + H $\longrightarrow$ H2 reaction are in yellow, those belonging to surface water molecules are in white, and O atoms are in red. Table 1: Results of the AIMD for the H2 formation on the crystalline (first row) and the ASW three positions (Fig. 1) (bottom three rows) ices. Second and third columns report the fraction of kinetic energy of the ice and the newly formed H2, respectively, at the end of the simulation (5 and 1 ps for crystalline and ASW ices, respectively). Fourth column reports the H2 binding energy (in kJ/mol) at the reaction site. Fifth and sixth columns report the vibrational state averaged over 0.2–1 ps and 1–5 ps, respectively (the latter can be computed only for the crystalline case). Last column reports the peak energy (in kJ/mol) of an ice-molecule neighbor (i.e. at $\sim 3$ Å) to where the reaction occurs.  Model | $\frac{K_{ice}}{K_{ice}+K_{H_{2}}}$ | $\frac{K_{H_{2}}}{K_{ice}+K_{H_{2}}}$ | BE${}_{H_{2}}$ | $\nu^{0.2--1}$ | $\nu^{1--5}$ | Epeak
---|---|---|---|---|---|---
|  |  | [kJ/mol] |  |  | [kJ/mol]
Crys. | 0.45 | 0.55 | 9.0 | 4–5 | 1–2 | 12.0
Pos1 | 0.65 | 0.35 | 9.7 | 1–2 | – | 43.0
Pos2 | 0.46 | 0.54 | 10.6 | 6 | – | 7.0
Pos3 | 0.47 | 0.53 | 12.4 | 5 | – | 7.0

Figure 1 shows the structure of the used crystalline and ASW ice models as
well as the starting positions of the center of mass of the two adsorbed H
atoms before the reaction.The starting H–H distance, as well as its evolution
during the AIMDs is shown in Supplementary Figure LABEL:h2_ene. The ASW model
presents a rugged morphology (e.g. holes, cavities, channels) in comparison to
the crystalline one and, therefore, we carried out simulations in three
positions, roughly representing different possible situations in terms of
energetics and surface morphology. Position marked as Pos1 is the deepest one,
within a cavity, while positions Pos2 and Pos3 are in the outermost parts of
the surface.

As a first step, we optimized the geometries of reactants (the two H atoms)
and product (the H2 molecule) on the water ice surface in order to obtain the
potential energy surface of the reaction. No search for transition state
structures has been performed because AIMDs indicated this reaction is
barrierless at 10 K. That is, the starting kinetic energy of the two H atoms
provided by the 10 K is high enough to overcome the eventually small potential
energy barrier. The total energy to be dissipated is around 435-440 kJ/mol,
depending on the surface and starting position.

Figure 2 provides a view of the two H atoms starting positions and H2 position
immediately after the reaction, on the crystalline and ASW (Pos1 as an
example) models, respectively.

###  3.2 Molecular hydrogen desorption

AIMDs results are summarised in Table 1 and shown in Fig. 3. First, at least
half of the kinetic energy is absorbed by the ice while the remaining one is
kept by the newly formed H2 molecule ($\sim$50–35%). These numbers are similar
for both for the crystalline and ASW ice positions, suggesting they do no
depend much on the surface structural details of the ice.

Second, after the H–H bond formation, the newly formed H2 molecule diffuses
over the surface, as it keeps a large kinetic energy, in both crystalline and
ASW ice models. However, despite the energetics of the two processes are
similar (Figs. 3 left panels) the H2 diffusion (Figs. 3 right panels) is
different.

On the crystalline surface, the diffusion of the newly formed H2 over the ice
surface is constrained in a specific direction, along the a-axis, parallel to
the ice surface, whereas in the b-axis and c-axis directions the starting
position does not change. This is because the crystalline model has alternated
and opposite electrostatic potentials along the b-axis direction (see
Supplementary Fig. 2, which constrains the H2 diffusion to the perpendicular
a-axis, within a channel-like structure created in between the two alternated
potential regions (see Figs. S2a and S2c). In contrast, on the ASW ice model,
H2 diffuses over all the three directions, also along the c-axis which
corresponds to direction perpendicular to the ice surface. However, the
movement is not the same in the three studied positions. In Pos1, H2 achieves
a maximum height of $\sim$10 Å over the surface in a timescale of 0.4 ps and
does not come back (Fig. 3-d), this way leaving definitely the surface. In
Pos2, on the contrary, H2 slides and hops on the surface within the 1 ps of
the simulation (Fig. 3-f). A similar behavior is observed for Pos3, but with
wider hops and jumps (Fig. 3-h). However, and this is the key point of the
simulations, the kinetic energy remains much larger ($\sim~{}150$ kJ/mol) than
the binding energy ($\sim 10$ kJ/mol) in all cases: therefore, the nascent H2
is likely doomed to leave the ice surface.

###  3.3 Energy dissipation towards the ice water molecules

Fig. 4 shows how the energy absorbed by the ice is first transmitted from
reaction site to a neighbor water(i.e. at $\sim$4 Å) in $\sim$200–800 fs and,
with time, to shells of water molecules with increasing distances. The
neighbor water molecule acquires from $\sim$3 to $\sim$14 kJ/mol and stays
with that energy for more than about 100–200 fs before the energy is
dissipated towards more external water molecules. These timescales are in
excellent agreement with previous studies of the sound speed in ices (see e.g.
Ruocco et al., 1996).

###  3.4 H2 vibrational state

The initial state of the H2 formed on the crystalline ice is $\nu$ = 4–5
(obtained averaging over 0.2–1 ps) and then it decreases to $\nu$ = 1–2 (1–5
ps average) because of the H2 energy dissipation. Conversely, at Pos1 of the
ASW model, H2 is formed with a low vibrational state ($\nu$ = 1–2, 0.2–1 ps
average).

Figure 3: Results of the AIMDs for the crystalline (top panels) and ASW Pos1,
Pos2, Pos3 (panels on the second, third and bottom rows) models, respectively.
Left panels: Evolution over time (in fs) of the most relevant energetic
components (in kJ/mol) of the H2/ice system for the crystalline (panel a) and
ASW Pos1, Pos2 and Pos3 (panel c, e and g) models, respectively. VTOT is the
total potential energy (green lines), K${}_{H_{2}}$ and KIce are the kinetic
energies of H2 (red lines) and ice (blue lines), respectively, while
BE${}_{H_{2}}$ (10 kJ/mol = 1200 K) is the binding energy of the H2 (black
lines). Right panels: Diffusion of the center of mass of the H2 molecule split
into the three cartesian components: c-axis is the direction perpendicular to
the ice surface, while a-axis and b-axis are along the ice surface (§ 2).

For the other two positions (Pos2 and Pos3), the vibrational state of the
newly formed H2 molecule lies in highly excited states, $\nu$ = 6 and $\nu$ =
5, respectively. However, our simulations stop at 1 ps and, likely, as on the
crystalline ice, the vibrational state of the newly formed H2 will lower to
$\nu$ = 1–2. The difference in the behavior between the Pos1 and Pos2/3 is
probably due to a transient chemical bond between one of the reacting H and
one O atom of a neighboring water molecule (see Supplementary Figure 3.

##  4 Discussion

###  4.1 H2 and other molecules desorption

While there were no doubts that, once formed on the icy surfaces, H2 molecules
will be released into the gas-phase, our new computations provide a
quantitative and atomistic-view support to this theory. The newly formed H2
molecules possess an energy much larger than their binding energy, by more
than a factor ten, so that very likely H2 will be injected into the gas-phase.
Even in the worst case of crystalline ice, where the simulations show H2
sliding over the perfect ice surface, the newly formed H2 will sooner or later
stumble on an imperfection of the ice so that its trajectory would be deviated
and the molecule escape into the gas. To have an estimation of the timescale
of this phenomenon one could consider the timescale for the newly H2 to scan
the entire surface of an interstellar grain. Assuming its radius equal to 0.1
$\mu$m and each icy site size equal to 3 Å, the number of sites to scan are
about 105. A first estimate of the time to scan all the sites can be directly
obtained by our simulations by considering that about 15 sites of our
crystalline ice are covered in less 1 ps: therefore, to scan the 105 it will
take less than 10 ns. However, this is strictly true if the H2 molecule does
not loose energy in other minor impacts, so $\sim 10$ ns can be considered a
lower limit. The upper limit can be computed by assuming that no residual
energy is left to the H2 molecule except the thermal one and, in this case,
the velocity to scan is given by the hopping rate and the timescale for
scanning the entire grain surface becomes:

$t_{scan}=N_{s}\nu_{0}^{-1}e^{\frac{E_{d}}{k_{B}T_{dust}}}$ (2)

where $\nu_{0}$ is about 1012 s-1 and Ed is assumed to be 0.3 times the
binding energy (about 400 K). Inserting the numbers, a rough estimate of the
time that H2 takes to leave a typical interstellar grain is $\leq 1000$ yr. In
laboratory analogues of interstellar grains, the timescale would be even
larger, and, consequently, not observable. In conclusion, the newly formed H2
will leave the surface in a timescale between a few ns to a max of about 1000
yr, which is still a very short lifetime with respect to the cloud life.

We want to highlight that the choice of using a proton order ice is just a
matter of convenience to test our simulations, before going to the more
realistic case of the amorphous surface. On a more realistic proton disordered
crystalline ice we expect a more anisotropic behavior of the H2 molecule, and,
as a consequence, a faster H2 desorption.

Our computations also show that the ice absorbs a significant fraction of the
reaction energy, 45–65% (Tab. 1). This energy propagates like a wave from the
point where the H + H reaction occurs (see Fig. 4 and also Pantaleone et al.,
2020). In ASW ice, a water molecule close to reaction site acquires 7–43
kJ/mol, (840–5160 K), for more than 100-200 fs. Within the first shell of
radius 4 Å, water molecules acquire energies from 3 to 14 kJ/mol (360 to 1680
K, average value counting all the water molecules within the first shell). and
farther away the energy acquired by the ice-water molecules decreases to less
than 1.6 kJ/mol. The energy acquired by the water ice molecules within a
radius of 4 Å from the reaction site could, potentially, be enough to release
into the gas-phase any molecule whose binding energy is lower than 3–14
kJ/mol. This could be the case of another H2 molecule frozen on the ice, a
possibility predicted by some astrochemical models (see e.g. Hincelin et al.,
2015). Indeed, since the binding energy of H2 is less than $\sim$5 kJ/mol
(e.g. Vidali, 2013; Ferrero et al., 2020; Molpeceres & Kästner, 2020), several
H2 molecules could be kicked into the gas-phase.

Figure 4: Kinetic energy (in kJ/mol) acquired by the water ice molecules as a
function of time for the crystalline (top first row panels) and ASW Pos1-3
(lower panels) ices, as marked by the left-hand labels. Left first column
panels: Energy acquired by a single neighboring water ice molecule (i.e. at
$\sim 3$ Å) from the H2 reaction site. Other column panels: Energy acquired by
the ice surface divided in concentric shells (normalized by the number of
water molecules per shell), centered at the reaction site. The labels on the
top mark the shells radii.

More interesting is the case of CO, the most abundant molecule after H2 in
galactic cold molecular clouds. It has a binding energy between 7 and 16
kJ/mol (e.g. He et al., 2016; Ferrero et al., 2020) so that frozen CO could
potentially be liberated into the gas-phase by the formation of H2. It has
been long recognized that, in absence of a non-thermal desorption process, CO
in molecular clouds should be entirely frozen onto the interstellar ices in
$\sim 3\times 10^{5}$ yr (see e.g. Léger, 1983). Several mechanisms have been
proposed to explain its presence in the gas (see e.g. Leger et al., 1985;
Schutte & Greenberg, 1991; Duley & Williams, 1993; Shen et al., 2004; Ivlev et
al., 2015). They are all based on the idea that the grain is locally heated
and that the frozen molecules are statistically desorbed on a time scale of a
few thousands seconds (see e.g. Hasegawa & Herbst, 1993; Herbst & Cuppen,
2006), except when faster photodesorption processes are involved. In
particular, Duley & Williams (Duley & Williams, 1993) and Takahashi & Williams
(Takahashi & Williams, 2000) focused on the CO desorption caused by the H2
formation. Using the results of classical MD simulations by Takahashi and
colleagues (Takahashi et al., 1999b; Takahashi, 1999; Takahashi & Williams,
2000) and a statistical approach, these authors found that the surface within
a radius of $\sim$4 Å from the reaction site is heated up to $\leq$40 K for a
too short time ($\sim 10$ fs) to allow CO to sublimate. However, despite the
Takahashi’s simulations are impressive considering that they were carried out
more than 20 years ago, old force fields (like the TIP3P used by Takahashi et
al.) have limitations in describing the actual intermolecular interactions and
the energy transfer from one molecule to another, compared to our simulations
treated at quantum mechanical level.

Our AIMDs show that the water molecules within a radius of 4 Å from the
reaction site can be excited from 3 to 14 kJ/mol (depending on the position)
for more than 100 fs (Fig. 4). Using the range of values in the literature for
the CO binding energies (7–16 kJ/mol: He et al., 2016; Ferrero et al., 2020),
and the usual Eyring equation to estimate the half-life time of CO desorption,
we obtain that a CO molecule within a radius of 4 Å from Pos1 would desorb in
33–62 fs, and 69–185 fs from Pos2. On the contrary, CO molecules close to Pos3
and on crystalline ice would not desorb. Therefore, based on this rough
argument, the formation of H2 can potentially desorb frozen CO molecules.
Whether this hypothesis is realistic or not, it depends on (i) how many sites
of the AWS ice water molecules are excited as in Pos1 and Pos2, (ii) the ratio
of H2 formation rate with respect to the CO freezing rate
($R=\frac{v(H_{2})}{v(CO)}$), and (iii) the probability that CO and H2 are
adjacent.

While a more realistic model for the ASW ice is needed to estimate the first
quantity (i), one can estimate the second one (ii) as follows. At steady
state, assuming that the H + H reaction has efficiency 1 (e.g. Vidali, 2013),
the ratio $R$ of the H2 formation rate with respect to the CO freezing rate is
equal to the ratio between the rate of H atoms over CO molecules landing on
the grain surface, divided by 2 (because two H atoms are needed for the H2
formation). Considering an average Milky Way molecular cloud with H2 density
of $10^{3}-10^{4}$ cm-3, temperature 10 K, a gaseous (undepleted) CO abundance
with respect to H2 equal to $2\times 10^{-4}$, and a cosmic ray ionisation
rate of $3\times 10^{-17}$ s-1, one obtains $\rm{n_{H}}\sim 2$ cm-3 and $R\sim
26-2.6$. That is, the H2 formation rate dominates over the CO freezing one.
Therefore, frozen CO molecules can potentially be desorbed by the energy
released by the H2 formation on the icy grain surfaces. This process, if
confirmed, might naturally explain the presence of gaseous CO in not too dense
(n${}_{H}\lesssim 10^{4}$ cm-3) molecular clouds and solve a decades long
mystery. However, to put this affirmation on a solid ground, specific AIMDs
showing the effective sublimation of the CO molecule as well as larger ASW ice
models and dedicated astrochemical modeling simulations that include this
effect are mandatory. They will be the focus of forthcoming works.

###  4.2 H2 vibrational state

Previous experimental and computational works on graphite surfaces show that
after its formation, the H2 molecule populates vibrational levels around $\nu$
= 3–4 (Latimer et al. (2008); Islam et al. (2010); Casolo et al. (2013)). On
crystalline and amorphous ice surfaces only a few studies were carried out.
Based on classical MD simulations, Takahashi et al. (1999a); Takahashi &
Uehara (2001) predicted that H2 formed on amorphous ice would be vibrationally
excited to $\nu$ = 7–8. In contrast, in some laboratory experiments Roser et
al. (2003); Hornekær et al. (2003); Hama et al. (2012), the authors did not
detect such excited states. Our new simulations confirm both findings:
depending on the site where the formation occurs, H2 can have large (up to 6)
or low (1–2) $\nu$ (Tab. 1). For example, the case of ASW Pos1, which is the
one in a cavity, shows the lowest $\nu$. In the ISM, the vibrational state of
the newly formed H2 molecules will depend on the probability of H2 leaving the
grain surface before being ”thermalised” by collisions with ice water
molecules, as also found and discussed in Watanabe et al. (2010).

The vibrational state of the newly formed H2 molecule, i.e. how much the
molecule is vibrationally excited, has important consequences on two major
astrophysical aspects: (i) to help gas-phase reactions with high activation
energy barriers, some of which are the starting points of the entire chemistry
in molecular clouds Agúndez et al. (2010); (ii) the possibility to detect
newly formed H2 with near-future JWST (James Webb Space Telescope)
observations and, consequently, measure the rate of H2 formation in molecular
clouds and put constraints to theories.

## Summary

We studied the energy dissipation of the $H_{2}$ formation reaction on both
crystalline and amorphous (ASW) ice models by means of state-of-art ab-initio
molecular dynamics simulations (AIMDs). In the ASW ice, we explored three
formation sites meant to represent different situations in terms of energetics
and surface morphology: one is in a cavity (Pos1) and two in the outermost
parts of the surface (Pos2 and Pos3).

In all simulations, we found that about at least 30% of the reaction energy
($\sim$440 kJ/mol) is acquired by the newly formed H2, namely more than ten
times the H2 binding energy ($\sim$10 kJ/mol), so that H2 is likely doomed to
leave the ice and to be injected into the gas. The remaining two thirds of the
reaction energy are absorbed by the water ice in less than 1 ps. The water
molecules nearby to the reaction site have energy peaks of 7–43 kJ/mol for
more than 100–200 fs, while those within a 4.0 Å radius of 3–14 kJ/mol.

We showed that it is in principle possible that frozen CO molecules close to
the H2 formation site sublimate. If confirmed, this will be a simple
explanation to the decades-long conundrum of why gaseous CO is present in cold
molecular clouds. In order to quantify this effect new focused AIMDs adopting
larger ASW ice models and dedicated astrochemical modeling will be necessary.

Finally, the nascent H2 molecules have large ($\nu=6$) vibrational states in
the first ps and later decay to 1–2. This high vibrational state could help
reactions with an activation barrier involving H2 to occur also in cold gas
and be observable by JWST.

The authors acknowledge funding from European Union’s Horizon 2020 research
and innovation program, the European Research Council (ERC) Project “the Dawn
of Organic Chemistry”, grant agreement No 741002, the European Research
Council (ERC) Project ”Quantum Chemistry on Interstellar Grains”, grant
agreement No 865657, and the Marie Skłodowska-Curie project “Astro-Chemical
Origins” (ACO), grant agreement No 811312. AR is indebted to the “Ramón y
Cajal” program. MINECO (project CTQ2017-89132-P) and DIUE (project
2017SGR1323) are acknowledged. BSC-MN and OCCIGEN HPCs are kindly acknowledged
for the generous allowance of super-computing time through the QS-2019-2-0028
and 2019-A0060810797 projects, respectively. SP, NB and PU acknowledge the
Italian Space Agency for co-funding the Life in Space Project (ASI N.
2019-3-U.O

## References

  * Agúndez et al. (2010) Agúndez, M., Goicoechea, J., Cernicharo, J., Faure, A., & Roueff, E. 2010, The Astrophysical Journal, 713, 662
  * Boogert et al. (2015) Boogert, A. C. A., Gerakines, P. A., & Whittet, D. C. B. 2015, Annual Review of Astronomy and Astrophysics, 53, 541, doi: 10.1146/annurev-astro-082214-122348
  * Casolo et al. (2013) Casolo, S., Tantardini, G. F., & Martinazzo, R. 2013, Proceedings of the National Academy of Sciences, 110, 6674
  * Duley & Williams (1993) Duley, W. W., & Williams, D. A. 1993, Mon. Not. R. Astron. Soc., 260, 37, doi: 10.1093/mnras/260.1.37
  * Ferrero et al. (2020) Ferrero, S., Zamirri, L., Ceccarelli, C., et al. 2020, The Astrophysical Journal, 904, 11
  * Goedecker et al. (1996) Goedecker, S., Teter, M., & Hutter, J. 1996, Phys. Rev. B, 54 (3), 1703, doi: 10.1103/PhysRevB.54.1703
  * Grimme et al. (2010) Grimme, S., Antony, J., Ehrlich, S., & Krieg, H. 2010, J. Chem. Phys., 132 (15), 154104, doi: 10.1063/1.3382344
  * Grimme et al. (2011) Grimme, S., Ehrlich, S., & Goerigk, L. 2011, J. Comput. Chem., 32, 1456, doi: 10.1002/jcc.21759
  * Hama et al. (2012) Hama, T., Kuwahata, K., Watanabe, N., et al. 2012, The Astrophysical Journal, 757, 185
  * Hasegawa & Herbst (1993) Hasegawa, T. I., & Herbst, E. 1993, Mon. Not. R. Astron. Soc., 261, 83, doi: 10.1093/mnras/261.1.83
  * He et al. (2016) He, J., Acharyya, K., & Vidali, G. 2016, The Astrophysical Journal, 825, 89, doi: 10.3847/0004-637X/825/2/89
  * Herbst & Cuppen (2006) Herbst, E., & Cuppen, H. M. 2006, Proceedings of the National Academy of Sciences, 103, 12257
  * Hincelin et al. (2015) Hincelin, U., Chang, Q., & Herbst, E. 2015, Astronomy & Astrophysics, 574, A24
  * Hollenbach & Salpeter (1970) Hollenbach, D., & Salpeter, E. 1970, The Journal of Chemical Physics, 53, 79
  * Hollenbach & Salpeter (1971) —. 1971, The Astrophysical Journal, 163, 155
  * Hornekær et al. (2003) Hornekær, L., Baurichter, A., Petrunin, V., Field, D., & Luntz, A. 2003, Science, 302, 1943
  * Hutter et al. (2014) Hutter, J., Iannuzzi, M., Schiffmann, F., & VandeVondele, J. 2014, WIREs Comput. Mol. Sci., 4, 15, doi: 10.1002/wcms.1159
  * Islam et al. (2010) Islam, F., Cecchi-Pestellini, C., Viti, S., & Casu, S. 2010, The Astrophysical Journal, 725, 1111
  * Ivlev et al. (2015) Ivlev, A., Röcker, T., Vasyunin, A., & Caselli, P. 2015, The Astrophysical Journal, 805, 59
  * Jorgensen et al. (1983) Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W., & Klein, M. L. 1983, J. Chem. Phys., 79, 926
  * Latimer et al. (2008) Latimer, E. R., Islam, F., & Price, S. D. 2008, Chemical Physics Letters, 455, 174
  * Léger (1983) Léger, A. 1983, Astronomy and Astrophysics, 123, 271
  * Leger et al. (1985) Leger, A., Jura, M., & Omont, A. 1985, Astronomy and Astrophysics, 144, 147
  * Lippert et al. (1997) Lippert, G., Hutter, J., & Parrinello, M. 1997, Mol. Phys., 92 (3), 477, doi: 10.1080/002689797170220
  * Molpeceres & Kästner (2020) Molpeceres, G., & Kästner, J. 2020, Physical Chemistry Chemical Physics (Incorporating Faraday Transactions), 22, 7552, doi: 10.1039/D0CP00250J
  * Pantaleone et al. (2020) Pantaleone, S., Enrique-Romero, J., Ceccarelli, C., et al. 2020, The Astrophysical Journal, 897, 56
  * Perdew et al. (1996) Perdew, J., Burke, K., & Ernzerhof, M. 1996, Phys. Rev. Lett., 77 (18), 3865, doi: 10.1103/PhysRevLett.77.3865
  * Roser et al. (2003) Roser, J., Swords, S., Vidali, G., Manico, G., & Pirronello, V. 2003, The Astrophysical Journal Letters, 596, L55
  * Ruocco et al. (1996) Ruocco, G., Sette, F., Bergmann, U., et al. 1996, Nature, 379, 521, doi: 10.1038/379521a0
  * Schutte & Greenberg (1991) Schutte, W., & Greenberg, J. 1991, Astronomy and Astrophysics, 244, 190
  * Shen et al. (2004) Shen, C., Greenberg, J., Schutte, W., & Van Dishoeck, E. 2004, Astronomy & Astrophysics, 415, 203
  * Smith et al. (1988) Smith, R. G., Sellgren, K., & Tokunaga, A. T. 1988, The Astrophysical Journal, 334, 209
  * Smoluchowski (1983) Smoluchowski, R. 1983, J. Phys. Chem., 87, 4229
  * Takahashi (1999) Takahashi, J. 1999, Earth, planets and space, 51, 1215
  * Takahashi et al. (1999a) Takahashi, J., Masuda, K., & Nagaoka, M. 1999a, The Astrophysical Journal, 520, 724
  * Takahashi et al. (1999b) —. 1999b, Monthly Notices of the Royal Astronomical Society, 306, 22
  * Takahashi & Uehara (2001) Takahashi, J., & Uehara, H. 2001, The Astrophysical Journal, 561, 843
  * Takahashi & Williams (2000) Takahashi, J., & Williams, D. A. 2000, Monthly Notices of the Royal Astronomical Society, 314, 273
  * van de Hulst (1946) van de Hulst, H. C. 1946, Recherches Astronomiques de l’observatoire d’Utrecht, 11, 2
  * Vidali (2013) Vidali, G. 2013, Chemical Reviews, 113, 8762, doi: 10.1021/cr400156b
  * Wakelam et al. (2017) Wakelam, V., Bron, E., Cazaux, S., et al. 2017, Mol. Astrophys., 9, 1
  * Watanabe et al. (2010) Watanabe, N., Kimura, Y., Kouchi, A., et al. 2010, The Astrophysical Journal Letters, 714, L233
  * Whittet (1993) Whittet, D. C. B. 1993, Observations of Molecular Ices, ed. T. J. Millar & D. A. Williams, 9

##  Appendix A Computational details

Proper reactants spin localisation is achieved by relying on the unrestricted-
DFT broken (spin) symmetry approach. This is mandatory when dealing with
diradical systems in singlet electronic state. In the H + H specific case,
despite the total spin moment on the global system is 0, the unpaired
electrons on the two H atoms are expected to have spin $\alpha$ (i.e. +1) and
$\beta$ (i.e. -1) each.

In order to check the quality of our results, in Table LABEL:h2_ene a
comparison between the H2 formation energy in the gas phase calculated at PBE
and CCSD(T) levels is presented Moreover, also the interaction between H2 and
one H2O molecule was evaluated at the same levels of theory, giving results in
good agreement. Moreover, also the interaction between H2 and one H2O molecule
was evaluated at the same levels of theory, giving results in good agreement.
The difference between the two methodologies is acceptable, considering the
energy released by the reaction.

As regards the AIMD simulations, we run an equilibration AIMD in the NVT
ensemble (using the CSVR thermostat, with a time constant of 20 femtoseconds)
at 10 K for 1 ps (with a time step of 1 fs) only for the bare ice surface.
This ensures an initial thermally equilibrated ice. Then the velocities of the
equilibrated ice surfaces are used as input for the NVE production runs. Those
velocities, correspond to the initial kinetic energy of the ice, as presented
in Figure 3. The H velocities, instead, were manually set to favor the H–H
bond formation respecting the 10 K limitation.

List of supptables 1 Adsorption and reaction (H + H $\longrightarrow$ H2) energy data. | Reaction energy (kJ mol-1)
---|---
Crystalline | -438.9
Pos1 | -436.2
Pos2 | -435.9
Pos3 | -440.3
H2 CP2Ka | -440.6
H2 Gaussianb | -457.4
| Binding energy (kJ mol-1)
CP2Ka | 4.1
Gaussianb | 2.6
a H2 gas phase reaction energy calculated at PBE-D3/TZVP level
b H2 gas phase reaction energy calculated at CCSD(T)/aug-cc-pv5z level
List of suppfigures 1 Evolution of the H–H bond distance during the AIMD
simulations. List of suppfigures 2 Top view of the electrostatic potential map
of the crystalline (a) and amorphous (b) ice models. Lateral view of the
crystalline ice model (c). White arrows represent the H2 diffusion direction.
The light blue circle in the bottom panel represents the cavity inside of
which the H2 diffuses over the crystalline ice model. Red and blue zones of
the electrostatic maps correspond to positive and negative potentials. O atoms
in red, H atoms in white. List of suppfigures 3 Evolution of the kinetic
energy of H2 (panel (a)) and a neighbor water molecule (panel (b)) for the
amorphous Pos1 case. Starting position of the AIMD simulation (panel (c)) and
formation of the H-/H3O+ adduct after 0.2 picoseconds of simulation (panel
(d)). O atoms in red, H atoms in white. The H atoms highlighted in green are
the reactants.