Department of Computer Science and Engineering, NYU Tandon, New York, USA and
https://cse.engineering.nyu.edu/~hstein/
lisa.hellerstein@nyu.eduhttps://orcid.org/ 0000-0002-3743-7965Management
Science and Information Systems Department, Rutgers Business School, New
Jersey,
USAtlidbetter@business.rutgers.eduhttps://orcid.org/0000-0001-6111-2899
Department of Computer Science and Engineering, NYU Tandon, New York, USA and
https://www.rtealwitter.com
rtealwitter@nyu.eduhttps://orcid.org/0000-0003-3096-3767 John Q. Public and
Joan R. Public <ccs2012> <concept>
<concept_id>10003752.10003809.10003636.10003812</concept_id>
<concept_desc>Theory of computation Facility location and
clustering</concept_desc> <concept_significance>300</concept_significance>
</concept> </ccs2012> [300]Theory of computation Facility location and
clustering

###### Acknowledgements.

# A Local Search Algorithm for the Min-Sum Submodular Cover Problem

Lisa Hellerstein  Thomas Lidbetter  R. Teal Witter

###### Abstract

We consider the problem of solving the Min-Sum Submodular Cover problem using
local search. The Min-Sum Submodular Cover problem generalizes the NP-complete
Min-Sum Set Cover problem, replacing the input set cover instance with a
monotone submodular set function. A simple greedy algorithm achieves an
approximation factor of 4, which is tight unless P=NP [Streeter and Golovin,
NeurIPS, 2008]. We complement the greedy algorithm with analysis of a local
search algorithm. Building on work of Munagala et al. [ICDT, 2005], we show
that, using simple initialization, a straightforward local search algorithm
achieves a $(4+\epsilon)$-approximate solution in time
$O(n^{3}\log(n/\epsilon))$, provided that the monotone submodular set function
is also second-order supermodular. Second-order supermodularity has been shown
to hold for a number of submodular functions of practical interest, including
functions associated with set cover, matching, and facility location. We
present experiments on two special cases of Min-Sum Submodular Cover and find
that the local search algorithm can outperform the greedy algorithm on small
data sets.

###### keywords:

Local search, submodularity, second-order supermodularity, min-sum set cover

###### category:

##  1 Introduction

We consider the Min-Sum Submodular Cover problem, defined as follows. The
input to the problem consists of an oracle for a monotone submodular function
$u:2^{[n]}\rightarrow\operatorname*{\mathbb{R}}_{\geq 0}$, and positive costs
$c_{1},\ldots,c_{n}\in\operatorname*{\mathbb{R}}_{>0}$, where
$[n]=\\{1,\ldots,n\\}$. Let
$c:2^{[n]}\rightarrow\operatorname*{\mathbb{R}}_{\geq 0}$ be such that for all
$S\subseteq[n]$, $c(S)=\sum_{i\in S}c_{i}$. We refer to $u$ as the “utility”
and $c$ as the “cost” function. The problem is to find the permutation of the
elements of $[n]$ that minimizes

$\displaystyle\sum_{i=1}^{n}c(S_{i})\left(u(S_{i})-u(S_{i-1})\right)$ (1)

where $S_{i}$ is the set containing the first $i$ elements of the permutation.
The Min-Sum Submodular Cover problem generalizes the NP-Complete Min-Sum Set
Cover problem introduced by Feige et al. [5]. It has a simple greedy algorithm
that achieves a 4-approximation [18, 9].111See the Appendix A.1 for comments
on Streeter and Golovin [18]. The 4-approximation is tight, assuming $P\neq
NP$ [5].

In this work, we analyze a local search algorithm for Min-Sum Submodular
Cover. Local search algorithms have been extensively applied to discrete
optimization problems [3, 14, 15] and offer several benefits over other types
of algorithms [2]. One advantage of local search algorithms compared to greedy
methods in practice is their ability to explore a diverse set of solutions. In
Section 4, we present the results of preliminary experiments which demonstrate
that this ability can yield improved solutions to Min-Sum Submodular Cover
instances.

The local search algorithm we consider works in iterations, starting with an
initial solution. In an iteration, the algorithm updates the current solution
to its best “neighbor”. Once a solution is locally optimal (or after a fixed
number of iterations), the algorithm returns the current solution.

Our analysis builds on previous work of Munagala et al. [16] for the Pipelined
Set Cover problem (defined below). Munagala et al. exploit the observation
that utility can be attributed to elements of the ground set for the set cover
instance, and use these elements as variables in their linear program. The
main challenge of generalizing their analysis is that the utility in general
submodular set functions is more abstract and cannot be attributed to
particular objects. In order to apply the linear program used in Munagala et
al., our analysis relies on an additional property of the utility function
called second-order supermodularity.

Second-order supermodularity was first studied by Korula et al. [12]. It can
be viewed as a natural extension of submodularity: If one considers the
multilinear extension $F:[0,1]^{n}\rightarrow\operatorname*{\mathbb{R}}_{\geq
0}$ of a set function
$f:\\{0,1\\}^{n}\rightarrow\operatorname*{\mathbb{R}}_{\geq 0}$ (which is a
way of interpolating the values of $f$ from the vertices of the Boolean
hypercube to points in its interior), the submodularity of $f$ is equivalent
to the property that the second partial derivatives of $F$ are non-positive.
As mentioned by Korula et al. [12] and Iyer et al. [11], the second-order
supermodularity of $f$ is equivalent to the property that the third partial
derivatives of $F$ are non-negative.

The second-order supermodularity property is not overly restrictive; there are
several classes of submodular functions that have this property including
weighted coverage functions, weighted matching functions, and facility
location [12, 10, 11]. Since this property was first defined, improved bounds
have been obtained for optimization problems by assuming the property [12,
10]. The related properties of second-order modularity and second-order
submodularity have also been used in analyzing local search algorithms for
constrained submodular maximization problems [6, 7].

We now define two special cases of Min-Sum Submodular Cover: the Pipelined Set
Cover problem and the Min-Sum Facility Location problem.

#### Pipelined Set Cover

The inputs to the problem consist of (i) $m$ “ground” elements
$\\{1,\ldots,m\\}=[m]$, (ii) $D_{1},\ldots,D_{n}$, a family of $n$ subsets of
the ground elements $[m]$ such that $\bigcup_{i\in[n]}D_{i}=[m]$, and (iii)
positive costs $c_{1},\ldots,c_{n}$ associated with each $D_{i}$. Let
$u:2^{[n]}\rightarrow\operatorname*{\mathbb{R}}_{\geq 0}$ be such that for all
$S\subseteq[n]$, $u(S)$ is the number of ground elements in $\bigcup_{i\in
S}D_{i}$. We call $u$ a “coverage” function. Let
$c:2^{[n]}\rightarrow\operatorname*{\mathbb{R}}_{\geq 0}$ be such that for all
$S\subseteq[n]$, $c(S)=\sum_{i\in S}c_{i}$. The problem is to find the
permutation of $[n]$ that minimizes the objective function in the Min-Sum
Submodular Cover problem,
$\sum_{i=1}^{n}c(S_{i})\left(u(S_{i})-u(S_{i-1})\right)$. Thus Pipelined Set
Cover problem is equivalent to the special case of the Min-Sum Submodular
Cover problem where the utility function $u$ is a coverage function. The Min-
Sum Set Cover problem is the special case of Pipelined Set Cover with unit
costs. (We note that Munagala et al. also present results for Weighted
Pipelined Set Cover, where the ground elements have weights.)

#### Min-Sum Facility Location

Consider the following problem facility location problem, studied by Krause
and Golovin [13]. There is a set $[n]$ of possible locations where facilities
could be opened, to serve a collection of $m$ customers. Opening a facility at
location $a$ provides a service of value $M_{a,b}$ to customer $b$, where
$M\in\operatorname*{\mathbb{R}}_{\geq 0}^{n\times m}$. The utility of opening
facilities in a subset $S$ of the locations is $u(S)$, where
$u(S)=\sum_{b=1}^{m}\max_{a\in S}M_{a,b}$. This corresponds to the total value
obtained by all the customers, assuming each customer chooses the open
facility with highest service value. The problem of Krause and Golovin is to
maximize the utility function $u$ subject to a constraint on the number of
facilities that can be opened.

We introduce a min-sum version of this facility location problem by
considering the Min-Sum Submodular Cover problem with the utility function $u$
just described, and with $c_{i}$ representing the time to open a facility $i$.
This problem corresponds to a situation where facilities will be opened in all
$n$ locations, but they can only be opened one at a time. $M_{a,b}$ represents
the estimated value facility $a$ will provide to customer $b$ per unit of
time, once facility $a$ is opened. Minimizing the objective value
$\sum_{i=1}^{n}c(S_{i})\left(u(S_{i})-u(S_{i-1})\right)$ corresponds to
finding the order to build facilities so as to minimize lost value as
facilities are built.

#### Our Contributions

We introduce the study of solving Min-Sum Submodular Cover using local search.
Building on work of Munagala et al. [16], who presented a local-search
algorithm for Pipelined Set Cover, we generalize their LP-based analysis by
redefining a key quantity in their proof and using second-order
supermodularity. We show that local search produces a
$(4+\epsilon)$-approximate solution for Min-Sum Submodular Cover in time
$O(n^{3}\log(\frac{d}{\epsilon}))$, assuming second-order supermodularity of
the utility function, when initialized with a $d$-approximate solution. We
prove that a permutation listing the items in non-decreasing cost order is an
$n$-approximate solution. Thus initializing local search with a non-decreasing
cost permutation enables us to reach a $(4+\epsilon)$-approximate solution in
time $O(n^{3}\log(\frac{n}{\epsilon}))$. Applying this result to Pipelined Set
Cover improves on the $O(n^{3}\log(\frac{mn}{\epsilon}))$ time bound from
Munagala et al., where $m$ is the size of the ground set of the set cover
instance, by eliminating the dependence on $m$. We also present results of
experiments on two types of Min-Sum Submodular Cover problems: Pipelined Set
Cover and Min-Sum Facility Location. Our empirical findings suggest that local
search can reliably produce better solutions than the natural greedy algorithm
on small data sets.

##  2 Preliminaries

Let $f(e|S):=f(S\cup\\{e\\})-f(S)$ be the marginal utility of adding element
$e$ to set $S$. With this notation in hand, we define several useful
properties of set functions.

######  Definition 2.1 (Set Function Properties).

Consider a positive integer $n$ and set function
$f:2^{[n]}\rightarrow\operatorname*{\mathbb{R}}_{\geq 0}$. We first define the
following properties of $f$, which hold if the inequality given below for that
property holds for all $S\subseteq[n]$ and all $i,j,k\in[n]\setminus S$,

  * •

monotone: $f(S\cup\\{i\\})\geq f(S),$

  * •

submodular (diminishing returns): $f(i|S)\geq f(i|S\cup\\{j\\}),$

  * •

second-order supermodular: $f(i|S)-f(i|S\cup\\{j\\})\geq
f(i|S\cup\\{k\\})-f(i|S\cup\\{k,j\\}).$

Note that the way in which we have written the above properties illustrates
the progression from monotonicity to submodularity and submodularity to
second-order supermodularity: we arrive at the ‘next’ property by subtracting
the left-hand side from the right-hand side. Another related property is
modularity: for all $S\subseteq[n]$, $f(S)=\sum_{i\in S}f(\\{i\\})$. In this
paper, the functions we consider will be monotone set functions that are
normalized, i.e., $f(\emptyset)=0$ unless otherwise stated.

The Min-Sum Submodular Cover problem is a special case of the Min-Sum
Permutation Problem, defined by Happach et al. [8]. That problem has the same
objective function as Min-Sum Submodular Cover, and minimization may be over
all permutations, or only over a subset of them. The only assumptions on $u$
and $c$ in [8] are that they are monotone and normalized.

##  3 A Local Search Algorithm for Min-Sum Submodular Cover

Munagala et al. [16] gave a local search algorithm for the special case of the
Min-Sum Submodular Cover problem where $u$ is a coverage function. Applying
the same approach to the general Min-Sum Submodular Cover problem, we have the
following local search algorithm: initialize the algorithm with a permutation
$\pi$ of $[n]$. Define a neighbor $\pi^{\prime}$ of $\pi$ to be a permutation
that can be be produced from $\pi$ by removing the element in some position
$i$ of $\pi$ and reinserting it in position $j$. Find the neighbor
$\pi^{\prime}$ of $\pi$ with lowest objective value (given by Equation 1). If
that value is less than the objective value of $\pi$, then replace $\pi$ by
$\pi^{\prime}$ and repeat. Otherwise, output $\pi$. Pseudocode for this
algorithm is given in Algorithm LABEL:alg:local.

The analysis of Munagala et al. [16] shows that in the special case where $u$
is a coverage function, Algorithm LABEL:alg:local achieves a
$(4+\epsilon)$-approximation to the optimal permutation. We generalize their
analysis to all utility functions $u$ that are submodular and second-order
supermodular (in addition to being monotone and normalized, which we assume is
the case for all utility functions in this paper).

Algorithm 1: Local search algorithm to produce a $(4+\epsilon)$-approximation.

⬇

Input: $\epsilon>0$, $n>0$, utility function
$u:2^{[n]}\rightarrow\operatorname*{\mathbb{R}}_{\geq 0}$,

cost function $c:2^{[n]}\rightarrow\operatorname*{\mathbb{R}}_{\geq 0}$,
$d$-approximate permutation $\pi$

Output: permutation $\pi$

for iteration in $\\{1,\ldots,2n^{3}\log(d/\epsilon)\\}$ do

$\pi^{*}\leftarrow\pi$

for $i,j\in[n]$ do

# $\pi^{\prime}$ is $\pi$ with $\pi(i)$ moved to position $j$

$\pi^{\prime}\leftarrow move(\pi,i,j)$

# objective$(u,c,\pi)$ is Equation (1)

if objective$(u,c,\pi^{\prime})<$ objective$(u,c,\pi^{*})$ do

$\pi^{*}\leftarrow\pi^{\prime}$

if $\pi^{*}=\pi$

# Algorithm converged:

# $\pi$ is locally optimal with respect to moves

return $\pi$ # $4$-approximation

$\pi\leftarrow\pi^{*}$

return $\pi$ # $(4+\epsilon)$-approximation

Let $\pi_{c}$ be a non-decreasing cost permutation, i.e.,
$c(\\{\pi_{c}(i)\\})\leq c(\\{\pi_{c}(i+1)\\})$ for $i\in[n-1]$. We prove the
following results.

######  Theorem 3.1.

Fix a positive integer $n$. Let
$u:2^{[n]}\rightarrow\operatorname*{\mathbb{R}}_{\geq 0}$ be a submodular and
second-order supermodular set function and let
$c:2^{[n]}\rightarrow\operatorname*{\mathbb{R}}_{\geq 0}$ be a modular set
function. If Algorithm LABEL:alg:local converges before terminating, then the
solution it returns is a $4$-approximation to Min-Sum Submodular Cover on $u$
and $c$.

Unfortunately, we cannot guarantee that Algorithm LABEL:alg:local will
converge before terminating. The next result guarantees a
$(4+\epsilon)$-approximation when the algorithm terminates.

######  Theorem 3.2.

Consider the positive integer $n$, utility function $u$, and cost function $c$
considered in Theorem 3.1. Fix $\epsilon>0$. Let $\pi$ be a $d$-approximate
permutation. If Algorithm LABEL:alg:local does not converge before
terminating, then the solution it returns (after $2n^{3}\log(d/\epsilon)$
iterations), is a $(4+\epsilon)$-approximation to Min-Sum Submodular Cover on
$u$ and $c$.

Assuming constant access query access to $u$ and $c$, Algorithm
LABEL:alg:local returns a $(4+\epsilon)$-approximation in
$2n^{3}\log(d/\epsilon)$ time.

As in Munagala et al. [16], in our analysis we consider a modified version of
local search based on “insertions” rather than “moves.” We find it easier to
analyze local search with insertions and the approximation result immediately
applies to local search with moves since a permutation that is locally optimal
with respect to moves is also locally optimal with respect to insertions. In
each iteration of local search with insertions, rather than considering the
set of neighbors $\pi^{\prime}$ of $\pi$, the modified algorithm considers a
set of what we will call pseudo-neighbors. Each is derived from $\pi$ by
taking an element appearing in some position $i$ of $\pi$, and inserting a
second copy of the element into some position $j<i$. Each pseudo-neighbor of
$\pi$ corresponds to a unique neighbor of $\pi$, produced from the pseudo-
neighbor by removing the original copy of the repeated element (which appears
closer to the end of the permutation).

Define the objective value of a pseudo-neighbor $\pi^{\prime}$ (which has
length $n+1$) to be $\sum_{i=1}^{n+1}c^{\prime}(S_{i})[u(S_{i})-u(S_{i-1})]$,
where here $S_{i}$ is the prefix of $\pi^{\prime}$ containing its first $i$
elements, $u(S_{i})$ is the value of $u$ for the set of distinct items in
$S_{i}$, and $c^{\prime}(S_{i})=\sum_{j=1}^{i}c(\\{s_{j}\\})$ where $s_{j}$ is
the element in position $j$ of $\pi^{\prime}$. That is, if both copies of the
repeated element appear within the first $i$ positions of $\pi^{\prime}$, then
$c^{\prime}(S_{i})$ charges for both copies.

If the objective value of $\pi$ is no greater than the value of its pseudo-
neighbors, then the modified algorithm outputs $\pi$. Otherwise, the algorithm
takes the pseudo-neighbor with lowest objective value, deletes the original
copy of its repeated element, and uses the resulting permutation as the new
value of $\pi$ in the next iteration.

The objective value of a pseudo-neighbor of $\pi$ is clearly greater than or
equal to the objective value of the corresponding neighbor. Therefore, if
$\pi$ has no neighbor with lower objective value, then it has no pseudo-
neighbor with lower objective value. It follows that the bounds we prove on
the modified local search algorithm (with insertions) also apply to the
original local search algorithm (with moves).

We prove Theorems 3.1 and 3.2 in the remainder of this section.

###  3.1 Proof of Theorem 3.1: $4$-approximation

Say a permutation $\pi$ is locally optimal if no pseudo-neighbor has lower
objective value. We begin by proving that a locally optimal solution satisfies
a certain inequality, expressed in terms of variables $b_{ij}$. This
inequality is taken from the analysis in Munagala et al. [16], but we define
the variables $b_{ij}$ differently here. We will use the following technical
Observation 3.3 to prove Lemma 3.5.

######  Observation 3.3.

Consider three sequences of non-negative real numbers, $X_{0},\ldots,X_{n}$,
$Y_{0},\ldots,Y_{n}$, and $C_{0},\ldots,C_{n}$. Let $j\in[n]$. Suppose that
the following hold: (i) $X_{0}\geq Y_{0}$, (ii) for all
$r\in\\{j,\ldots,n\\}$, $C_{0}\leq C_{r}$ and $X_{r}\leq Y_{r}$, and (iii)
$X_{0}+\sum_{r=j}^{n}X_{r}=Y_{0}+\sum_{r=j}^{n}Y_{r}$. Then

$\displaystyle C_{0}X_{0}+\sum_{r=j}^{n}C_{r}X_{r}\leq
C_{0}Y_{0}+\sum_{r=j}^{n}C_{r}Y_{r}.$

######  Proof 3.4.

Rewriting the final assumption yields

$\displaystyle X_{0}-Y_{0}=\sum_{r=j}^{n}(Y_{r}-X_{r})\iff
C_{0}(X_{0}-Y_{0})=\sum_{r=j}^{n}C_{0}(Y_{r}-X_{r})$ $\displaystyle\Rightarrow
C_{0}(X_{0}-Y_{0})\leq\sum_{r=j}^{n}C_{r}(Y_{r}-X_{r})$

where the second implication follows from the non-negativity of $Y_{r}-X_{r}$
and the assumption that $C_{0}\leq C_{r}$. Observation 3.3 follows
immediately.

######  Lemma 3.5.

Suppose $u:2^{[n]}\rightarrow\operatorname*{\mathbb{R}}_{\geq 0}$ is a
submodular and second-order supermodular set function and
$c:2^{[n]}\rightarrow\operatorname*{\mathbb{R}}_{\geq 0}$ is a modular set
function. Let $L_{j}$ denote the first $j$ elements of the locally optimal
permutation and $O_{i}$ denote the first $i$ elements of the optimal
permutation. Similarly, we use $l_{j}$ to represent the $j$th element of the
local permutation and $o_{i}$ to represent the $i$th element of the optimal
permutation. Then

$\displaystyle\sum_{r=j}^{n}c(L_{r})\sum_{s=1}^{n}b_{sr}\leq[c(o_{i})+c(L_{j-1})]\sum_{r=j}^{n}b_{ir}+\sum_{r=j}^{n}[c(o_{i})+c(L_{r})],\sum_{\begin{subarray}{c}s=1\\\
s\neq i\end{subarray}}^{n}b_{sr}$ (2)

where

$\displaystyle b_{ij}=u(o_{i}|O_{i-1}\cup L_{j-1})-u(o_{i}|O_{i-1}\cup
L_{j-1}\cup\\{l_{j}\\})$ $\displaystyle=u(l_{j}|O_{i-1}\cup
L_{j-1})-u(l_{j}|O_{i-1}\cup L_{j-1}\cup\\{o_{i}\\}).$

######  Proof 3.6.

When $u$ is a coverage function, as in the analysis in Munagala et al. [16],
$b_{ij}$ represents the number of ground elements covered in the optimal
permutation by subset $o_{i}$ (and not by $o_{1},\ldots,o_{i-1}$) and in the
local permutation by subset $l_{j}$ (and not by $l_{1},\ldots,l_{j-1})$. Our
definitions of $b_{ij}$ generalize this intuition to functions where the
utility is more abstract. In particular, we can use telescoping sums to derive
the following identities:

$\displaystyle\sum_{i=1}^{n}b_{ij}=u(l_{j}|L_{j-1})\hskip
20.00003pt\sum_{r=j}^{n}b_{ir}=u(o_{i}|O_{i-1}\cup L_{j-1})\hskip
20.00003pt\sum_{\begin{subarray}{c}s=1\\\ s\neq
i\end{subarray}}^{n}b_{sr}=u(l_{r}|L_{r-1})-b_{ir}$ (3)

Then Equation (2) is equivalent to

$\displaystyle\sum_{r=j}^{n}c(L_{r})$ $\displaystyle
u(l_{r}|L_{r-1})\leq[c(o_{i})+c(L_{j-1})]{u(o_{i}|O_{i-1}\cup
L_{j-1})}+\sum_{r=j}^{n}[c(o_{i})+c(L_{r})][u(l_{r}|L_{r-1})-b_{ir}].$ (4)

Notice that Equation (4) is not quite equivalent to the property of local
optimality with respect to insertions. Instead, we know the following (very
similar) inequality holds by the property that $L$ is locally optimal:

$\displaystyle\sum_{r=j}^{n}c(L_{r})$ $\displaystyle
u(l_{r}|L_{r-1})\leq[c(o_{i})+c(L_{j-1})]{u(o_{i}|L_{j-1})}+\sum_{r=j}^{n}[c(o_{i})+c(L_{r})]{u(l_{r}|L_{r-1}\cup\\{o_{i}\\})}.$
(5)

We will now show that the right-hand side of Equation (5) lower bounds the
right-hand side of Equation (4). Then Equation (4) follows from Equation (5).
We do this through the following four conditions combined with Observation
3.3:

$\displaystyle c(L_{j-1})\leq c(L_{r-1})$ (6) $\displaystyle
u(o_{i}|L_{j-1})\geq u(o_{i}|O_{i-1}\cup L_{j-1})$ (7) $\displaystyle
u(l_{r}|L_{r-1}\cup\\{o_{i}\\})\leq u(l_{r}|L_{r-1})-b_{ir}$ (8)
$\displaystyle
u(o_{i}|L_{j-1})+\sum_{r=j}^{n}{u(l_{r}|L_{r-1}\cup\\{o_{i}\\})}={u(o_{i}|O_{i-1}\cup
L_{j-1})}+\sum_{r=j}^{n}[u(l_{r}|L_{r-1})-b_{ir}]$ (9)

for all $i,j,r\in[n]$ with $j\leq r\leq n$. Equation (6) holds by the
monotonicity of $c$ and Equation (7) holds by the submodularity of $u$.
Equation (8) holds because $u$ is second-order supermodular. (This is where we
use second-order supermodularity.)

Finally, Equation (9) holds by the following argument: Notice that the left-
hand side is equal to $u([n])-u(L_{j-1})$ since we sequentially sum the
marginal utility of adding the next element to our current chain. Similarly,
the right-hand side simplifies to $u([n])-u(L_{j-1})$ after cancellation using
Equation (3).

Using the above, we apply Observation 3.3 with $C_{0}=c(o_{i})+c(L_{j-1})$,
$C_{r-j+1}=c(o_{i})+c(L_{r-1})$, $X_{0}=u(o_{i}|L_{j-1})$,
$X_{r-j+1}=u(l_{r}|L_{r-1}\cup\\{o_{i}\\})$, $Y_{0}=u(o_{i}|O_{i-1}\cup
L_{j-1})$, and $Y_{r-j+1}=u(l_{r}|L_{r-1})-b_{ir}$. This yields

$\displaystyle[c(o_{i})+c(L_{j-1})]u(o_{i}|L_{j-1})+\sum_{r=j}^{n}[c(o_{i})+c(L_{r-1})][u(l_{r}|L_{r-1}\cup\\{o_{i}\\})]$
$\displaystyle\leq[c(o_{i})+c(L_{j-1})]u(o_{i}|O_{i-1}\cup
L_{j-1})+\sum_{r=j}^{n}[c(o_{i})+c(L_{r-1})][u(l_{r}|L_{r-1})-b_{ir}].$

Then the right-hand side of Equation (4) must be greater than or equal to the
right-hand side of Equation (5). Therefore Lemma 3.5 follows by Equation (5).

######  Proof 3.7 (Proof of Theorem 3.1).

We will now prove that a locally optimal permutation with respect to
insertions is (at worst) a 4-approximation of the optimal permutation. To do
this, we will show that the same linear program used by Munagala et al. [16]
to prove their 4-approximation for coverage functions also applies to all
utility functions that are monotone, submodular and second-order supermodular.
The linear program in Munagala et al. [16] is:

$\displaystyle\textrm{maximize }\sum_{j=1}^{n}c(L_{j})\sum_{i=1}^{n}b_{ij}$
(10) $\displaystyle\textnormal{ subject to
}\sum_{i=1}^{n}c(O_{i})\sum_{j=1}^{n}b_{ij}\leq 1\hskip
10.00002pt\textnormal{and}$
$\displaystyle\sum_{r=j}^{n}c(L_{r})\sum_{s=1}^{n}b_{sr}\leq(c(o_{i})+c(L_{j-1}))\sum_{r=j}^{n}b_{ir}+\sum_{r=j}^{n}(c(o_{i})+c(L_{r}))\sum_{\begin{subarray}{c}s=1\\\
s\neq i\end{subarray}}^{n}b_{sr}\hskip 10.00002pt\forall i,j\in[n]$ (11)

where the variables $b_{ij}$ are non-negative for all $i,j\in[n]$.

The first constraint scales the variables in the optimal permutation so that
the objective is the ratio of the local permutation to the optimal
permutation. By Lemma 3.5, the second constraint must hold for locally optimal
$L$, with the given values for $b_{ij}$. As shown in Munagala et al. [16],
there is a feasible solution to the dual of the LP with objective value 4 so,
by strong duality, the locally optimal permutation achieves a 4-approximation
to the optimal permutation.

###  3.2 Proof of Theorem 3.2: $(4+\epsilon)$-approximation

In this section, we show that local search reaches a
$(4+\epsilon)$-approximation in reasonable time. The proof is based on the
following fact, which bounds the progress made in each “round” of the local
search (i.e., in each iteration of the while loop in Algorithm
LABEL:alg:local).

######  Fact 1.

Let $\pi$ be a permutation of $[n]$ that is an $M$-approximation to Min-Sum
Submodular Cover. Then applying one round of local search to this permutation
will either establish that it is a local optimum or yield a new permutation
that is an $(M-\frac{M-4}{2n})$-approximation.

######  Proof 3.8.

In Munagala et al. [16], they prove that Fact 1 holds for a Pipelined Filter
Ordering instance, i.e., when $u$ is a coverage function. Consider an
$M$-approximate solution. Let $A$ be a variable representing the reduction in
approximation factor after making the local search step. As in Equation (11),
the variables $b_{ij}$ are non-negative real numbers and correspond to a
concept of shared utility formalized in Lemma 3.5. Then the following linear
program is what Munagala et al. [16] use to lower-bound the improvement in
approximation ratio:

$\displaystyle\textrm{minimize }A$ $\displaystyle\textrm{ subject to
}\sum_{i=1}^{n}c(O_{i})\sum_{j=1}^{n}b_{ij}\leq 1\textrm{ and }$
$\displaystyle\sum_{r=j}^{n}c(L_{r})\sum_{s=1}^{n}b_{sr}\leq
A+(c(o_{i})+c(L_{j-1}))\sum_{r=j}^{n}b_{ir}+\sum_{r=j}^{n}(c(o_{i})+c(L_{r}))\sum_{\begin{subarray}{c}s=1\\\
s\neq i\end{subarray}}^{n}b_{sr}\textrm{ and }$
$\displaystyle\sum_{j=1}^{n}c(L_{j})\sum_{i=1}^{n}b_{ij}\geq M\hskip
10.00002pt\forall i,j\in[n].$

For the more general case where $u$ is an abstract submodular and second-order
supermodular utility function, the first and third inequalities trivially hold
for our generalized definition of $b_{ij}$. The second inequality follows from
Lemma 3.5. By taking the dual, Munagala et al. [16] show that the objective of
the primal and therefore the reduction in approximation ratio is at least
$\frac{M-4}{2n}$ which gives Fact 1. When $u$ is an arbitrary submodular and
second-order supermodular function, it follows from Lemma 3.5 that the same
linear program still lower-bounds the improvement.

To achieve a good bound on the time required for the local search for Min-Sum
Submodular Cover, we want to begin the local search from a permutation that is
not too far from optimal. The following theorem gives such a permutation.

######  Theorem 3.9 (Non-Decreasing Cost).

A permutation that orders the elements $i\in[n]$ in non-decreasing order of
$c(\\{i\\})$ is an $n$-approximate solution to Min-Sum Submodular Cover.

######  Proof 3.10.

Suppose, without loss of generality, that $c(\\{1\\})\leq\cdots\leq
c(\\{n\\})$.

Suppose an optimal solution is given by the sets
$\emptyset=S_{0},S_{1}\ldots,S_{n}$, and for each $j\in[n]$, let
$S^{\prime}_{j}=[\max S_{j}]$, where $\max S_{j}$ is the maximum integer in
$S_{j}$. For instance, for $n=4$ and $S_{1}=\\{2\\}$, $S_{2}=\\{12\\}$,
$S_{3}=\\{124\\}$, $S_{4}=\\{1234\\}$, we have $S^{\prime}_{1}=\\{12\\}$,
$S^{\prime}_{2}=\\{12\\}$, $S^{\prime}_{3}=\\{1234\\}$,
$S^{\prime}_{4}=\\{1234\\}$. Setting $[0]$ and $S^{\prime}_{0}$ equal to
$\emptyset$, observe that

$\displaystyle\sum_{j=1}^{n}c([j])(u([j])-u([j-1]))\leq\sum_{j=1}^{n}c(S^{\prime}_{j})(u([j])-u([j-1]))=\sum_{j=1}^{n}c(S^{\prime}_{j})(u(S^{\prime}_{j})-u(S^{\prime}_{j-1})),$

where the inequality follows by the monotonicity of $c$ and the equality
follows (not term-wise but in total) by considering when utility is accrued by
each side.

Now, for each $j\geq 0$, we have $c(S_{j})\geq c(\\{\max S_{j}\\})$, by the
monotonicity of $c$. Also, since $c(\\{\max S_{j}\\})\geq c(\\{i\\})$ for
every $i\in S_{j}$ by the indexing assumption,

$n\cdot c(S_{j})\geq n\cdot c(\\{\max S_{j}\\})=n\cdot c(\\{\max
S^{\prime}_{j}\\})\geq c(S^{\prime}_{j})$

where the equality follows from the definition of $S^{\prime}_{j}$ and the
second inequality follows since there are at most $n$ elements in
$S^{\prime}_{j}$. Then the objective value of the optimal permutation is at
least ${\frac{1}{n}\sum_{j=1}^{n}c(S^{\prime}_{j})(u(S_{j})-u(S_{j-1}))}$, or
equivalently by charging utility to each increase in cost,

$\frac{1}{n}\sum_{j=1}^{n}(u([n])-u(S_{j-1}))(c(S^{\prime}_{j})-c(S^{\prime}_{j-1})).$

By the monotonicity of $u$, this sum is at least

$\frac{1}{n}\sum_{j=1}^{n}(u([n])-u(S^{\prime}_{j-1}))(c(S^{\prime}_{j})-c(S^{\prime}_{j-1}))\Leftrightarrow\frac{1}{n}\sum_{j=1}^{n}c(S^{\prime}_{j})(u(S^{\prime}_{j})-u(S^{\prime}_{j-1})).$

This is at least $1/n$ times the objective value of the increasing cost
permutation, by our earlier observation.

We note that Theorem 3.9 also holds for Min-Sum Permutation Problems
minimizing over all permutations where $u$ only satisfies monotonicity and $c$
only satisfies monotonicity and subadditivity (both are still normalized).

We can now prove that the output of Algorithm LABEL:alg:local is a
$(4+\epsilon)$-approximation to Min-Sum Submodular Cover. The time bound
assumes constant time oracle queries.

######  Proof 3.11 (Proof of Theorem 3.2).

The improvement in the quality of the solution in each round of local search,
guaranteed by Fact 1, implies that a $(4+\epsilon)$-approximation is achieved
within $O(n\log(\frac{d}{\epsilon}))$ rounds of Algorithm LABEL:alg:local,
when it is initialized with a $d$-approximate permutation. This implication
was stated without proof by Munagala et al. [16], in proving the same bounds
for Pipelined Set Cover. For completeness, we present a proof of the
implication in Appendix A.2. The $O(n^{3}\log(\frac{d}{\epsilon}))$ time bound
is achieved by spending $O(n^{2})$ per round. To accomplish this, we do not
recompute all $n$ terms of the objective function for each of the
$\theta(n^{2})$ neighbors of the current solution. Instead, by considering
‘adjacent’ neighbors sequentially, we can compute the objective function value
for the next neighbor from the value obtained for the previous neighbor in
constant time, by recomputing only two terms of the objective function. The
time bound for the non-decreasing cost permutation follows from Theorem 3.9.

##  4 Experiments

Assuming constant-time oracle access to the utility functions, the greedy
algorithm runs in $O(n^{2})$ total time, while our local search algorithm
spends time $O(n^{2})$ in each round. In our experiments, with no oracle, we
had to compute the utility function.

The greedy algorithm is certainly faster than the local search, but it only
explores one type of solution, where cost effective elements appear earlier in
the permutation. Local search initialized with random permutations can sample
from the entire solution space. Our experiments compare the greedy solution to
local search solutions from four random initial permutations, and from a non-
decreasing cost permutation. We run each local search for $n$ steps (rather
than running it to convergence, or until it is guaranteed to find a
$(4+\epsilon)$-approximate solution). We see empirically that the best of the
5 local search solutions tends to be better than the worst, and also better
than the greedy solution. In applications where computing is cheap, $n$ is not
large, and the quality of the solution is crucial, using local search may be
preferable to using greedy.

Our experiments compare greedy and local search on 100 random instances of two
problems: Pipelined Set Cover and Min-Sum Facility Location. For our random
instances, we set $n=30$ and each cost $c_{i}$ to be a uniform random value
between 0 and 1. Figure 1 shows the results of our experiments, with objective
values given relative to the best of the 6 solutions (1 greedy and 5 local
search). Local search finds the best solution in almost $100\%$ of the 100
instances whereas greedy finds the best in roughly $50\%$.

Figure 1: Histograms of greedy and local search performance with $n=30$ in
Pipelined Set Cover and Min-Sum Facility Location. The relative objective
values are normalized with respect to the best of the 6 generated solutions (1
greedy and 5 local search). Frequency is reported from the 100 random
instances generated for each dataset.

#### Pipelined Set Cover

We perform experiments on synthetic, randomly generated instances of the
(unweighted) Pipelined Set Cover problem with correlated subsets, following an
approach of Babu et al. [4]. Recall that an instance of (unweighted) Pipelined
Set Cover consists of a finite ground set $[m]$ and a collection of subsets
$D_{i}\subseteq V$ for $i\in[n]$. The utility of a set $S\subseteq[n]$ is
$u(S)=|\bigcup_{i\in S}D_{i}|$. The $n$ subsets $D_{i}$ in our random instance
are divided into $\lceil n/\Gamma\rceil$ groups, where $\Gamma$ is a
“correlation factor.” The instance has the following properties, for each
element $j\in[m]$. For all $i\in[n]$, $\Pr[j\in D_{i}]$ is a fixed value $p$.
For two subsets $D_{i}$ and $D_{i^{\prime}}$ in different groups, membership
of $j$ in $D_{i}$ is independent of its membership in $D_{i^{\prime}}$. For
two subsets $D_{i}$ and $D_{i^{\prime}}$ in the same group, the probability
that $j$ has the same membership status in $D_{i}$ and $D_{i^{\prime}}$ (i.e.,
is either in both subsets, or in neither), is a fixed value $\rho$. In our
experiments, $m=2n$, $\Gamma=4$, $p=0.3$, and $\rho=0.7858$. In Appendix A.3,
we describe in detail how we generated the instance.

#### Min-Sum Facility Location

We use the locations of $n$ Citi Bike stations [1] in New York City as the
facilities for our facility location data set. For calculating the utility to
customers, we uniformly generate customer locations within the range of
latitude and longitude of the stations. The value $M_{a,b}$ for customer $b$
and station $a$ is the inverse of the Euclidean distance between them.

## References

  * [1] Citi bike system data.  https://ride.citibikenyc.com/system-data, 2021.  Accessed: 2021-12-1.
  * [2] Emile Aarts, Emile HL Aarts, and Jan Karel Lenstra.  Local search in combinatorial optimization.  Princeton University Press, 2003.
  * [3] Daniel Antunes, Claire Mathieu, and Nabil H. Mustafa.  Combinatorics of Local Search: An Optimal 4-Local Hall’s Theorem for Planar Graphs.  In Kirk Pruhs and Christian Sohler, editors, 25th Annual European Symposium on Algorithms (ESA 2017), volume 87 of Leibniz International Proceedings in Informatics (LIPIcs), pages 8:1–8:13, Dagstuhl, Germany, 2017. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik.  URL: http://drops.dagstuhl.de/opus/volltexte/2017/7829, doi:10.4230/LIPIcs.ESA.2017.8.
  * [4] Shivnath Babu, Rajeev Motwani, Kamesh Munagala, Itaru Nishizawa, and Jennifer Widom.  Adaptive ordering of pipelined stream filters.  In Proceedings of the 2004 ACM SIGMOD International Conference on Management of Data, pages 407–418, 2004.
  * [5] Uriel Feige, László Lovász, and Prasad Tetali.  Approximating min sum set cover.  Algorithmica, 40(4):219–234, 2004.
  * [6] Mehrdad Ghadiri, Richard Santiago, and Bruce Shepherd.  Beyond submodular maximization.  arXiv preprint arXiv:1904.09216, 2019.
  * [7] Mehrdad Ghadiri, Richard Santiago, and Bruce Shepherd.  A parameterized family of meta-submodular functions.  arXiv preprint arXiv:2006.13754, 2020.
  * [8] Felix Happach, Lisa Hellerstein, and Thomas Lidbetter.  A general framework for approximating min sum ordering problems.  INFORMS Journal on Computing, 2021.
  * [9] Satoru Iwata, Prasad Tetali, and Pushkar Tripathi.  Approximating minimum linear ordering problems.  In Proceedings of Approx-Random, 2012.
  * [10] Rishabh Iyer, Ninad Khargoankar, Jeff Bilmes, and Himanshu Asanani.  Submodular combinatorial information measures with applications in machine learning.  In Algorithmic Learning Theory, pages 722–754. PMLR, 2021.
  * [11] Rishabh Iyer, Ninad Khargonkar, Jeff Bilmes, and Himanshu Asnani.  Generalized submodular information measures: Theoretical properties, examples, optimization algorithms, and applications.  IEEE Transactions on Information Theory, 2021.
  * [12] Nitish Korula, Vahab Mirrokni, and Morteza Zadimoghaddam.  Online submodular welfare maximization: Greedy beats 1/2 in random order.  SIAM Journal on Computing, 47(3):1056–1086, 2018.
  * [13] Andreas Krause and Daniel Golovin.  Submodular function maximization.  In Tractability: Practical Approaches to Hard Problems. Cambridge University Press, 2014.
  * [14] Gilad Kutiel and Dror Rawitz.  Local Search Algorithms for Maximum Carpool Matching.  In Kirk Pruhs and Christian Sohler, editors, 25th Annual European Symposium on Algorithms (ESA 2017), volume 87 of Leibniz International Proceedings in Informatics (LIPIcs), pages 55:1–55:14, Dagstuhl, Germany, 2017. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik.  URL: http://drops.dagstuhl.de/opus/volltexte/2017/7821, doi:10.4230/LIPIcs.ESA.2017.55.
  * [15] Wenjun Li, Jianer Chen, and Jianxin Wang.  Deeper local search for better approximation on maximum internal spanning trees.  In European Symposium on Algorithms, pages 642–653. Springer, 2014\.
  * [16] Kamesh Munagala, Shivnath Babu, Rajeev Motwani, and Jennifer Widom.  The pipelined set cover problem.  In International Conference on Database Theory, pages 83–98. Springer, 2005.
  * [17] Matthew Streeter and Daniel Golovin.  An online algorithm for maximizing submodular functions.  Technical report, Carnegie-Mellon University School of Computer Science, 2007.
  * [18] Matthew J. Streeter and Daniel Golovin.  An online algorithm for maximizing submodular functions.  In Neural Information Processing Systems (NeurIPS), pages 1577–1584, 2008.

##  Appendix A Appendix

###  A.1 Prior Use of the Term ‘Min-Sum Submodular Cover’

The first proof that a greedy algorithm for Min-Sum Submodular Cover yields a
4-approximation was given by Streeter and Golovin. The proof appears both in a
Technical Report [17] and an associated conference paper [18]. They actually
gave their proof for a more general problem than the one we considered in this
paper, in which
$u:2^{[n]\times\operatorname*{\mathbb{R}}}\rightarrow\operatorname*{\mathbb{R}}$,
and the output is a sequence of pairs of the form
$(i,\tau)\in[n]\times\operatorname*{\mathbb{R}}$. In their Technical Report,
Streeter and Golovin used the name Min-Sum Submodular Cover to refer to the
more general problem, but they did not use this name (nor any other) to refer
to the problem in their conference paper. We opted to use the name Min-Sum
Submodular Cover to refer to the problem we defined in this paper, as we
believe this usage of the name is natural given the connection to Min-Sum
Submodular Cover.

We note that the definition Streeter and Golovin gave for the more general
problem is problematic as written. The greedy algorithm may not be well-
defined for functions $u$ that are non-zero for subsets of
$[n]\times\operatorname*{\mathbb{R}}$ that include pairs $(i,\tau)$, where
$\tau$ is infinitesimally small. However, the results in the paper are not
dependent on allowing such $u$, and the problem with the definition can be
fixed by restricting the domain of $u$.

###  A.2 Bounding the number of rounds in the proof of Theorem 3.2

We show using Fact 1 that Algorithm LABEL:alg:local yields a
$(4+\epsilon)$-approximation in at most $O(n\log(\frac{d}{\epsilon}))$ rounds
from a $d$-approximate permutation.

We introduce a recurrence relation $T(\ell)=aT(\ell-1)-b$ where $a=2n/(2n-1)$
and $b=4/(2n-1)$. We can derive this recurrence by setting $T(\ell)=M$ and
$T(\ell-1)=M-(M-4)/2n$. Intuitively, $\ell$ is the number of iterations until
we reach $(4+\epsilon)$ and so we set $T(0)=4+\epsilon$. By repeatedly
expanding $T(\ell)$, we get

$\displaystyle
T(\ell)=-b\sum_{i=0}^{\ell-1}a^{i}+a^{\ell}(4+\epsilon)=-b\frac{a^{\ell}-1}{a-1}+a^{\ell}(4+\epsilon)=\epsilon\left(\frac{2n}{2n-1}\right)^{\ell}+4$

where the last equality follows by plugging in the values of $a$ and $b$. We
claim that a $d$-approximate permutation is at most
$2n\log(\frac{d}{\epsilon})$ rounds from a $(4+\epsilon)$-approximation. We
can verify this by evaluating

$\displaystyle
T\left(2n\log\left(\frac{d}{\epsilon}\right)\right)=\epsilon\left(\frac{2n}{2n-1}\right)^{2n\log(\frac{d}{\epsilon})}+4=d^{2n\log\left(\frac{2n}{2n-1}\right)}+4.$

Notice that

$\displaystyle x\log\left(\frac{x}{x-1}\right)\geq 1\iff
e^{x}\left(\frac{x}{x-1}\right)\geq e$

which is certainly true for $x>2$. It follows that
$T(2n\log(\frac{d}{\epsilon}))\geq d$ so local search converges in at most
$2n\log(\frac{d}{\epsilon})$ rounds.

###  A.3 Generation of the Pipelined Set Cover Data Set

We generate the $n$ subsets $D_{1},\ldots,D_{n}$ of $[m]$ randomly as follows.
Initially, for each group $G$ of subsets, we generate an advice bit $a_{G,j}$
for each element $j\in[m]$, which is True with probability $p$, and False with
probability $1-p$. Then for each element $j$, and for each $D_{i}$ in group
$G$, we do the following: with probability $p^{\prime}$ we use advice bit
$a_{G,j}$ to determine whether or not to include $j$ in $D_{i}$ (if $a_{G,j}$
is True, we include $j$ in $D_{i}$, else we do not). With probability
$(1-p^{\prime})$, we ignore the advice bit, and instead, we include $j$ in
$D_{i}$ with probability $p$, and exclude it with probability $(1-p)$. The
probability that $j$ is in a given set $D_{j}$ is clearly $p$.

For $i\neq i^{\prime}$, if subsets $D_{i}$ and $D_{i^{\prime}}$ are in
different groups, then membership of an element $j$ in $D_{i}$ is clearly
independent of its membership of $D_{i^{\prime}}$. If $D_{i}$ and
$D_{i^{\prime}}$ are subsets in the same group $G$, then the probability that
$j$ has the same membership status in both subsets can be calculated by noting
that this can happen in two ways: either the advice bit $a_{G,j}$ was used to
determine membership in both subsets, or $a_{G,j}$ was ignored for one or both
of the two subsets and membership ended up being the same in both subsets. The
first event happens with probability $p^{\prime}\cdot p^{\prime}$. The second
happens with probability $(1-p^{\prime}\cdot p^{\prime})(p\cdot
p+(1-p)\cdot(1-p))$. Because $p=0.3$ and $p^{\prime}=0.7$, the probability
that the membership status of $j$ is the same for both subsets is $.7858$.

We note that it is possible that this process results in subsets $D_{i}$ and
$D_{i^{\prime}}$ where $i\neq i^{\prime}$ and $D_{i}=D_{i^{\prime}}$. We do
not eliminate such duplicates.


# Critical dehydrogenation steps of perhydro-$N$-ethylcarbazole on Ru(0001)
surface

Chunguang Tang1,2<EMAIL_ADDRESS>Preetham Permude1 Shunxin Fei3
Terry J. Frankcombe4 Sean C. Smith2,5 Yun Liu1,2<EMAIL_ADDRESS>1 Research
School of Chemistry, The Australian National University, Canberra, Australia;
2 Institute of Climate, Energy and Disaster Solutions, The Australian National
University, Canberra, Australia;
3 School of Materials Science and Engineering, Anhui University of Technology,
Maanshan 243002, People’s Republic of China;
4 School of Physical, Environmental and Mathematical Sciences, University of
New South Wales, Canberra, Australia.
5 Research School of Physics, The Australian National University, Canberra,
Australia

###### Abstract

Understanding of the critical atomistic steps during the dehydrogenation
process of liquid organic hydrogen carriers (LOHCs) is important to the design
of cost-efficient, high-performance LOHC catalysts. Based on the density
functional theory (DFT) we studied the thermodynamics and kinetics of the
complete dehydrogenation path of perhydro-$N$-ethylcarbazole (12H-NEC) on
Ru(0001) surface, involving the adsorption of 12H-NEC, the discharge of H ions
onto Ru surface, and the desorption of H2 and hydrogen-lean NEC. It was found
that the bonding of $n$H-NEC is significantly strengthened for $n$$\leq$4
because of the flat aromatic ring. Although the whole dehydrogenation process
is endothermic, the release of H from $n$H-NEC, with H adsorbed onto the Ru
surface, was found to be exothermic. The desorption of flat, hydrogen-lean
NEC, which costs $\sim$255 kJ/mol, was identified as the most energy demanding
step. In addition, the effect of surface morphology on adsorption was studied
based on an amorphous surface model. Overall, the results imply more efficient
dehydrogenation could be achieved from relatively weak bonding of NEC to
catalysts, either through engineering catalyst surface (such as surface
defects or smaller catalyst particles) or different catalyst materials. Our
calculations also revealed possible dealkylation at elevated temperatures.

Published version at Computational Materials Science, 229, 112373 (2023)

##  I Introduction

Liquid organic hydrogen carriers (LOHCs) have attracted extensive research
interests sotoodeh_overview_2013 ; preuster_liquid_2017 ;
modisha_prospect_2019  as a potential alternative hydrogen storage approach to
the traditional approaches such as liquid H2, compressed H2 gas and even
competitive circular carriers abdin_large-scale_2021 . Nevertheless, the
large-scale application of LOHCs faces challenges such as abundant, cost-
effective LOHCs as well as relevant catalysts. Currently a number of LOHCs,
including benzene, toluene, naphthalene and $N$-ethylcarbazole (NEC), have
been identified as potential LOHCs and their properties were discussed in some
review or perspective articles modisha_prospect_2019 ; niermann_liquid_2019 ;
abdin_large-scale_2021 . Exploration of new LOHCs, such as those from natural
products, has also been carried out recently tang_natural_2020 . For a given
LOHC material, catalysts significantly impact on its (de)hydrogenation process
and hence represent a key parameter for LOHC performance. To date, a number of
catalysts have been explored for the above-mentioned LOHCs, as summarized in
recent reviews modisha_prospect_2019 . For example, for dehydrogenation of
cyclohexane, the performance of Ni-based biniwale_dehydrogenation_2005 ;
xia_cyclohexane_2017 , Ag-based pande_catalytic_2012 , Pt-based
kariya_efficient_2002 ; kariya_efficient_2003  catalysts with various supports
has been investigated and the synergistic effect from the addition of a second
metal has been reported.

As an ideal tool for investigating atomistic process and complementary to
experiments, first principles modelling based on the density functional theory
(DFT) has been often used to provide fundamental insights into the LOHC-
catalyst reactions shin_thermodynamic_2018 ; yook_density_2020 ;
ouma_insight_2018 . For example, DFT calculations have been used to examine
the structural features of tetrahydrocarbazole adsorption on Pd surfaces and
its preferred dehydrogenation pathways crawford_understanding_2007 . The
calculations of dodecahydro(12H)-NEC, 12H-carbazole, and 12H-fluorene on
Pd(111) surface sotoodeh_dehydrogenation_2012  have revealed a linkage between
the dehydrogenation rate and the adsorption strength. Some researchers have
explored the adsorption of NEC and its hydrogenated states on different Ru
surfaces eblagon_study_2010  and suggested that intermediate 8H-NEC, which is
kinetically stable, prefers surface sites of low coordination (such as edges).
To date most studies have focused on the adsorption of LOHCs on catalyst
surface, but the atomistic-scale (de)hydrogenation process involves multiple
subprocesses, such as (for dehydrogenation) the adsorption of hydrogen-rich
LOHCs on catalyst surface, the release of H onto catalyst surface, the
recombination of H ions into H2 molecules and the desorption of H2 molecules
and hydrogen-lean LOHCs. Therefore, understanding of the full picture of the
dehydrogenation process is necessary for identifying the critical steps, which
is important for designing cost-efficient, high-performance catalysts for
LOHCs. In this work we compare the adsorption/desorption energetics as well as
kinetics of both NEC-based molecules and hydrogen species on Ru(0001) in an
attempt to identify the energetically critical steps for dehydrogenating
12H-NEC. We also examine the impact of imperfect surfaces on the adsorption of
the molecules based on an amorphous Ru structure. In addition, a critical
challenge for controlling dehydrogenation is the potential fluctuation of
temperature in a big reactor, which may cause the decomposition of LOHCs
gleichweit_dehydrogenation_2013 , especially when the LOHCs contain weak
bonds. In this work we also address the dealkylation of NEC and compare with
the other events in energetics.

Table 1: Calculated dehydrogenation energy (kJ/mol-H2) based on PBE and hybrid
HSE03 functionals, with and without van der Waals (VDW) and thermal
corrections. Thermal corrections and experimental data are for the standard
condition (1 atm and 298 K). More details are in work tang_natural_2020  and
references therein. $\delta E=\Delta E_{\rm{PBE}}-\Delta
E_{\rm{HSE03}}^{\rm{thermal}}$. We note $\Delta E_{\rm{PBE}}$ for perhydro-NEC
in this table (63.9 kJ/mol-H2) is slightly different from that (401.48/6=66.9
kJ/mol-H2) in Table 2. The difference may arise from different parameters used
in the previous work tang_natural_2020  and this study, such as the perhydro-
NEC isomer chosen and the planewave cutoff energy.

Dehydrogenation energy averaged over reaction steps
---
Hydrogen-rich system | $\Delta H^{\circ}_{\rm{expt}}$ | $\Delta E_{\rm{PBE}}$ | $\Delta E_{\rm{HSE}}$ | $\Delta E_{\rm{HSE+VDW}}$ | $\Delta E_{\rm{PBE}}^{\rm{thermal}}$ | $\Delta E_{\rm{HSE}}^{\rm{thermal}}$ | $\Delta E_{\rm{HSE+VDW}}^{\rm{thermal}}$ | $\delta E$
Cyclohexane | 68.7 | 79.36 | 92.29 | 96.26 | 54.62 | 67.56 | 71.52 | 11.8
Decalin | 66.55 | 84.08 | 97.88 | 102.52 | 57.13 | 70.93 | 75.57 | 13.2
Perhydro-NEC | 53.2 | 63.94 | 77.82 | 83.13 | 38.14 | 52.02 | 57.33 | 11.9
Methylcyclohexane | 68.3 | 87.28 | 100.85 | 105.25 | 61.25 | 74.82 | 79.22 | 12.5
Perhydro-trisphaeridine |  | 66.40 | 81.11 |  | 39.83 | 54.54 |  | 12.8

Stepwise dehydrogenation energy for perhydro-trisphaeridine
---
Energy |  step 1 |  step 2 |  step 3 |  step 4 |  step 5 |  step 6 |  step 7
$\Delta E_{\rm{PBE}}$ |  60.94 |  106.98 |  -7.34 |  103.52 |  42.54 |  129.10 |  29.08
$\Delta E_{\rm{HSE}}^{\rm{thermal}}$ |  47.77 |  96.24 |  -21.41 |  94.67 |  31.12 |  117.14 |  16.25
$\delta E$ |  13.17 |  10.74 |  14.07 |  8.85 |  11.42 |  11.96 |  12.83

Figure 1: Comparison of the accuracy of various approximation methods used for
dehydrogenation energy calculation. (top) Calculated dehydrogenation energy
(kJ/mol-H2) of various molecules with different levels of approximation, as
compared with experimental values. The D3 method with Becke-Jonson damping was
used for van der Waals (VDW) correction. (bottom) Stepwise dehydrogenation
energy for perhydro-trisphaeridine. Data from Table 1.

##  II Methods

Geometry optimization of all the relevant structures in this work was
performed using the code VASP kresse_efficient_1996  with the GGA-PBE
exchange-correlation functional and a plane wave basis set with an energy
cutoff of 300 eV. For molecule-only structures, a simulation cell separating
the molecules from their periodic images by more than 10 Å was used, and the
$k$-point sampling was performed at the $\Gamma$ point. For bulk Ru, an
orthogonal cell ($\sim$2.7$\times$4.69$\times$4.28 Å3) with
11$\times$7$\times$7 $k$-point sampling was used for structural optimization.
To represent the catalyst surface, we built an orthogonal Ru(0001)
(7$\times$8) slab out of the optimized bulk Ru structure with 3 atomic layers
in thickness, in view of the balance of accuracy and computational cost. The
molecule/slab structure was then contained in a simulation cell of
$\sim$18.8$\times$18.9$\times$25.5 Å3, which separates the structure with its
periodic images along the surface normal by $\sim$15 Å and is coupled with
$k$-point sampling at the $\Gamma$ point. During the calculations the bottom
Ru layer was fixed. The geometry optimization of the system was stopped after
the force on atoms is below 0.01 eV/Å.

To examine the effect of an imperfect Ru surface on adsorption, we first
constructed an amorphous bulk Ru based on classical molecular dynamics, with
the atomic interactions described by an embedded atom method (EAM) potential
fortini_asperity_2008 . To this end, a bulk Ru of 224 atoms (similar to the
slab structure mentioned above but with 4 layers) was well liquefied at 2500 K
and then quenched to 100 K at 50 K/ps. The simulations were carried out under
an NPT (constant particle number, pressure, and temperature) ensemble using
code LAMMPS plimpton_fast_1995  with Nos$\acute{\rm{e}}$-Hoover thermostat and
barostat, and the timestep was set as 2 fs. The obtained amorphous structure,
after being optimized using DFT, was used to create a slab by inserting a
vacuum layer of $\sim$15 Å into the simulation cell, with the final cell size
being $\sim$19.4$\times$18.7$\times$25.5 Å3. The $n$H-NEC molecules were then
put to 9 different surface sites of the slab for geometry optimization, with
the bottom Ru layer of thickness $\sim$3.5 Å being fixed.

A previous study on natural LOHCs tang_natural_2020  indicates that DFT
calculations using a PBE-based hybrid functional (HSE03) plus thermal
corrections, which include the zero point energy and the enthalpy of the
molecules, can predict the dehydrogenation enthalpy reasonably well. However,
the thermal corrections result in extra computational cost and, more
importantly, sometimes we experienced difficulties for the hybrid functional
calculations to converge, especially for metals. To examine the feasibility of
studying the dehydrogenation reactions using the PBE functional alone, we
tested the performace of PBE and HSE03, with or without thermal and van der
Waals corrections, on several example molecules following the parameters in
reference tang_natural_2020 . As shown in Table 1 and Fig. 1, although the PBE
functional alone overestimates the average and stepwise dehydrogenation
enthalpies by $>$10 and $\sim$9$-$14 kJ/mol-H2, respectively, it provides very
good estimates for the relative energetics. Hence, in this work we address the
dehydrogenation energetics only based on the PBE functional at zero
temperature.

Table 2: System energy (kJ/mol) calculated with the PBE functional at each step ($n$) of the four reactions shown in Fig. 2. For reaction 1, we also computed the energy with the wB97XD functional and cc-pVTZ basis set for comparison, implemented with code Gaussian. $E_{d}$ represents the desorption energy of $n$H-NEC. reaction |  $n=12$ |  $n=10$ |  $n=8$ |  $n=6$ |  $n=4$ |  $n=2$ |  $n=0$
---|---|---|---|---|---|---|---
1 |  0 |  97.55 |  120.22 |  248.55 |  258.59 |  392.61 |  401.48
1 (wB97XD) |  0 |  122.18 |  157.39 |  309.08 |  332.70 |  488.99 |  509.77
2 |  -66.45 |  51.07 |  76.55 |  197.93 |  91.11 |  225.52 |  146.50
3 |  -66.45 |  -5.76 |  -37.12 |  27.43 |  -136.21 |  -58.63 |  -254.98
4 |  -66.45 |  -71.56 |  -168.92 |  -170.27 |  -399.81 |  -388.14 |  -589.89
$E_{d}^{n\rm{H-NEC}}$ |  66.45 |  46.48 |  43.67 |  50.62 |  167.48 |  167.09 |  254.98

Figure 2: Dehydrogenation of 12H-NEC on Ru(0001) surface. (top left) Schematic
illustration of four reaction paths related to the dehydrogenation. (top
right) System energy along the four reactions as a function of $n$, the number
of H atoms charged to NEC, with the energies referenced to state A. For
reaction 1, we also computed the energy with the wB97XD functional using code
Gaussian for comparison. The energy values are in Table 2. (bottom) Structures
of 12-NEC isomers (a-f) with energy referenced to the most stable isomer and
the structures of $n$H-NEC considered during dehydrogenation.

For the dehydrogenation process, we assume the initial state (denoted as state
A) as a 12H-NEC molecule far from the Ru surface and the final state (state F)
as a NEC molecule plus 6 H2 molecules far from the Ru surface, as
schematically shown in Fig. 2. For energetics study, we consider the following
intermediate states: 12H-NEC being adsorbed onto the Ru surface (state B), all
of the 12 H ions being released (state C), and the H ions combining into 6 H2
molecules on the surface (state D), the H2 molecules being desorbed (state E),
and the system reaches state F upon NEC desorption. Alternative to the process
C→D→E, direct desorption of H ions into gaseous H2 (C→E) is also possible, of
which the energetics is straightforward. We note that the above order of
states is assumed for convenience (for example, experimentally H2 desorption
may start while some remaining H ions are not released yet, or NEC desorption
may occur before H2 desorption is finished) and it does not affect the
conclusions of this work. To identify the effects of various reaction
processes and the catalyst on the dehydrogenation energetics, we computed the
system energies along four reaction paths as detailed below.

Figure 3: Electronic structures of 8H-NEC and 4H-NEC on Ru(0001). (left)
Structure and isosurface of the differential charge density, or the charge
density minus the superposition of atomic charge densities. The isosurfaces
represent charge density level of 10 $\mu e/$Å3 (red) and $-10~{}\mu e/$Å3
(blue). The brown and green spheres represent different Ru layers, and C, N,
and H atoms are shown in grey, white, and blue, respectively. (right)
Projected density of states of H and C atoms as well as three surface Ru atoms
(labelled as 1-3 in the left).

Reaction 1 (state A to state F) was considered for a reference scenario where
Ru is not involved in the dehydrogenation process. In this case we treated the
molecules and Ru slab separately and at each reaction step the system energy
is the sum of the separate parts, which can be written as

$E(n)=E_{n{\rm{H\mathchar 45\relax
NEC}}}+\frac{12-n}{2}E_{\rm{H}_{2}}+E_{\rm{Ru}}$ (1)

where $n$ decreases from 12 for the initial state A to 0 for the final state
F.

For reaction 2, we considered the transition from state B to state E, with the
released H ions directly forming gas H2 molecules far from Ru surface. This
hypothetical case is similar to reaction 1 except the $n$H-NEC molecule is
adsorbed on Ru surface, and so it allows for identifying the effect of
$n$H-NEC adsorption alone on dehydrogenation energetics. In this case, the
system energy is

$E(n)=E_{n{\rm{H\mathchar 45\relax NEC/Ru}}}+\frac{12-n}{2}E_{\rm{H}_{2}}$ (2)

Similarly, we considered transitions from B to D (reaction 3) and from B to C
(reaction 4), with the released H ions forming adsorbed H2 molecules and
adsorbed H ions, respectively. This allows for identifying the adsorption
effects of H2 and H ions. Here we treat H2 (or H ions) and the $n$H-NEC
separately, and hence for reaction 3, the system energy can be written as

$E(n)=E_{n{\rm{H\mathchar 45\relax
NEC/Ru}}}+\frac{12-n}{2}\mu_{{}_{\rm{H_{2}}}}$ (3)

where the chemical potential
$\mu_{{}_{\rm{H_{2}}}}$=$(E_{\rm{H}_{2}/Ru}$$-$$E_{\rm{Ru}})$ represents the
energy of H2 adsorbed on Ru surface. $E(n)$ for reaction 4 can be similarly
defined.

The structures for $n$H-NEC we considered at various steps are shown at the
bottom of Fig. 2, which are suggested by a few previous studies fei_study_2017
; dong_dehydrogenation_2016 ; yu_bimetallic_2020 . For 12H-NEC, there are six
possible isomerseblagon_study_2010 , labelled as a-f in Fig. 2, and we chose
the one with the lowest calculated energy at 0 K. The chosen isomer is also
one of the experimentally observed species eblagon_study_2010 .

Figure 4: Effect of ethyl group orientation on the adsorption of 2H-NEC on
Ru(0001). The preferred ethyl-up state enables a better match of the aromatic
ring to surface Ru.

##  III Results and discussions

###  III.1 Dehydrogenation thermodynamics

The computed system energies for the four reactions, all referenced to that of
state A, are listed in Table 2 and shown in Fig. 2(right). For reaction 1, the
calculated system energy increases as dehydrogenation proceeds, with the
average energy increase to be $\sim$67 kJ/mol-H2. For reaction 2, the system
energies are reduced because of the adsorption of $n$H-NEC. The desorption
energy ($E_{d}$) for $n$H-NEC, which is the $E$ difference between reactions 1
and 2, ranges from $\sim$44 to $\sim$66 kJ/mol for $n$$\geq$6 and from
$\sim$167 to $\sim$255 kJ/mol for $n$$\leq$4\. We attribute the significantly
higher $E_{d}$ for $n$$\leq$4 to the formation of the aromatic ring, of which
the flat structure allows for stronger C-Ru bonding with the flat Ru(0001)
surface. To illustrate this we compare the the structures of 8H-NEC and 4H-NEC
as an example. The former molecule adsorbs to Ru mainly via H atoms and hence
the carbon rings are not as well aligned with the surface Ru atoms as the
aromatic carbon ring in 4H-NEC (top view in Fig. 3). The charge density
isosurfaces in Fig. 3 show that more electrons transfer from surface Ru
towards 4H-NEC than towards 8H-NEC because of the Ru-C bonding. To
quantitatively illustrate this, we performed Bader charge analysis. As shown
in Table 3, the three Ru atoms beneath the aromatic ring of 4H-NEC lose
$\sim$0.6 more electrons as compared to their counterparts beneath 8H-NEC.
Fig. 3 also compares their projected density of states (PDOS) for H and C
atoms as well as the three surface Ru atoms beneath the aromatic ring of
4H-NEC and the corresponding ring of 8H-NEC. By integrating over the PDOS, one
can compute the electron energy of the corresponding orbitals (per atom)
according to

$E_{e}=\frac{1}{N_{a}}\int^{E_{F}}_{-\infty}n(\varepsilon)\varepsilon
d\varepsilon$ (4)

where $n(\varepsilon)$ is the PDOS density, $E_{F}$ is the Fermi energy, and
$N_{a}$ is the number of H, C, or Ru atoms contributing to the PDOS. We found
$E_{e}$ values of $s_{{}_{\rm{H}}}$, $p_{{}_{\rm{C}}}$ and $d_{\rm{Ru}}$ of
4H-NEC/Ru are lower than those of 8H-NEC/Ru by $\sim$0.3, $\sim$1.2, and
$\sim$1.3 eV, respectively. This indicates the stronger bonding between 4H-NEC
and Ru, consistent with the charge density analysis.

Table 3: Bader charges of selected surface Ru atoms (refer to Fig. 3) and average Bader charges of the top surface Ru atoms. The charges, in unit $e$, are referenced to neutral Ru atoms. The analysis show that the top Ru layer gain electrons from the beneath Ru layer, resulting in a positive average charge. System | atom 1 | atom 2 | atom 3 | Surface Ru
---|---|---|---|---
8H-NEC/Ru | +0.11 | 0.00 | +0.04 | +0.07
4H-NEC/Ru |  $-$0.11 |  $-$0.18 |  $-$0.15 | +0.06
Figure 5: Adsorption of H2 molecule and 2 H atoms on Ru(0001) surface. In (a)
and (b), the red squares schematically represent metastable adsorption sites.
In (a), the adsorption energies
($E_{\rm{H_{2}/Ru}}-E_{\rm{Ru}}-E_{\rm{H_{2}}}$) were calculated to be
$-$56.82, $-$2.47, and +0.78 kJ/mol for atop${}_{{}_{\rm{F}}}$,
atop${}_{{}_{\rm{V}}}$, and hcp${}_{{}_{\rm{V}}}$ site, respectively. The
subscript F or V indicates a flat or vertical position of the H2. In (b), the
adsorption energies ($E_{\rm{2H/Ru}}-E_{\rm{Ru}}-2$$\times$$E_{\rm{H}}$) are
$-$548.73, $-$528.06, and $-$557.05 kJ/mol for hcp-hcp, hcp-fcc, and fcc-fcc
site, respectively, which are close to previous calculations
chou_theoretical_1989 . The adsorption energies become $-$114.36, $-$93.69,
and $-$122.69 kJ/mol, respectively, if a H2 molecule is assumed in the initial
state (i.e., 2$\times$$E_{\rm{H}}$ is replaced by $E_{\rm{H_{2}}}$ in the
formula). In (c) for the unstable sites, the squares/arrows represent the
initial/final positions of hydrogen.

For $n$H-NEC adsorption we tested the effect of the ethyl group state
(pointing down to the surface or pointing up) by rotating the ethyl group
around the neighboring N-C bond. We found that for $n\leq 4$ the up state is
highly favorable (by energy difference more than 70 kJ/mol) while for the
other $n$H-NEC structures the two states are nearly equally favorable (with
energy difference less than 10 kJ/mol). For 4H- and 2H-NEC, pointing-down of
the ethyl group reduces the bonding of the molecule with Ru by displacing the
aromatic ring from the favorable position (Fig. 4). On the other hand, for
$n$=0, the pointing-down ethyl group relaxes into a flat position, allowing
the aromatic rings to stay parallel to the Ru surface, but we found that a
pointing-up position is energetically more favorable. For $n$$\geq$6, the
charged H atoms increase the space between the carbon rings and the Ru
surface, making the up/down state of the ethyl group relatively insignificant.
For gaseous $n$H-NEC molecules we did not find significant energy differences
between the ethyl group states. We note that our results are different from a
previous study eblagon_study_2010  which reported negligible adsorption energy
for 8H-NEC with the ethyl group facing down and $\sim$30 kJ/mol for facing up.
However, both studies predict that the adsorption of 8H-NEC on a flat Ru
surface is relatively weak, which could result in the desorption of 8H-NEC
molecule. It was found that 8H-NEC prefers low-coordination sites, such as
surface steps, to a flat surface eblagon_study_2010 .

Figure 6: Activation energies for dehydrogenation. (a) Hydrogen discharge from
$n$H-NEC, with the reaction coordinate indicating the initial (0),
intermediate, and final (1) states, and for each discharge event, the system
energy is referenced to the initial state. The structure pictures are for the
initial/final state of each discharge event. For 2H-NEC to 1H-NEC, the picture
also shows further structural relaxation after the discharged H diffuses away.
(b) Activation energies for other relevant surface events.

If the released hydrogen are adsorbed onto the Ru surface, their chemical
potential is lower than that of gas H2 molecule and hence reduces the system
energy. We studied the adsorption of hydrogen molecules and ions on various Ru
surface sites (Fig. 5). It was found that hydrogen molecules are stable on the
atop sites while two hydrogen ions could stay on hcp sites, fcc sites, or
their combination, with the fcc sites being favorable. Assuming hydrogen
chemical potential based on the most favorable hydrogen molecule/ion
adsorptions, we obtained the system energies for reactions 3 and 4,
respectively, as shown in Fig. 2. These two reactions reveal an exothermic
hydrogen release process, with the path of hydrogen ions being more favorable.

Fig. 2 reveals not only the most probable path for 12-NEC dehydrogenation but
also the contributions of different factors ($n$H-NEC desorption, H2
desorption, and H ion recombination) to dehydrogenation energetics, which
allows for identifying the energetically critical step. For example, for the
dehydrogenation path
A$\rightarrow$B$\rightarrow$C$\rightarrow$D$\rightarrow$E$\rightarrow$F, the
first two steps are exothermic and the final three steps cost 335, 401, and
255 kJ/mol, respectively. In view that C$\rightarrow$D and D$\rightarrow$E are
for 6 H2 molecules while E$\rightarrow$F is the desorption of a single NEC
molecule, it is clear that NEC desorption is the energetically critical step
for the dehydrogenation. Because of the reversibility of (de)hydrogenation, it
is easy to see that, for NEC hydrogenation (assuming the same catalyst), the
endothermic step C$\rightarrow$B costs 523 kJ/mol, or on average $\sim$44
kJ/mol for charging each H atom, and the critical desorption of 12H-NEC costs
$\sim$66 kJ/mol. As the lowest energy point in Fig. 2, state C is also
critical since too strong adsorption of H ions on the catalyst surface could
be an energy trap that slows the reactions down. From Fig. 2 one also can
predict the kinetically stable intermediate products, 8H-NEC and 4H-NEC, as
observed in experiments dong_dehydrogenation_2016 ; jiang_experimental_2019 .
Along path B$\rightarrow$C, the system energies of $n$=8 and $n$=4 are close
to those of $n$=6 and $n$=2, respectively, which makes 8H-NEC and 4H-NEC
kinetically stable in view of the activation energy for the reactions. On the
other hand, the big thermodynamic driving force makes the other intermediate
products ($n$=10, 6, 2) much less stable. Especially, we note that the
formation of aromatic sextets (transitions $n$=6→$n$=4 and $n$=2→$n$=0) is
associated with a big driving force ($\sim$200 kJ/mol). Similar phenomena were
observed in our previous studies for bio-based LOHCs tang_natural_2020 .

Due to the strong NEC-Ru bonding and the relatively weak nitrogen-alkyl
bonding, it is possible for dealkylation to occur before the desorption of
NEC. Here we considered the dealkylation reaction by comparing the energy of
NEC/Ru with two additional H atoms adsorbed on two fcc sites and that of a
carbazole and an ethane adsorbed on Ru. It was found the latter is higher in
energy by $\sim$98.5 kJ/mol, much smaller than the desorption energy of NEC,
implying the possibility of dealkylation under suitable situations. Further,
we note that the desorption of the resulting carbazole and ethane costs
$\sim$212 kJ/mol, while the desorption energies of NEC and 2H (from two fcc
sites to gaseous H2) sum to $\sim$255+123=378 kJ/mol. This indicates that the
dealkylation approach is thermodynamically favorable by
378$-$212$-$98.5$\approx$67 kJ/mol, which essentially is energy difference
between (NEC+H2) and (carbazole+ethane) in their gaseous state.
Experimentally, dealkylation of NEC on Pt(111) surface
gleichweit_dehydrogenation_2013  was reported for temperatures above 390 K.

###  III.2 Dehydrogenation kinetics

The above discussions reveal the thermodynamics during the (de)hydrogenation
process of NEC on Ru. We also studied the kinetics of the process by computing
the energy barrier of relevant steps on the Ru surface using the climbing
image nudged elastic band method henkelman_climbing_2000 . The steps
considered here include H discharge from $n$H-NEC, H diffusion, the
combination of 2 H into a H2 molecule, and the desorption of H2 and NEC. For
hydrogen discharge, although it is proper to consider H pairs in the above
thermodynamics calculations and in theory it is possible for H to be
discharged in pairs, here for the kinetics we consider a single H atom,
consistent with experimental observations crawford_understanding_2007 . The H
atoms on $n$H-NEC can be classified into three types, with the C-H bond
roughly pointing down towards the Ru surface, parallel to the surface, or
pointing up. We computed the activation energy for discharging a H atom from
12H-NEC, 4H-NEC, and 2H-NEC for example, with the corresponding C-H bond
either pointing down or relatively parallel to the surface. As shown in Fig.
6(a), the activation energies for these cases range from $\sim$50 to $\sim$100
kJ/mol, which are similar to the calculated barriers ($\sim$39 to $\sim$92
kJ/mol) for dehydrogenating tetrahydrocarbazole to carbazole on palladium
surface crawford_understanding_2007 . We note that for
4H-NEC$\rightarrow$3H-NEC and 2H-NEC$\rightarrow$1H-NEC, the final state is
higher in energy than the initial state (Fig. 6(a)), partially because the
relaxation of the system is limited when the discharged H atom sits below the
$n$H-NEC. For 2H-NEC$\rightarrow$1H-NEC as an example, we found that once the
discharged H atom diffuses away, the system decreases in energy by $\sim$130
kJ/mol with stronger bonding to the Ru surface.

For most 12H-NEC isomers (see Fig. 2, bottom) adsorbed on the Ru surface,
there exist some C-H bonds of the five-membered ring pointing up. The
activation energy for discharging an upward H atom to Ru surface was found to
be higher than 300 kJ/mol. Alternatively, we considered rotating the upward
C-H bonds downwards, which allows easier H discharge later, and found similar
activation energy. In view of the relatively low desorption energy of $n$H-NEC
for $n$$\geq$6, a more probable mechanism for discharging these upward H atoms
is that the $n$H-NEC molecule flips over the Ru surface and hence turns them
downwards. In the presence of a solvent that accepts hydrogen, it is also
possible to discharge the upward H atoms via the solvent
crawford_understanding_2007 .

We also computed the activation energies for other relevant events on Ru
surface. For hydrogen migrating from a fcc site to a neighboring hcp site, the
activation energies was found to be $\sim$14.9 kJ/mol only, which is
consistent with the high mobility of H on Ru surface. For the combination of
two neighboring fcc-site H atoms into an atop H2 molecule, $\sim$68.9 kJ/mol
is necessary to activate the process. Finally, we examined the desorption of
H2 and NEC. For these events, we set the final states to be local energy
minima corresponding to the H2 and NEC molecules being $\sim$4 and $\sim$6 Å,
respectively, above the Ru surface. We note that these final states are
slightly lower in energy than the states with the molecules infinitely away
from the surface (by $\sim$4.2 and $\sim$5.5 kJ/mol for H2 and NEC,
respectively). Based on the selected final states, we computed the activation
energies for H2 and NEC desorption to be 54.3 and 249.5 kJ/mol, respectively,
which are close to the desorption energies.

In the above discussion of H2 desorption, we assumed the path of two H atoms
recombining into an atop H2 before the desorption. The intermediate atop H2
tends to decompose back into H atoms due to the small activation energy.
Hence, we also studied the direct desorption of two H atoms from neighboring
fcc sites into a gaseous H2 molecule, as shown in Fig. 6(b). We found the
activation energy to be $\sim$122 kJ/mol, close to the desorption energy. The
calculated desorption energy in this work compares well to the experimental
value, $\sim$120 kJ/mol, reported for low H coverage on Ru(0001) surface
feulner_adsorption_1985 .

Overall, the kinetics calculations also identify the desorption of NEC as the
critical step for dehydrogenation. On the other hand, the discharged H ions
are found to be very mobile on Ru surface.

###  III.3 Effect of surface morphology

In the above we have discussed the dehydrogenation thermodynamics and kinetics
of perhydro-$N$-ethylcarbazole on perfect Ru(0001) surface. In practice, a
crystal surface may have structural defects and these defects may impact on
the behavior of adsorbed molecules. For example, a previous study
eblagon_study_2010  has shown that surface steps enables more stable
adsorption of 8H-NEC as compared to a flat surface. Here, in a more general
way, we use an amorphous Ru surface to represent an imperfect surface in view
of its rich local structural features.

Fig. 7 shows the average system energy for reaction 2 on the amorphous
surface, as compared to that on Ru(0001). It can be seen that the adsorption
energies of 12H- and 10H-NEC are similar for the two cases, but the energies
for 8H- and 6H-NEC on the amorphous surface are about 50 kJ/mol lower than on
Ru(0001). This highlights the impact of surface morphology and is
qualitatively consistent with the stabilization effect of a surface step on
8H-NEC eblagon_study_2010 . The adsorption of 4H- and 2H-NEC may prefer a
relatively flat surface due to the formation of the flat aromatic ring, and
consequently the system energy for the amorphous case is slightly higher. This
effect is more significant for the NEC case, with its adsorption energy on the
amorphous surface is about 100 kJ/mol higher than on Ru(0001). Overall,
however, the energy trends for reaction 2 in the two cases are similar, and
the energy difference observed here does not alter the conclusions based on
the crystalline phase.

Figure 7: System energy for reaction 2 on amorphous Ru surface and Ru(0001),
with the energies referenced to state A for both cases. The energy for the
amorphous case is averaged over nine different adsorption sites.

##  IV Summary

Based on first principles calculations we studied the thermodynamics and
kinetics of the complete dehydrogenation path of perhydro-$N$-ethylcarbazole
(12H-NEC) on Ru(0001) surface, involving the adsorption of 12H-NEC, the
release of H ions, H ion diffusion and recombination into H2, and the
desorption of H2 and hydrogen-lean NEC. During the dehydrogenation process,
the adsorption of $n$H-NEC on Ru(0001) is significantly strengthened upon the
formation of aromatic ring, of which the flat structure allows for stronger
C-Ru bonding. Although the whole dehydrogenation process is endothermic, the
release of H from $n$H-NEC was found to be exothermic because of H adsorption
onto Ru surface. The desorption of flat, hydrogen-lean NEC, which costs
$\sim$255 kJ/mol, was identified as the most energy demanding step. Based on
an amorphous model, we also showed that the surface morphology has a great
impact on the stability of molecule adsorption. Overall, the results imply
more efficient dehydrogenation could be achieved by weakening the bonding of
NEC to catalysts, either through engineering catalyst surface (such as surface
defects or smaller catalyst particles) or different catalyst materials. Our
calculations also revealed possible dealkylation at elevated temperatures.

Acknowledge CT thanks the financial support from the Australian National
University Grand Challenge program (Zero-Carbon Energy for the Asia-Pacific).
This work was supported by computational resources provided by the Australian
Government through the National Computational Infrastructure (NCI) under the
ANU Merit Allocation Scheme.

Data Availability The raw/processed data required to reproduce these findings
cannot be shared at this time as the data also forms part of an ongoing study.

## References

  * (1) Sotoodeh F and Smith KJ, An overview of the kinetics and catalysis of hydrogen storage on organic liquids, Can. J. Chem. Eng., 91, 1477, (2013).
  * (2) Preuster P, Papp C, and Wasserscheid P, Liquid Organic Hydrogen Carriers (LOHCs): Toward a Hydrogen -free Hydrogen Economy, Accounts Chem. Res., 50, 74, (2017).
  * (3) Modisha PM, Ouma CNM, Garidzirai R, Wasserscheid P, and Bessarabov D, The Prospect of Hydrogen Storage Using Liquid Organic Hydrogen Carriers, Energy Fuels, 33, 2778, (2019).
  * (4) Abdin Z, Tang C, Liu Y, and Catchpole K, Large-scale stationary hydrogen storage via liquid organic hydrogen carriers, iScience, 24, 102966, (2021).
  * (5) Niermann M, Drünert S, Kaltschmitt M, and Bonhoff K, Liquid organic hydrogen carriers (LOHCs) – techno-economic analysis of LOHCs in a defined process chain, Energy Environ. Sci., 12, 290, (2019).
  * (6) Tang C, Fei S, Lin GD, and Liu Y, Natural liquid organic hydrogen carrier with low dehydrogenation energy: A first principles study, International Journal of Hydrogen Energy, 45, 32089, (2020).
  * (7) Biniwale RB, Kariya N, and Ichikawa M, Dehydrogenation of cyclohexane over Ni based catalysts supported on activated carbon using spray-pulsed reactor and enhancement in activity by addition of a small amount of Pt, Catal. Lett., 105, 83, (2005).
  * (8) Xia Z, Lu H, Liu H, Zhang Z, and Chen Y, Cyclohexane dehydrogenation over Ni-Cu/SiO2 catalyst: Effect of copper addition, Catalysis Communications, 90, 39, (2017).
  * (9) Pande JV, Shukla A, and Biniwale RB, Catalytic dehydrogenation of cyclohexane over Ag-M/ACC catalysts for hydrogen supply, International Journal of Hydrogen Energy, 37, 6756, (2012).
  * (10) Kariya N, Fukuoka A, and Ichikawa M, Efficient evolution of hydrogen from liquid cycloalkanes over Pt-containing catalysts supported on active carbons under ”wet-dry multiphase conditions”, Appl. Catal. A-Gen., 233, 91, (2002).
  * (11) Kariya N $et~{}al$., Efficient hydrogen production using cyclohexane and decalin by pulse-spray mode reactor with Pt catalysts, Appl. Catal. A-Gen., 247, 247, (2003).
  * (12) Shin BS, Yoon CW, Kwak SK, and Kang JW, Thermodynamic assessment of carbazole-based organic polycyclic compounds for hydrogen storage applications via a computational approach, Int. J. Hydrog. Energy, 43, 12158, (2018).
  * (13) Yook H, Kim K, Park JH, Suh Y, and Han JW, Density functional theory study on the dehydrogenation of 1,2-dimethyl cyclohexane and 2-methyl piperidine on Pd and Pt catalysts, Catal. Today, 352, 345, (2020).
  * (14) Ouma CNM, Modisha P, and Bessarabov D, Insight into the adsorption of a liquid organic hydrogen carrier, perhydro-i-dibenzyltoluene (i = m, o, p), on Pt, Pd and PtPd planar surfaces, RSC Adv., 8, 31895, (2018).
  * (15) Crawford P $et~{}al$., Understanding the dehydrogenation mechanism of tetrahydrocarbazole over palladium using a combined experimental and density functional theory approach, J. Phys. Chem. C, 111, 6434, (2007).
  * (16) Sotoodeh F, Huber BJM, and Smith KJ, Dehydrogenation kinetics and catalysis of organic heteroaromatics for hydrogen storage, Int. J. Hydrog. Energy, 37, 2715, (2012).
  * (17) Eblagon KM $et~{}al$., Study of Catalytic Sites on Ruthenium For Hydrogenation of N-ethylcarbazole: Implications of Hydrogen Storage via Reversible Catalytic Hydrogenation, J. Phys. Chem. C, 114, 9720, (2010).
  * (18) Gleichweit C $et~{}al$., Dehydrogenation of Dodecahydro-N-ethylcarbazole on Pt(111), ChemSusChem, 6, 974, (2013).
  * (19) Kresse G and Furthmuller J, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B, 54, 11169, (1996).
  * (20) Fortini A, Mendelev MI, Buldyrev S, and Srolovitz D, Asperity contacts at the nanoscale: Comparison of Ru and Au, Journal of Applied Physics, 104, 074320, (2008).
  * (21) Plimpton S, Fast Parallel Algorithms for Short-Range Molecular-Dynamics, J. Comput. Phys., 117, 1, (1995).
  * (22) Fei S $et~{}al$., A study on the catalytic hydrogenation of N-ethylcarbazole on the mesoporous Pd/MoO3 catalyst, Int. J. Hydrog. Energy, 42, 25942, (2017).
  * (23) Dong Y, Yang M, Mei P, Li C, and Li L, Dehydrogenation kinetics study of perhydro-N-ethylcarbazole over a supported Pd catalyst for hydrogen storage application, Int. J. Hydrog. Energy, 41, 8498, (2016).
  * (24) Yu H $et~{}al$., Bimetallic Ru-Ni/TiO2 catalysts for hydrogenation of N-ethylcarbazole: Role of TiO2 crystal structure, Journal of Energy Chemistry, 40, 188, (2020).
  * (25) Chou MY and Chelikowsky JR, Theoretical study of hydrogen adsorption on Ru(0001): Possible surface and subsurface occupation sites, Phys. Rev. B, 39, 5623, (1989).
  * (26) Jiang Z, Gong X, Wang B, Wu Z, and Fang T, A experimental study on the dehydrogenation performance of dodecahydro-N-ethylcarbazole on M/TiO2 catalysts, International Journal of Hydrogen Energy, 44, 2951, (2019).
  * (27) Henkelman G, Uberuaga BP, and Jonsson H, A climbing image nudged elastic band method for finding saddle points and minimum energy paths, J. Chem. Phys., 113, 9901, (2000).
  * (28) Feulner P and Menzel D, The adsorption of hydrogen on ruthenium (001): Adsorption states, dipole moments and kinetics of adsorption and desorption, Surface Science, 154, 465, (1985).