Fatih Tas

Berkay Demireller
email: bdemireller@sabanciuniv.edu

Kamer Kaya
email: kaya@sabanciuniv.edu

Bora Uc ¸ar
email: bora.ucar@ens-lyon.fr

Streaming Hypergraph Partitioning Algorithms on Limited Memory Environments

Keywords: Hypergraph partitioning, streaming hypergraphs, single-board computers

Many well-known, real-world problems involve dynamic, interrelated data items. Hypergraphs are powerful combinatorial structures that are frequently used to model such data. Many of today's data-centric have streaming data; new items arrive continuously, and the data grow over time. With paradigms such as Internet of Things and Edge Computing, such applications become more natural and more practical. In this work, we assume a streaming model where the data items and their relations are modeled as a hypergraph, which is generated at the edge. This hypergraph is then partitioned, and the parts are sent to remote nodes via an algorithm running on a memory-restricted device, such as a single board computer. Such a partitioning is usually performed by taking a connectivity metric into account to minimize the communication cost of later analyses that will be performed in a distributed fashion. Although there are many offline tools that can partition static hypergraphs effectively, algorithms for the streaming settings are rare. We analyze a well-known algorithm from the literature and significantly improve its run time by altering its inner data structure. On a medium-scale hypergraph, the new algorithm reduces the run time from 17800 seconds to 10 seconds. We then propose sketch-and hash-based algorithms, as well as ones that can leverage extra memory to store a small portion of the data to enable the refinement of partitioning when possible. We experimentally analyze the performance of these algorithms and report their run times, connectivity metric scores, and memory uses on a high-end server and four different single-board computer architectures.

I. INTRODUCTION

Real-world data can be complex, multi-model and multidimensional with irregular relations among the data items. Most of the models such as column-or row-oriented tabular representation fail in capturing the essence of knowledge contained in such data. Hypergraphs, which are generalizations of graphs, are highly flexible and appropriate for modeling and analysing such data. Therefore, they are used in various areas: DNA sequencing [START_REF] Venkatraman | Unimodular hypergraph for DNA sequencing: A polynomial time algorithm[END_REF], scientific computing [START_REF] Aykanat | Hypergraph-partitioning based decomposition for parallel sparse-matrix vector multiplication[END_REF], VLSI design [START_REF] Karypis | Multilevel hypergraph partitioning: applications in vlsi domain[END_REF], citation recommendation [START_REF] Küc ¸üktunc | Fast recommendation on bibliographic networks[END_REF], finding semantic similarities between documents [START_REF] Menezes | Semantic hypergraphs[END_REF], finding descriptor similarities between images [START_REF] Skiker | The representation of semantic similarities between object concepts in the brain: a hypergraph-based model[END_REF], and classification [START_REF] Sun | Hypergraph spectral learning for multilabel classification[END_REF].

Distributed graph and hypergraph stores became popular in today's applications. A good partitioning of the data among the compute nodes in a distributed framework is necessary to reduce the communication in the upstream applications. With the increasing popularity of data-centric paradigms such as Internet of Things and Edge Computing, the data that is fed to these stores started to have a streaming fashion. In this work, we assume that the data is generated/processed at the edge of a network and partitioned on a memory-restricted device, such as a single-board computer (SBC). There are popular algorithms to partition streaming graphs [START_REF] Zhang | Akin: A streaming graph partitioning algorithm for distributed graph storage systems[END_REF], [START_REF] Tsourakakis | Fennel: Streaming graph partitioning for massive scale graphs[END_REF], [START_REF] Stanton | Streaming graph partitioning for large distributed graphs[END_REF], and two recent benchmarks to evaluate the performance of such algorithms [START_REF] Pacaci | Experimental analysis of streaming algorithms for graph partitioning[END_REF], [START_REF] Abbas | Streaming graph partitioning: An experimental study[END_REF]. Although hypergraphs have a more expressive power capability, and there are fine-tuned, optimized, offline hypergraph partitioning tools, e.g., [START_REF] Andre | Memetic multilevel hypergraph partitioning[END_REF], [START_REF] Aykanat | PaToH (Partitioning Tool for Hypergraphs)[END_REF], [START_REF] Deveci | Hypergraph partitioning for multiple communication cost metrics: Model and methods[END_REF], [START_REF] Karypis | Multilevel hypergraph partitioning: applications in vlsi domain[END_REF], the hypergraph partitioning problem in the streaming setting is not addressed thoroughly.

Hypergraphs generalize graphs. In a graph, the edges represent pairwise connections, whereas in a hypergraph hyperedges (or nets for short) represent multi-way connections. That is, a net connects more than two vertices. In a streaming setting, this difference makes the hypergraph partitioning problem much harder than the graph partitioning problem. For graph partitioning, when a vertex appears with its edges, the endpoint vertex IDs are implicit; hence, just the part information of the vertices is sufficient to judiciously decide on the part of the vertex at hand. However, in a hypergraph, a vertex appears with its nets and the neighbor vertices are not implicit. Hence, one needs to keep track of the connectivity among the nets and the parts to effectively assign the current vertex to a part.

We assume a streaming setting where the vertices of a hypergraph appear in some order along with the id's of their nets. The contribution of the study is three-fold:

• We take one of the existing and popular algorithms from the literature [START_REF] Alistarh | Streaming min-max hypergraph partitioning[END_REF] and make it significantly faster by modifying its inner data structure used to store the part-to-net connectivity.

• We propose techniques to refine the existing partitioning at hand with the help of some extra memory to store some portion of the hypergraph. • We propose and experiment with various algorithms and benchmark their run times, memory requirements, and partitioning quality on a high-end server and multiple SBCs.

The rest of the paper is organized as follows. Section II presents the notation and background on streaming hypergraph partitioning. The proposed algorithms are presented in Section III. The related work is summarized in Section IV.

Section V presents the experiments, and Section VI concludes the paper.

II. NOTATION AND BACKGROUND

A hypergraph H = (V, N ) consists of a set V of vertices and a set N of nets. A net n ∈ N contains a set of vertices, where the vertices in n are called its pins. The size of n ∈ N is the number its pins, and the degree of v ∈ V is the number of nets containing v. The notation pins[n] and nets[v] represent the pin set of a net n, and the set of nets containing a vertex v, respectively. We assume unit weighted vertices and nets.

A K-way partition of a hypergraph H, which is denoted as

Π = {V 1 , V 2 , . . . , V K }, is a vertex partition where • parts are pairwise disjoint, i.e., V k ∩ V = ∅ for all 1 ≤ k < ≤ K, • each part V k is a nonempty subset of V, i.e., V k ⊆ V and V k = ∅ for 1 ≤ k ≤ K, • the union of K parts is equal to V, i.e., K k=1 V k = V.
In the streaming setting, the vertices in V appear one after another. The elements of the stream are pairs (v, nets[v]). For each stream element, v will be partitioned, i.e., the part vector entry part[v] will be set by the partitioning algorithm. In the strict streaming setting, each stream element (v, nets[v]) is forgotten after part[v] is decided. Besides, none of the partitioning decisions can be revoked. In a more flexible streaming setting, a buffer with a capacity B is reserved to store some of the net sets. These vertices can then be reprocessed and re-partitioned. In this setting, the cost of storing

(v, nets[v]) in the buffer is |nets[v]|.
At any time point of the partitioning, the partition must be kept balanced by limiting the difference between the number of vertices in the most loaded and the least loaded parts. Let s be the slack denoting this difference. We say that the partition is balanced if and only if

abs(|V i | -|V j |) ≤ s for all 1 ≤ i < j ≤ K , (1)
where abs(x) is the absolute value of x. For a K-way partition Π, a net that has at least one pin (vertex) in a part is said to connect that part. The number of parts connected by a net n is called its connectivity and denoted as λ n . A net n is said to be uncut (internal) if it connects only one part (i.e., λ n = 1), and cut (external), otherwise (i.e., λ n > 1). Given a partition Π, if a vertex is in the pin set of at least one cut net, it is called a boundary vertex. We use parts[n] to denote the set of parts net n is connected to. Let Λ(n, p) = |pins[n] ∩ V p | be the number of pins of net n in part p. Hence, Λ(n, p) > 0 if and only if p ∈ parts [n]. There are various metrics to measure the quality of a partitioning in terms of the connectivity of the nets [START_REF] Lengauer | Combinatorial Algorithms for Integrated Circuit Layout[END_REF]. The one which is widely used in the literature and shown to accurately model the total communication volume of many data-processing kernels is called the connectivity-1 metric. This cutsize metric is defined as:

χ(Π) = n∈N (λ n -1) . (2)
In this metric, each cut net n contributes (λ n -1) to the cut size. The hypergraph partitioning problem can be defined as the task of finding a balanced partition Π with K parts such that χ(Π) is minimized. This problem is NP-hard even in the offline setting [START_REF] Lengauer | Combinatorial Algorithms for Integrated Circuit Layout[END_REF].

III. PARTITIONING STREAMING HYPERGRAPHS

In describing the algorithms, we assume that the stream elements come from a hypergraph H, where a vertex of H is revealed with its nets. We use the notation |H| to denote the number of pins in this hypergraph while analyzing algorithms.

The simplest partitioning algorithm one can use in the streaming setting is random partitioning, RANDOM, which assigns the vertex in a stream element to a random part, while keeping the partitioning always balanced as shown in Algorithm 1. In the algorithm, p is the candidate part, p min is the part ID having the least number of vertices, and rand(1, K) chooses a random integer in between [1, K]. When the difference between the number of vertices is equal to s, v cannot be assigned to p since this decision makes the partitioning unbalanced.

Algorithm 1: RANDOM Input: K, s Output: part[•] for all (v, nets[v]) in streaming order do p ← rand(1, K) while |V p | -|V pmin | = s do p ← rand(1, K) part[v] ← p if p = p min then Update p min return part A. Min-Max Partitioning
MINMAX is a well-known approach proposed for streaming hypergraph partitioning [START_REF] Alistarh | Streaming min-max hypergraph partitioning[END_REF]. The approach, whose pseudocode is given in Algorithm 2, keeps track of net connectivity, i.e., which net is connected to which part, by keeping a net set p2n[i] for each 1 ≤ i ≤ K. Each streaming vertex v is assigned to the part p with the largest intersection set p2n[p] ∩ nets[v] that does not violate the balance constraint. After setting part[v] to p, there will not be an additional cost for each net in p2n[p] ∩ nets [v]. Hence, the maximum intersection cardinality is a good greedy decision. The downside of this approach is that there is no way of knowing if any of v's nets are connected to a part i without checking p2n[i]. This approach requires many unnecessary checks since all parts need to be considered even if most of them are not connected to the vertex. This problem is exacerbated when K is large, which is the case for many real-life applications. Overall, even when the cost of the intersection computation is O(1) per net, the algorithm takes O(K × |H|). This is not acceptable, since K can be in the order of thousands, and |H| can be huge for streaming, massive-scale hypergraphs. We observe that the performance problem in Algorithm 2 arises, because the connectivity among the nets and the parts is stored from the parts' perspective. However, if this information had been organized from the perspective of the nets it would be possible to identify the active parts that are connected to at least one net of the current vertex v at hand. Algorithm 3 describes the pseudocode of this approach. For an efficient computation, it uses four auxiliary arrays, save, mark, pids, and indx. Each of these arrays is of size K. However, they are only initialized once, and no expensive reset operation with complexity Θ(K) is performed after a vertex is processed. Thanks to these arrays, when a part is first identified to be connected to one of the nets in nets[v], it is marked to save a single net and placed into the active part array. Once it is placed, the next access to the same part (but for a different net) will only increase the savings of this part. Both these operations can be performed in constant time and no loop over all the parts is required.

Algorithm 2: MINMAX

Input: K, s Output: part[•] for i from 1 to K do p2n[i] ← ∅ for all (v, nets[v]) in streaming order do saved ← -1 for i from 1 to K do if |V i | -|V pmin | < s then if |p2n[i] ∩ nets[v]| > saved then saved ← |p2n[i] ∩ nets[v]| p ← i part[v] ← p p2n[p] ← p2n[p] ∪ nets[v] if p = p min then Update p min return part Algorithm 3: MINMAX-N2P Input: K, s Output: part[•] save[•] ← an array of size K mark[•] ← an array of size K with all -1s pids[•] ← an array of size K indx[•] ← an array of size K for all (v, nets[v]) in streaming order do active ← 0 for n ∈ nets[v] do if n appears for the first time then n2p[n] ← ∅ for i ∈ n2p[n] do if mark[i] = v then mark[i] ← v active ← active + 1 pids[active] ← i save[active] ← 1 indx[i] ← active else save[indx[i]] ← save[indx[i]] + 1 saved ← -1 for j from 1 to active do i ← pids[j] if |V i | -|V pmin | < s then if save[j] > saved then saved ← save[j] p ← i part[v] ← p for n ∈ nets[v] do n2p[n] ← n2p[n] ∪ {p} ; if p = p min
MINMAX as proposed in [START_REF] Alistarh | Streaming min-max hypergraph partitioning[END_REF], and as described in Algorithm 2, has been used in the literature to benchmark novel algorithms for scalable hypergraph partitioning [START_REF] Mayer | Hype: Massive hypergraph partitioning with neighborhood expansion[END_REF], [START_REF] Epple | Partitioning Billionscale Hypergraphs[END_REF]. For instance, it is reported that a hypergraph with 0.43M vertices and 180M pins is partitioned into K = 128 parts in around 1000 seconds [START_REF] Mayer | Hype: Massive hypergraph partitioning with neighborhood expansion[END_REF]. On the other hand, the variant in Algorithm 3 can partition a hypergraph with 1.1M vertices and 228M pins into K = 2048 parts in around 200 seconds.

Considering that K is at most in the order of tens of thousands, the extra memory due to the four additional arrays is not prohibitive. On the other hand, both MINMAX and MINMAX-N2P use approximately the same amount of memory to store the connectivity information. That is the total number of entries in n2p[•] and p2n[•] arrays are the same, and exactly n more than the (connectivity -1) metric given in Eq. 2. This being said, MINMAX-N2P uses slightly more memory since unlike K, |N | is not known beforehand in the streaming setting, and a dynamic data structure with more overhead is required to organize the connectivity as in n2p.

C. MINMAX Variants Using Less Memory

Depending on the hypergraph and the number of parts, the memory requirements of the two previous approaches can be too large. For streaming data, there is no upper limit, which is obviously a problem for SBCs. In this subsection, we propose variants of MINMAX that can use a fixed amount of memory.

1) MINMAX-L : As explained above, the total memory spent for n2p grows as long as the data is streaming. To avoid this, while working similar to Algorithm 3, MINMAX-L restricts the maximum length of each n2p

[•] to . When a new part p is being added to a n2p[n] for a net n ∈ N , if |n2p[n]| = , a random part id from n2p[n]
is chosen and replaced with p. Although the connectivity information is only kept in a lossy fashion, we expect it to guide the partitioning decisions for sufficiently large values.

2) MINMAX-BF: Bloom filters (BF) are memory-efficient and probabilistic data structures used to answer set membership queries [START_REF] Bloom | Space/time trade-offs in hash coding with allowable errors[END_REF]. Compared to the traditional data structures

Algorithm 4: MINMAX-L Input: K, s, Output: part[•] • • • same as Algorithm 3 for all (v, nets[v]) in streaming order do • • • same as Algorithm 3 for n ∈ nets[v] do if |n2p[n]| < then n2p[n] ← n2p[n] ∪ {p} else if p / ∈ n2p then idx ← rand(1, ) n2p[idx] ← p if p = p min then
Update p min return part such as arrays, sets, or hash tables used for the same task, a BF occupies much less space while allowing false positives.

If the item is a member of the set, BF always returns true. If the item is not in the set, BF most likely answers correctly but can also return a false positive.

A Bloom Filter, which employs an m-bit sequence, uses k hash functions to find the indices of bits and sets them to 1 to mark the existence of a new. To answer a query for an item x, it simply checks the corresponding k hash functions on x, and answers positively if each of the corresponding bits is 1. An important parameter for a BF is its false positive probability. Assuming the hash functions are independent of each other, when n items are inserted into a BF, kn bits are altered. Hence, the probability of a bit staying zero is

(1 -1/m)
kn ≈ e -kn/m , and the false-positive probability can be computed as 1 -e -kn/m k .

The BF variant of MINMAX is similar to Algorithm 2. However, instead of p2n, it leverages a Bloom Filter BF to store connectivity tuples (n, p) which means that the net n has a pin at part p. For a given (v, nets[v]), it goes over all the parts, and for each net in nets[v], it queries the corresponding tuple within the BF. Then it chooses the part with the most number of positive answers. The pseudocode of this approach is given in Algorithm 5. As in MINMAX-L , using a BF limits and fixes the amount of memory that will be used to store the connectivity among the nets and the parts.

3) MINMAX-MH: An approach fundamentally different from MINMAX-L and MINMAX-BF is to completely discard the connectivity information and try to cluster similar vertices with similar nets[•] sets into the same parts. A natural tool for this task is hashing. We use a MinHash-based approach [START_REF] Broder | On the resemblance and containment of documents[END_REF], MINMAX-MH, to find the part id for a given vertex. For implementation, we use k hash functions in the form of h i (x) = a i x + b i mod q where 1 ≤ i ≤ k, q is a prime number, a i and b i are random integers chosen from [0, q) per hash function. Given (v, nets[v]), this approach first computes (α 1 , α 2 , . . . , α k ) where

α i = min n∈nets[v] {h i (n)}. Then the Algorithm 5: MINMAX-BF Input: K, s, m Output: part[•] BF ← ∅ (creates an m-bit, all zero sequence) for all (v, nets[v]) in streaming order do saved ← -1 for i from 1 to K do if |V i | -|V pmin | < s then saved i ← 0 for n ∈ nets[v] do if BF.query((n, i)) then saved i ← saved i + 1 if saved v > saved then saved ← saved i p ← i part[v] ← p for n ∈ nets[v] do BF.insert((n, p)) if p = p min then Update p min return part part id for v is computed as p = k i=1 α i mod K.
If v cannot be assigned to p due to the balance constraint, the approach tries the next part in the natural order, (p+1 mod K) until a suitable part is found.

It is intuitive to think that vertices with similar net sets will end up with closer hash values which will result in assigning them in the same part. The fact that there is no need for additional memory to keep net-part connectivity unlike the previous algorithms makes this approach suitable for low memory environments.

D. Buffering and Refining

Revoking the partitioning decisions is impossible for the strict streaming setting. This is so as the net sets are forgotten after the corresponding vertex is put to a part. We can overcome this by using an additional buffer to keep the nets[•] sets. With such a buffer, one can revisit the buffered vertices and reassign them to a different part if such a reassignment improves the cutsize. A high-level description of this approach is given in Algorithm 6.

Algorithm 6 works along the same lines of MINMAX-N2P. In addition to finding the part id for a vertex v, after v is processed it can be chosen to be inserted to the buffer BUF. Once the buffer is full, the algorithm goes over all the vertices in the buffer passes times. For each vertex u, first the leaveGain of u is computed which is the change in the (connectivity -1) metric when u is removed from part[u]. Then if u is decided to be movable, it is removed from part[u]. A new part p is then found and u is moved to p. We desgined three strategies, To compute leaveGain, one also needs to keep track of the number of pins of each net residing at each part. This almost doubles the memory requirement of the refinement heuristics compared to MINMAX, since for every connectivity entry stored in n2p , an additional positive integer is required. That is the original entry shows that "net n is connected to part p", and the additional entry required for refinement shows that "net n has k pins in part p". The refinement-based algorithms require this information since when k = 1, they can detect a gain on the connectivity.

IV. RELATED WORK

There exist excellent offline hypergraph partitioners in the literature such as PaToH [START_REF] Aykanat | PaToH (Partitioning Tool for Hypergraphs)[END_REF] and HMetis [START_REF] Karypis | Multilevel hypergraph partitioning: applications in vlsi domain[END_REF]. Recently, more tools are developed focusing on different aspects and using different approaches. Deveci et al. [START_REF] Deveci | Hypergraph partitioning for multiple communication cost metrics: Model and methods[END_REF] focus on handling multiple communication metrics at once, Mayer et al. [START_REF] Mayer | Hype: Massive hypergraph partitioning with neighborhood expansion[END_REF] focus on the speed, and Schlag et al. [START_REF] Schlag | k-way hypergraph partitioning via n-level recursive bisection[END_REF] focus more on the quality by using a more advanced refinement mechanism. Faraj et al. [START_REF] Faraj | Buffered streaming graph partitioning[END_REF] have recently proposed a streaming graph partitioning algorithm which yields high quality partitioning on streaming graphs utilizing buffering model. Jafari et al. [START_REF] Jafari | Fast shared-memory streaming multilevel graph partitioning[END_REF] have proposed a fast parallel algorithm which processes the given graph in batches.

The streaming setting has not been analyzed thoroughly on hypergraphs except the work by Alisarth et al. [START_REF] Alistarh | Streaming min-max hypergraph partitioning[END_REF] which proposes the original MINMAX algorithm at once. We improved this algorithm significantly by altering its inner data structure used to store the part-to-net connectivity. Furthermore, we proposed techniques to refine the existing partitioning at hand with the help of some extra memory.

V. EXPERIMENTAL RESULTS

We tested the proposed algorithms on hypergraphs created from six matrices downloaded from the SuiteSparse Matrix Collection (https://sparse.tamu.edu/). For each n × n matrix, we create a column-net hypergraph H = (V, N ) where the vertices (nets) in V (in N ) correspond to the rows (columns) of the matrix. The properties of the resulting hypergraphs are given in Table I. Moreover, we tested the algorithms on a cutting-edge server and multiple SBCs. Since the variants studied in this work have different time, memory, and quality tradeoff characteristics, we used a set of single-board computers to observe their performance. These SBCs have different number of cores and different amounts of memory. The specifications of these architectures are as follows:

• Server: Intel Xeon Gold 6140, 2.30GHz, 256GB RAM, gcc 5.4.0, Ubuntu 16.04.6.

• LattePanda: Intel Atom x5-Z8350, 1.44GHz, 4096MB
DDR3L @ 1066 MHz RAM, gcc 5.4.0, Ubuntu 16.04.2 • Pi: Broadcom BCM2837, 1.2GHz, 1024MB LPDDR2 @ 400 MHz RAM, gcc 6.3.0, Raspbian 9 • Odroid: ARM Cortex-A15, 2GHz and Cortex-A7 @1.3GHz, 2048MB LPDDR3 @ 933 MHz, 2048MB LPDDR3 @ 933 MHz, gcc 7.5.0, Ubuntu 18.04.1 We implemented algorithms in C++ and compiled on each device separately with the gcc version available in each SBC. For each hypergraph, we created three random streams with different vertex orderings. Although it is an offline partitioner and is not suitable for the streaming setting, we used PaToH v3.3 [START_REF] Aykanat | PaToH (Partitioning Tool for Hypergraphs)[END_REF] to evaluate the quality and performance of the streaming algorithms. As balance constraint, we used a dynamic value s in Eq. 1 computed by using a constant allowed imbalance ratio β = 0.1, for which values in the range 5%-10% are common in the partitioning literature. That is, while processing the rth stream element, the allowed slack is set to s = max(1, β × (r/K)). This translates to obtaining balance by settign the slack s to β times the average part weight at any time during partitioning. The first two rows show the effectiveness of the modification on MINMAX. The modified version MINMAX-N2P is ≈ 1700x faster on this hypergraph while using 16-18MB more memory. This is due to the reduced number of unnecessary operations on unrelated parts while computing the part savings. The MINMAX-L3 and MINMAX-L5 variants obtain partitions comparable to MINMAX in terms of quality (i.e., with respect to connectivity-1) while using less memory.

For MINMAX-BF and MINMAX-MH, we use 4 and 16 hash functions, and for the former, we use 20M bits. Although being fast, the partitioning quality of the MinHash variant is half of the RANDOM. On the contrary, the Bloom Filter variant seems to work well with a comparable partitioning quality. However, it is slow since when BF is used instead of n2P, one needs to go over all the parts to compute their savings. Still, it is very promising in terms of memory/quality trade-off and enables a scenario in which a device with a small memory is used in partitioning on a network edge.

Figure 1 compares MINMAX-N2P's performance to those of a hypothetical streaming tool based on PaToH. That is the hypergraphs are partitioned by PaToH (with SPEED and DEFAULT configurations) and the run times, cuts (connectivities), and memory requirements are normalized with respect to those of MINMAX-N2P for hypergraphs coPapersDBLP, holywood2009, socLiveJournal1, com-Orkut, and uk2005. The experiments show that on these hypergraphs, the DEFAULT configuration can be 10-30% better in terms of partitioning quality. However, it can also be more than 100× slower (see the run time ratio bar for SocLiveJournal1). Furthermore, PaToH uses 12-18× more memory compared to MINMAX-N2P. Note that PaToH, or any other offline partitioner, is not suitable for the streaming setting. The results are only given to demonstrate that there is room for improvement on MinMax-n2p especially in terms of partitioning quality. This in turn justifies the attempts for refining the partitioning throughout streaming. III presents the results of these experiments with a buffer capacity B = 0.15× |H| for each hypergraph |H|. The results show that although refinement can be useful for reducing the connectivity, its overhead is significant. Furthermore, after 2 passes, there is only a minor improvement on the partitioning quality. Using passes = 4 reduces the cut size for a large number of experiments, and therefore we use this setting in the rest of the experiments.

Figure 2 shows the run times (in seconds), connectivity scores (normalized with respect to MINMAX-N2P), and memory usages (in MBs) of the refinement heuristics, MINMAX-N2P, and MINMAX-L5 on all hypergraphs in Table I. The experiments are executed on the Server. An imbalance ratio of 0.1 and passes = 4 are used for the experiments. As the results show, the refinement based heuristics improve the partitioning quality in between 5-20% depending on the hypergraph. Furthermore, when the buffer size is increased, these heuristics tend to improve the quality better. Besides, for the two largest hypergraphs in our experiments, REF RLX SV is much faster than the other two refinement heuristics REF and REF RLX with a similar improvement over the partitioning quality and the same memory requirement. Hence, it can be a good replacement to the original MINMAX-N2P heuristic with no refinement, if the partitioning quality has the utmost importance.

Table IV shows the run time performance of the proposed algorithms on different architectures and for K = 32, 256, 2048 and 16384. The hypergraph socLiveJournal1 is used for these experiments. The results are similar to the ones in the Server. Yet additionally, the slow-down values in the last three columns show that using much less memory, MINMAX-L5 stays more scalable compared to other algorithms when K is increased. Furthermore, the overhead of the refinement heuristics degrades the scaling behavior when they are added on top of the MINMAX-N2P. However, their negative impact tends to decrease when an SBC with more memory is used. This also shows the importance of streaming hypergraph algorithms with low-memory requirements in practice.

VI. CONCLUSION AND FUTURE WORK

We focused on the streaming hypergraph partitioning problem. The problem has unique challenges compared to similar problem of streaming graph partitioning. We significantly improved the run time performance of a well-known streaming algorithm and proposed several variants on top of it to reduce the memory footprint and improve the partitioning quality. The experiments show that there is still room for improvement for these algorithms. As future work, we plan to devise more advanced techniques that can overcome the trade-off among the run time, memory requirements, and partitioning quality.

  then Update p min ; return part B. Using Net-to-Part Information and Finding Active Parts

Fig. 1 :Fig. 2 :

 12 Fig. 1: Run times, cuts and memory usages of PaToH normalized with respect to those of MINMAX-N2P. The experiments are executed on the Server.

TABLE I :

 I Hypergraphs used for the experiments.

	Matrix	Size	Pins	Max. Deg Avg. Deg	Deg. Var.
	coPapersDBLP	0.5M	30.5M	3.2K	56.41	66.23
	hollywood2009	1.1M 227.8M	11.4K	99.91	271.69
	soc-LiveJournal1	4.8M	68.9M	13.9K	14.23	42.30
	com-Orkut	3.0M 234.3M	33.3K	76.28	153.92
	uk-2005	39.4M 936.3M	1.7M	23.72	1654.56
	webbase2001	118.1M	1.0B	816.1K	8.63	141.79

TABLE II :

 II Comparison of the proposed streaming hypergraph partitioning algorithms and the state-of-the-art on the Server and coPapersDBLP.

	Parts
	256	Algorithm	Run Time(sec)	Cut(×10 6 ) Memory(MB)
		MINMAX	13952.12	1.731	18.60
		MINMAX-N2P	7.40	1.731	36.78
		MINMAX-L3	2.20	2.543	15.00
		MINMAX-L5	2.60	1.913	25.41
		MINMAX-BF(4)	496.90	3.280	2.64
		MINMAX-BF(16)	1479.19	3.273	2.54
		MINMAX-MH(4)	1.49	11.214	0.18
		MINMAX-MH(16)	5.37	13.464	0.19
		RANDOM	0.04	23.817	0.00
	2048 Algorithm	Run Time(sec)	Cut(×10 6 ) Memory(MB)
		MINMAX	17823.45	3.096	28.65
		MINMAX-N2P	10.64	3.096	44.10
		MINMAX-L3	2.61	5.217	15.00
		MINMAX-L5	3.25	4.168	25.10
		MINMAX-BF(4)	3664.48	4.508	2.59
		MINMAX-BF(16)	11882.30	4.655	2.57
		MINMAX-MH(4)	1.51	15.134	0.12
		MINMAX-MH(16)	5.38	18.101	0.15
		RANDOM	0.07	28.992	0.00
	Table II reports the run times, cut, and memory us-
	age of the proposed MINMAX variants on the hypergraph
	CoPapersDBLP.

TABLE III :

 III Effect of the number of passes on refinement algorithms; the results are averaged for K = 256 and K = 2048. The experiments are executed on the Server. All of the values are normalized with respect to those of MINMAX on the same experiment. To find a good value for passes in the refinement-based algorithms, we performed 2, 4, and 8 passes over the buffered vertices and measured the run time and partitioning quality of REF, REF RLX, and REF RLX SV. Table

					Run Time			Cut
	Matrix	Buf.	Passes	R	RR	RRS	R	RR	RRS
			2	2.8	4.5	2.1	0.84	0.82	0.87
	coPapers	0.15	4	3.8	6.6	2.8	0.83	0.81	0.87
			8	6.0	10.8	4.2	0.83	0.79	0.86
			2	4.0	4.4	3.8	0.97	0.96	0.97
	hollywood	0.15	4	5.8	6.4	5.6	0.97	0.95	0.96
			8	9.6	10.2	8.8	0.96	0.95	0.97
			2	3.8	3.8	3.4	0.96	0.94	0.95
	soc-Live	0.15	4	5.5	5.7	4.9	0.96	0.94	0.94
			8	9.0	9.3	7.8	0.96	0.94	0.94

TABLE IV :

 IV Run times of the proposed algorithms on single board devices and socLiveJournal1 for different K values. The allowed imbalance is β = 0.1, and the buffer capacity B = 0.15×|H| is used. The results on the K = 32 column are given in seconds, and for K = 256, 2048 and 16384, the results are normalized with respect to K = 32. That is the results in the last three columns show the decrease in the run time when K is increased from 32 to the corresponding value.

				Parts
	Device	Algorithm	32	256	2048	16384
		MINMAX-N2P	142.6	1.43×	2.16× 3.73×
	Pi-1GB	REF REF RLX	593.1 1.64× 650.0 1.50×	2.93× 4.38× 2.89× 4.16×
		REF RLX SV	532.3 1.53×	2.89× 4.73×
		MINMAX-L5	55.6 1.17×	1.30× 1.74×
		MINMAX-N2P	72.4	1.32×	1.81× 2.80×
	Odroid-2GB	REF REF RLX	335.2 1.41× 362.8 1.33×	2.32× 3.62× 2.27× 3.60×
		REF RLX SV	287.5 1.47×	2.37× 3.92×
		MINMAX-L5	36.1 1.11×	1.19× 1.40×
		MINMAX-N2P	71.5	1.34×	1.89× 2.85×
	LattePanda-4GB	REF REF RLX	308.3 1.45× 339.8 1.35×	2.46× 3.73× 2.43× 3.62×
		REF RLX SV	264.8 1.52×	2.57× 4.05×
		MINMAX-L5	30.4 1.11×	1.18× 1.33×

ACKNOWLEDGEMENTS

This work is funded by The Scientific and Technological Research Council of Turkey (T ÜBITAK) under the grant number 117E249.


R Priam
email: rpriam@gmail.com

Latent block principal component analysis for binary tables

Keywords: Bernoulli latent block mixture model, co-clustering, pca, binary matrix, BEM algorithm

published or not. The documents may come

Latent block principal component analysis for binary tables

Introduction

For exploratory data analysis of binary tables, several methods have been proposed in the literature from diverse points of view. We are mainly interested in matricial approaches related to the principal component analysis (pca) [START_REF] Lebart | Multivariate Descriptive Statistical Analysis[END_REF] for binary tables for reducing the two dimensions of these tables. When the data matrix is defined from the set I of objects (rows, observations) and the set J of variables (columns, attributes), the dimensions are these two sets or their sizes. Reducing the dimensionality is able to show the rows or the columns in a subspace dramatically smaller than the former space, a visualization becomes possible but also a post analysis such as a discriminant one. This is how pca is able to visualize the rows and the columns on the plane thanks to this property. Several linear methods exist in the literature for binary data, such as correspondence analysis [START_REF] Lebart | Multivariate Descriptive Statistical Analysis[END_REF][START_REF] Benzecri | Correspondence Analysis Handbook[END_REF] or sparse logistic pca [START_REF] Lee | Sparse logistic principal components analysis for binary data[END_REF][START_REF] Lee | A coordinate descent MM algorithm for fast computation of sparse logistic pca[END_REF] and also [START_REF] Tipping | Probabilistic visualisation of high-dimensionnal binary data[END_REF][START_REF] Schein | A Generalized Linear Model for Principal Component Analysis of Binary Data[END_REF][START_REF] De Leeuw | Principal component analysis of binary data by iterated singular value decomposition[END_REF]. Here, a nonlinear version is of main interest because the linearity alone is not able to find the true relations betweens the objects from I or from J in a synthetic way. For improving the limits due to the linear setting, clustering is a powerful choice able to better discover the underlying shape of the data cloud, such was proposed in [START_REF] Tipping | Mixtures of probabilistic principal component analyzers[END_REF] with a mixture of probabilistic pca. A family of methods has been defined for visualization, it is based on a clustering procedure with vicinity constraints for the clusters. These methods are seen in the literature as a nonlinear pca with a projection on a surface instead of the linear projection on the plane in pca. This family is named self-organizing map, it generalizes a former algorithm called the Kohonen's map [START_REF] Kohonen | Self-organizing maps[END_REF] which is a k-means method which integrates the constraints at the level of the mean centers. For a sparse binary table, this idea was introduced with an alternative clustering model for co-clustering, the latent block model [START_REF] Govaert | Clustering with block mixture models[END_REF], by a method called herein blockgtm [START_REF] Priam | Topographic bernoulli block mixture mapping for binary tables[END_REF] for visualization of binary tables. Generally a coclustering is dramatically more efficient than a clustering of the two sets separately because the number of parameters is parsimonious and the clustering of one dimension is involved [START_REF] Bock | Simultaneous clustering of objects and variables[END_REF][START_REF] Govaert | Classification croisée[END_REF][START_REF] Govaert | Simultaneous clustering of rows and columns[END_REF] in the clustering of the other dimension. In order to get a reduction method instead of a visualization method, we propose to change the model by adding another set of latent vectors to be estimated. By this way the symmetry of the learning is retrieved as explained next after. In this paper, we present in section 2 the proposed model for binary data with its clustering foundation followed by the parameterization. In Section 3, an algorithm is given for the inference of the parameters with the new constraints. The last section is the concluding part with perspectives. such reducing purpose, we propose the new method named blockpca or latent block principal component analysis by rewriting the parameters from the Bernoulli block mixture model with these new latent vectors and adding a penalization for avoiding the empty clusters. The inference of the parameters is based on a generalized expectation maximization (EM) algorithm such as blockgtm. In this section, the model is presented before the inference of the parameters.

Latent block mixture model

In the following, the n × d data table is defined by x = {(x i j ); i ∈ I and j ∈ J} where x i j ∈ {0, 1}. The purpose of a method of block clustering is to partition the table by blocks with a distribution of the cells for each block. In [START_REF] Govaert | Clustering with block mixture models[END_REF][START_REF] Govaert | An EM algorithm for the block mixture model[END_REF], the methods from Govaert [START_REF] Govaert | Classification croisée[END_REF][START_REF] Govaert | Simultaneous clustering of rows and columns[END_REF] are embedded in a probabilistic mixture [START_REF] Symons | Clustering criteria and multivariate normal mixture[END_REF]. The resulting model is called the latent block mixture model (LBM) [START_REF] Govaert | Clustering with block mixture models[END_REF][START_REF] Govaert | An EM algorithm for the block mixture model[END_REF] which is a mixture model which sums over the crossproduct of two sets of assignments Z × W . Here Z and W denote the sets of all possible assignments z of I and w of J. The observed x i j cells are generated by n × d random variables which are assumed to be independent if z and w are known. Let also define p = (p 1 , . . . , p g ) (resp. q = (q 1 , . . . , q m )) for the vectors of probabilities p k (resp. q ) that a row (resp. a column) belongs to the k th component (resp. th component). Hence, this leads to the pdf:

f (x; θ) = ∑ (z,w) ∏ i,k p z ik k ∏ j, p z j ∏ i, j,k, ϕ(x i j ; α k ) z ik w j .
Here, ϕ(.; α k ) is the pdf in a block defined for x i j ∈ {0, 1} and the unknown parameter α k . They are aggregated into the vectorial parameter θ which is written from α = (α 11 , . . . , α gm ) with a vectorial form, p and q. For binary data the pdf is a Bernoulli law which leads to the Bernoulli latent block model [START_REF] Govaert | Clustering with block mixture models[END_REF]. In a block, the pdf of a cell is defined as follows,

ϕ(x i j ; α k ) = (α k ) x i j (1 -α k ) 1-x i j .
Next, the parameters α k are written with a function of latent vectors for reduction purpose, and a learning algorithm is given for this new model.

Parameterization with latent spaces

Herein, the parameters α k of the block mixture model are modeled with two sets of vectors. To keep the dependence on and k of α k , when h ∈ N + is a positive integer, it is introduced inner products defined as: v T u k from the latent variables {v } where v ∈ R h and {u k } where u k ∈ R h . Each product is mapped into a probability with a sigmoid function σ(.):

α k = exp(v T u k ) 1 + exp(v T u k )
.

This leads to suitable parameters for the block model. With this formulation, the latent parameter α, in a matricial form of size g × m, is replaced by the two latent matrices U and V ,

V = [v 1 |v 2 | • • • |v m ] T U = [u 1 |u 2 | • • • |u g ] T .
This model remains parsimonious because h is small in practice.

Algorithm BEM

In this section, the estimation of the parameter θ is studied via the likelihood maximization of the latent block mixture model. For this purpose an algorithm called BEM [START_REF] Govaert | Clustering with block mixture models[END_REF][START_REF] Govaert | An EM algorithm for the block mixture model[END_REF] generalizes the usual EM algorithm from mixture model for one partition to two partitions (or more) instead of just one. The new completed data is the vector (x, z, w) where the unobservable vectors z and w are the row and column labels. It is denoted the posterior probabilities c ik that a data row i is attributed to a row cluster k and d j that a data column j is attributed to a column cluster . The algorithm has two steps, like EM, which are iterated until convergence.

-E-step: The posterior probabilities are found by solving the same problem than in the unconstrained case but with the new parameterization, such that:

c (t) ik ∝ ∏ (σ(v T u k )) u i (1 -σ(v T u k )) d -u i , (1)
d (t) j ∝ ∏ k (σ(v T u k )) v jk (1 -σ(v T u k )) c k -v jk . (2)
-M-step: Here the proportions are supposed not equal. This results into the two following criteria. For the matrix U, let rewrite the objective function as the sum:

Q(θ|θ (t) ) = ∑ k Qk (u k |θ (t) ) Qk (u k |θ (t) ) = ∑ y (t) k v T u k -c (t) k d (t) log(1 + e v T u k ) . (3)
For the matrix V , let rewrite the objective function as the sum:

Q(θ|θ (t) ) = ∑ Q (v |θ (t) ) Q (v |θ (t) ) = ∑ k y (t) k v T u k -c (t) k d (t) log(1 + e v T u k ) . (4)
The optimization problem for maximizing the criteria (3) or (4) depending on U and V can be performed also after their corresponding a posteriori probabilities.

Here, it has been denoted the quantities y

(t) k = ∑ i, j c (t) ik d (t) j x i j , u (t) i = ∑ j d (t) j x i j , d (t) = ∑ j d (t) j and v (t) jk = ∑ i c (t) ik x i j , c (t)
k = ∑ i c (t)
ik . A limit with the objective function is the problem of empty clusters which comes with a larger number of blocks than the natural one, which is not suitable with our model. Hence, it is added a penalization with η a constant value in order to insure that the clusters have all different positions in the space,

Qk (u k |θ (t) ) = Qk (u k |θ (t) ) + η 2 ∑ k ||u k -u k || 2 Q (v |θ (t) ) = Q (v |θ (t) ) + η 2 ∑ ||v -v || 2 .
(

) 5
Note that cancelling η allows to retrieve the former form for the function to optimize, in case of the penalization is not required. Alternative penalizations may be possible but not tested herein, without loss of generality. For the estimation of the parameters, the proposed algorithm is further presented with the analytical expressions of the updates at the M-step next section.

M-step and detailed algorithm

An algorithm for the inference of the parameters is presented in this section with first the approxima-tions in stake for a matricial expression of the updates at each iteration, followed by all the steps of the algorithm.

Two quadratic approximations

At the M-step a nonlinear function is optimized because of the sigmoid functions. The problem to be solved is the following maximization:

θ (t+1) = argmax θ Q(θ|θ (t) ) .
For a solution to this problem, a second order approximation of Q is presented. The symmetry of the two original mirrored formulas in BEM for each side of the table is not anymore lost on the contrary to blockgtm where only rows were mapped, hence only one side is explained here for V , while the other side for U is found by transposition. Two approaches are considered next, a Newton-Raphson approximation (NR) and a variational approximation.

-M-step with a 2 nd order approximation: A NR algorithm is considered by using the first order and second order derivatives of the criterion. This approach is generally able to increase the log-likelihood at each step of EM. There is an interaction between the vectors v because of the penalization, hence the problem cannot be solved for each separately. Without the penalization, the Hessian matrix for the whole set of parameters is block-diagonal H = diag (H ) with a block H for each . With the penalization, the gradient vectors q v are not anymore defined for each v independently. It is denoted

C = (c (t) ik ) a g × n matrix of posterior proba- bilities, V = (v (t) jk ) a g × d matrix of sufficient statistics, G = (c (t) k ) and F = (α (t) k ) are g × g diagonal matrices, α = (α (t) k ) a g × 1 vector, u = (u (t) i ) a n × 1 vector, d = (d (t) j ) a d × 1 vector, and d ( ) = ∑ j d (t) j is a scalar. Let's b (t)
and H (t) stand respectively for the first and second order derivatives of Q at time t. A local quadratic approximation of the objective function in a vicinity of the current solution denoted v (t) is as follows when keeping the penalization

      diag                  H (t) 1 -2ηI h H (t) 2 -2ηI h . . . H (t) m -2ηI h                  + η       I h I h • • • I h I h I h . . . . . . . . . . . . . . . I h I h . . . I h I h                  v 1 v 2 . . . v m      =       b (t) 1 -H (t) 1 v(t) 1 b (t) 2 -H (t) 2 v(t) 2 . . . b (t) m -H (t) m v(t) m       . (6)
Figure 1: Linear regression for the M-step.

which is already quadratic:

Q (v |θ (t) ) ≈ Q (v (t) |θ (t) ) + (v - v(t) ) T b (t) + 1 2 (v - v(t) ) T H (t) (v - v(t) ) + η 2 ∑ ||v -v || 2 . (7)
The NR algorithm performs the maximization of a local quadratic approximation instead of the former objective function. The matrix and the vectors are as follows:

H (t) = -d ( ) U T GF (I h -F )U b (t) = U T y -d ( ) F g k v(t) = v (t) . (8)
Note that the dimensions of the matrices H remain small even when the data matrix is large because only the size h is involved. Another way to perform the optimization is to approximate the sigmoid function with a quadratic function as explained next.

-M-step with a variational bound: For a variational approach (VR), the function ϕ may be approximated for instance. In this case, a bound on the sigmoid function [START_REF] Saul | Mean field theory for sigmoid belief networks[END_REF] might be relevant. By convexity, this is written as a upper bound of σ(a) as an exponential of a quadratic function of a and ε multiplied with σ(ε) where ε is an unknown parameter to be computed. Here a is a scalar in R, while ε is the variational parameter, and λ(ε) = 1 4ε tanh ε 2 . As previously the objective function cannot break into m different criteria which can be maximized independently. Let's denote ε = (ε 1 , ε 2 , • • • , ε g ) T where at the iteration time each component is computed

ε (t) k = |(u (t) k ) T v (t) |.
Then, by using the bound on the sigmoid, the criterion to optimize may be written as follows:

Q (v |θ (t) ) ≥ +(v - v(t) ) T b(t) -1 2 (v - v(t) ) T H(t) (v - v(t) ) +d (t) ∑ k c (t) k log σ(ε k ) -d (t) ( ) ∑ k c (t) k λ(ε k )(-ε 2 k ) + 1 2 ε k + η 2 ∑ ||v -v || 2 = Q (v , ε |v (t) , ε (t) ) . (9)
Here, g k denotes the column vector compound of the elements ∑ i c

(t)
ik and y is the column vector compound of the elements y (t) k . Let's also have Λ for the diagonal matrix with diagonal elements equal to λ k = λ(ε k ). The expression is as for the NR procedure with new corresponding quantities:

b(t) = U T y - d ( ) 2 g k H(t) = 2 d ( ) U T GΛ U v(t) = 0 h . (10)
As for NR it is obtained a new maximization problem without a multidimensional nonlinear function involved but instead several quadratic problems which can be performed analytically. The dimensionality of the problem is only of order h instead of the size of the table.

Algorithm for learning the parameters

The two algorithms lead to the same form of quadratic approximation, such as finally, the gradient w.r.

t. v = vect(V ) = (v T 1 , • • • , v T m )
T is written by the derivative of the criterion. The two version of the M-step (NR and VR) have been unified under a similar notation and a similar regularization.

-

Initialization: Initialize {c (0) ik }, {d (0) j }, U (0) , V (0) . -E-Step: Update {c (t)
ik } by (1), and {d (t) j } by (2), alternatively.

-M-Step: Update {ε k },{ε } if needed and U,V from (6), Here, equating the gradient with zero leads to m equations. Matricially, this leads to a linear regression of the form Ω v v = b v with the expression given at figure 1. The same expression is available for the variational approximation by replacing b (t) by b(t) and H (t) by H(t) while v(t) depends on the algorithm. As in blockgtm, the linear regression may induce an eventual overfitting when the dimensions h, g or m increases. Recent developments in the literature have shown the good behavior of a parsimonious solution via a penalization of the criterion: the LASSO regression [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] adds a penalizing term -η v |v| to the sum of squares in order to induce automatically the sparsity of the regression coefficients. At the convergence, all the quantities computed at the final step are denoted with a hat. The algorithm is summarized at the figure 2, by giving all the steps of the algorithm, from the initialization to the stopping function at the convergence with ε a small constant.

V ← argmin V ||Ω v v -b v || 2 + η v |v| U ← argmin U ||Ω u u -b u || 2 + η u |u| . -End: If |v (t+1) -v (t) | + |u (t+1) -u (t) | < ε then stop. Else return to E-Step.

Conclusion

Herein it is explained how to model a pca method with a co-clustering for large sparse binary tables. The resulting reductions for the two dimensions of the table are written with the latent matrices Û, V and the posterior probabilities as follows: xi = ∑ k ĉik ûk x j = ∑ d j v , where xi is the reduced i th row and x j is the reduced j th column, as h-dimensional vectors where h n and h d. A perspective is other types of variables [START_REF] Priam | Generalized topographic block model[END_REF] such as in contingency tables for textual data [START_REF] Priam | Negative binomial latent block model with generalized constraints[END_REF], more constraints such as orthogonality or alternative reducing parameterization.

Figure 2 :

 2 Figure 2: Algorithm for learning the parameters in the binary blockpca with a Newton-Raphson or variational algorithm for the M-step.

Proposed model and criteriaFor principal component analysis, the clusters depend on latent vectors, from the two dimensions (rows and columns) in order to reduce the data table from the two sides. Thus in order to share


A Gloystein

N Nilius

J Goniakowski
email: jacek.goniakowski@insp.jussieu.fr

C Noguera

The Nanopyramidal Reconstruction of Cu 2 O(111): a Long-standing Surface Puzzle Solved by STM and DFT

Despite the importance of cuprous oxide for industrial applications, the precise structure of its most stable reconstructed surface remains unknown. Here, by a combination of highresolution STM measurements and simulations based on a DFT-HSE approach, we challenge the previous interpretations of the Cu 2 O(111) ( √ 3× √ 3)R30 • surface reconstruction, based on the removal of 1/3 of the under-coordinated surface oxygen atoms. We show that these models do not satisfy the shamrock-type features found by STM, i.e. their position in the Cu-O ring centers, the orientation towards the surface O cus atoms and the topographic height in the STM images. Moreover, they are thermodynamically unstable with respect to moving the vacancy subsurface, a property not recognized before which seems specific to the Cu 2 O (111) surface. We propose a nanopyramidal model of the ( √ 3× √ 3)R30 • surface reconstruction which is free from all these shortcomings. Here, the shamrock protrusions seen in STM are formed by three copper adatoms located at the center of the Cu-O rings and capped by an oxygen atom. This structure profoundly differs from existing models of the Cu 2 O(111) surface, and would change common perceptions on its reactivity.

Introduction

Cuprous oxide Cu 2 O is a widespread oxide formed in the first stages of corrosion of copper. It is a p-type semiconductor, with applications in photovoltaics and photocatalysis. As a catalyst it is widely used in water-gas shift reactions, methanol synthesis and partial oxidation. Precise knowledge of its surfaces properties is thus of paramount importance to foster its applications. This is why Cu 2 O surfaces have been extensively studied. [START_REF] Schulz | Photoemission and low-energy-electron-diffraction study of clean and oxygen-dosed Cu 2 O (111) and (100) surfaces[END_REF][START_REF] Nygren | Theoretical models of the polar Cu 2 O (100) Cu + -terminated surface[END_REF][START_REF] Mcclenaghan | A density functional theory study of the surface relaxation and reactivity of Cu 2 O (100)[END_REF][START_REF] Le | Reactivity of the Cu 2 O (100) surface: Insights from first principles calculations[END_REF][START_REF] Soldemo | The surface structure of Cu 2 O (100)[END_REF][START_REF] Zhang | Atomistic determination of the surface structure of Cu 2 O (111): experiment and theory[END_REF][START_REF] Jensen | Two new oxygen induced reconstructions on Cu (111)[END_REF][START_REF] Sträter | Growth and surface properties of cuprous oxide films on Au (111)[END_REF][START_REF] Ly | Growing ultrathin Cu 2 O films on highly crystalline Cu (111): a closer inspection from microscopy and theory[END_REF][START_REF] Soon | Surface oxides of the oxygencopper system: Precursors to the bulk oxide phase?[END_REF][START_REF] Islam | First principles investigation on the stabilization mechanisms of the polar copper terminated Cu 2 O (111) surface[END_REF][START_REF] Li | Stability and electronic properties of the O-terminated Cu 2 O (111) surfaces: Firstprinciples investigation[END_REF][START_REF] Nilius | Incorrect DFT-GGA predictions of the stability of non-stoichiometric/polar dielectric surfaces: The case of Cu 2 O (111)[END_REF][START_REF] Önsten | Atomic structure of Cu 2 O (111)[END_REF][START_REF] Bendavid | Firstprinciples predictions of the structure, stability, and photocatalytic potential of Cu 2 O surfaces[END_REF] Cu 2 O surface structures, mainly along the (100) and (111) orientations have been studied by a variety of surface science methods. [START_REF] Schulz | Photoemission and low-energy-electron-diffraction study of clean and oxygen-dosed Cu 2 O (111) and (100) surfaces[END_REF] Several LEED patterns have been recorded at the (100) surface, featuring (1 × 1), c(2 × 2) and (3,0;1,1) periodicities, and some model structures based on STM measurements and DFT simulations have been proposed. [START_REF] Nygren | Theoretical models of the polar Cu 2 O (100) Cu + -terminated surface[END_REF][START_REF] Mcclenaghan | A density functional theory study of the surface relaxation and reactivity of Cu 2 O (100)[END_REF][START_REF] Le | Reactivity of the Cu 2 O (100) surface: Insights from first principles calculations[END_REF][START_REF] Soldemo | The surface structure of Cu 2 O (100)[END_REF] Cuprous oxide Cu 2 O crystallizes in a cubic structure in which Cu and O atoms are coordinatively saturated with two and four first neighbors, respectively (hereafter referred to as Cu csa and O csa ). The oxygen ions form a bcc lattice and are surrounded by a tetrahedron of copper ions. Along the (111) orientation, the structure may be seen as a stacking of neutral and dipolefree O-4Cu-O trilayer repeat units. The stoichiometric O-terminated surface (referred to as ST in Figure 1) exposes an equal number of under-coordinated oxygen O cus and copper Cu cus ions, characterized by the loss of one first neighbor compared to the bulk crystal. The O cus ions are positioned in the corners of the Cu 6 O 6 rings, and the Cu cus ions, bonded to a sub-surface oxygen ion, are located in the ring centers. A Cu cus -deficient termination (CuD in Figure 1) obtained by removing all copper Cu cus ions, also displays a (1 × 1) periodicity.

Both (1 × 1) and ( √ 3 × √ 3)R30 • surface periodicities have been experimentally recorded along this orientation, which were initially attributed to the stoichiometric and an oxygendeficient configuration, respectively. 1 Later on, the (1 x 1)-Cu 2 O(111) surface has been explored by STM on bulk samples 6 and ultra-thin metal-supported films. [START_REF] Jensen | Two new oxygen induced reconstructions on Cu (111)[END_REF][START_REF] Sträter | Growth and surface properties of cuprous oxide films on Au (111)[END_REF][START_REF] Ly | Growing ultrathin Cu 2 O films on highly crystalline Cu (111): a closer inspection from microscopy and theory[END_REF] Its hexagonal pattern exhibits symmetry and spacing consistent with both Cu cus and O cus atoms, which does not allow to discriminate between the stoichiometric and the Cu-deficient termination. [START_REF] Soon | Surface oxides of the oxygencopper system: Precursors to the bulk oxide phase?[END_REF][START_REF] Islam | First principles investigation on the stabilization mechanisms of the polar copper terminated Cu 2 O (111) surface[END_REF][START_REF] Li | Stability and electronic properties of the O-terminated Cu 2 O (111) surfaces: Firstprinciples investigation[END_REF][START_REF] Nilius | Incorrect DFT-GGA predictions of the stability of non-stoichiometric/polar dielectric surfaces: The case of Cu 2 O (111)[END_REF] The reconstructed ( √ 3 × √ 3)R30 • surface, which is the subject of the present work, arouse much interest in the community. Relying on LEED 1 and STM, [START_REF] Zhang | Atomistic determination of the surface structure of Cu 2 O (111): experiment and theory[END_REF][START_REF] Önsten | Atomic structure of Cu 2 O (111)[END_REF] it has been interpreted as a configuration in which one third of the surface O cus atoms are missing, whether in the absence or presence of the unsaturated Cu cus atoms (configurations CuD-O s or ST-O s in Fig. 1, respectively). The former was predicted to be the most stable in the whole range of attainable oxygen chemical potentials in DFT+U approach. [START_REF] Bendavid | Firstprinciples predictions of the structure, stability, and photocatalytic potential of Cu 2 O surfaces[END_REF] However, despite the importance of this issue, the actual atomic structure of this surface reconstruction remains ambiguous, since many other surface configurations combining added or removed copper or oxygen atoms on the surface or sub-surface could be compatible with the ( √ 3 × √ 3)R30 • periodicity. Moreover, GGA approaches (whether pure GGA or GGA+U) without explicit inclusion of exchange terms have already been proved to provide wrong predictions on polar surface stability in Cu 2 O, [START_REF] Nilius | Incorrect DFT-GGA predictions of the stability of non-stoichiometric/polar dielectric surfaces: The case of Cu 2 O (111)[END_REF] so that simulations based on more advanced exchange-correlation functionals are needed to determine the configurations of lowest energy.

For these reasons, the present work is devised as a combination of high-resolution STM measurements on the Cu 2 O(111) (

√ 3× √
3)R30 • reconstructed surface, and simulations based on a DFT-HSE approach. The experimental results irrevocably invalidate the vacancy models by disclosing the location and orientation of the adatoms responsible for shamrock-shaped STM contrast. Among an extensive set of surface configurations, numerical simulations reveal a new configuration consistent with all the topographic characteristics derived from our high resolution STM images and being thermodynamically the most stable in the whole range of oxygen chemical potentials. This nanopy-ramidal model for the ( √ 3 × √ 3)R30 • surface reconstruction involves Cu 4 O entities formed by three copper adatoms located at the center of the Cu-O rings on-top a Cu cus atom and capped by an oxygen atom. The structural and electronic characteristics of the nanopyramidsatom coordination, absence of dangling bonds, charge state -strongly differ from those present in the vacancy models, and should have an incontestable impact on the surface reactivity.

Methods

Experimental Methods

The experiments have been carried out in an ultrahigh-vacuum chamber (p ≈ 2 * 10 -10 mbar) equipped with a liquid-nitrogen cooled scanning tunneling microscope (STM) and common surface-science tools for sample preparation and analysis. The Cu 2 O films were prepared by depositing ≈10 ML of high-purity copper on an Au(111) single crystal, followed by a highpressure oxidation in 50 mbar oxygen and 450 K. An oxidation pressure in the mbar regime turned out to be crucial to prepare oxide films thick enough to develop bulk-like properties. Low-pressure films, on the other hand, exhibit terminations that are unknown for bulk crystals, as discussed in Refs. [START_REF] Jensen | Two new oxygen induced reconstructions on Cu (111)[END_REF][START_REF] Sträter | Growth and surface properties of cuprous oxide films on Au (111)[END_REF][START_REF] Nilius | Incorrect DFT-GGA predictions of the stability of non-stoichiometric/polar dielectric surfaces: The case of Cu 2 O (111)[END_REF] To promote longrange ordering, the samples have been postannealed in front of an oxygen nozzle (p ≈ 10 -4 mbar) at about 650 K for 20 min. The stoichiometry of so prepared oxide layers was determined with X-ray photoelectron spectroscopy to be compatible with cuprous oxide (Cu 2 O). Also, the sharp, hexagonal spot pattern seen in low-energy electron diffraction (LEED) was found to match the (1×1) lattice of Cu 2 O(111). In addition, a ( √ 3 × √ 3)R30 • superstructure became discernable in LEED, pointing to a surface reconstruction (Fig. 2, inset). A similar reconstruction was earlier found on the (111) surface of bulk Cu 2 O crystals. 1,14

Computational Methods

Calculations have been performed within the DFT framework implemented in VASP. [START_REF] Kresse | Ab initio molecular dynamics for liquid metals[END_REF][START_REF] Kresse | Efficient iterative schemes for ab initio total energy calculations using a plane-wave basis set[END_REF] Beyond standard GGA (PW91), [START_REF] Perdew | Generalized gradient approximation made simple[END_REF] the hybrid (HSE06) [START_REF] Heyd | Hybrid functionals based on a screened Coulomb potential[END_REF][START_REF] Heyd | Hybrid functionals based on a screened Coulomb potential[END_REF] approximation to the exchangecorrelation functional has been used as to improve the description of the Cu 2 O band gap. The calculated HSE gap width (1.96 eV) increases significantly over the GGA one (0.4 eV) and is close to the experimental value of 2.15 eV. [START_REF] Meyer | Binary copper oxide semiconductors: From materials towards devices[END_REF] Let us note that the HSE gap is also much larger than those obtained within GGA+U approximation (0.5-0.8 eV, for U = 2-8 eV, respectively). [START_REF] Bendavid | Firstprinciples predictions of the structure, stability, and photocatalytic potential of Cu 2 O surfaces[END_REF] While HSE significantly improves the calculated electronic structure, it impacts only little the bulk and surface structural characteristics. We find bulk Cu 2 O lattice parameters of 4.31 Å and 4.29 Å in HSE and GGA, respectively, compared to the experimental value of 4.27 Å. Finally, it is worth reminding that an isolated O 2 molecule does not suffer in HSE from the overbinding typical for GGA functionals (calculated O 2 dissociation energy is equal to 4.9 and 6.4 eV in HSE and GGA, respectively).

The interaction of valence electrons with ionic cores is described within the projector augmented wave method. [START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmented-wave method[END_REF] Standard copper, and soft oxygen (energy cutoff of 270 eV) pseudopotentials provided by VASP were used in all bulk and surface calculations, enabling a complete structural relaxation at the hybrid level for all considered systems. GGA results obtained with the soft and the full (energy cutoff of 400 eV) oxygen pseudopotential were confronted and showed very satisfactory agreement (difference of bulk Cu 2 O lattice parameter smaller than 0.01 Å, difference of (1 × 1) surface energies below 0.002 eV/ Å2 ).

The Cu 2 O(111) surfaces are modelled with symmetric slabs composed of four O-4Cu-O trilayers, at the calculated in-plane lattice parameter of bulk Cu 2 O (6.10/6.07 Å in GGA/ HSE, respectively). We use asymmetric slabs with one bare (labeled B in the following) and one fully hydroxylated termination (labeled H) to improve the convergence with respect to the slab thickness. The latter is the stoichiometric Cu are estimated with usual ab-initio thermodynamics expression:

23 σ B (µ O ) + σ H (µ O ) = (E a-slab Cu 2 O -E bulk Cu 2 O -3E H 2 O -∆N O ( 1 2 E gas O 2 +∆µ O ))
S where E a-slab Cu 2 O is the energy of the asymmetric slab, ∆N O is the number of excess oxygen atoms with respect to the Cu 2 O stoichiometry, and E H 2 O is the total energy of a water molecule. S is the surface area of the ( √ 3 × √ 3)R30 • unit cell. Since the evaluation of σ B requires the knowledge of σ H , the total energy of a symmetric slab with both terminations equivalently hydroxylated with 3 water molecules each was also calculated.

In all calculations, slabs are separated by at least 10 Å of vacuum and dipole correction is applied to slabs with two non-equivalent terminations. For the most relevant configurations, we have additionally checked that increasing the vacuum thickness up to 15 Å changes the total slab energy by less than 0.01 eV. Such supercells with more vacuum were actually used in simulations of STM images. The Brillouin zone of the ( [START_REF] Monkhorst | Special points for Brillouin-zone integrations[END_REF] Atomic coordinates of all ions in the oxide film are relaxed until residual forces dropped below 0.01 eV/ Å.

√ 3 × √ 3)R30 • surface unit cell is sampled with a (2 × 2 × 1) Monkhorst-Pack k- point grid.
Simulations of STM images rely on the Tersoff-Hamann approximation. We have determined the computational settings which provide the best matching of the simulated and experimental STM images of the Cu-deficient termination. As shown in the SI Figure S1, the conventional value of 2.5*10 -8 e/ Å3 for the isodensity surface provides the most satisfactory matching for both positive and negative biases and is systematically used.

Experimental Results

STM topographic images of the Cu 2 O(111) surface reveal wide, flat terraces that are entirely covered by nanometer-sized protrusions of ≈1.5 Å height (Fig. 2a,b). The protrusions, being referred to as shamrocks in the following, are arranged on an incomplete hexagonal lattice with 10.5 Å periodicity. This distribution is readily assigned to the ( √ 3 × √ 3)R30 • superstructure seen in LEED (Fig. 2 inset), although only 40-50% of the lattice sites are actually filled with the shamrock units. In larger openings of the ( resolved. Two distinct contrast schemes are distinguished, both showing the porous nature of the bulk-cut Cu 2 O(111) surface, as composed of Cu-O six rings (Fig. 3a,b). In the first scheme, each surface pore is delimited by three atomsized triangles comprising the Cu ions in the Cu 6 O 6 rings. The respective imaging scheme has thus been labelled "copper contrast" and occurs mostly at negative sample bias. The oxygen contrast, on the other hand, reveals six protrusions on the Cu 6 O 6 rings that are readily assigned to the surface O ions. Hereby, chemically saturated (O csa ) and unsaturated (O cus ) O ions alternate along the ring, with the O cus species appearing slightly brighter in the STM. In both contrast schemes, the center of Cu-O six rings remains dark, suggesting that the Cu cus ions that would be expected at those positions are not available in the surface. This conclusion is consistent with the characteristics of simulated STM images of the stoichiometric and Cu-deficient (1 × 1) surfaces (see Supporting Information (SI), Figure S1).

√ 3 × √ 3)R30 • reconstruction, the fundamental (1 × 1) lattice of Cu 2 O(111) is
By correlating the atomically-resolved STM images of (1 × 1) patches with a model of bulk-cut Cu 2 O(111), the topographic position of the shamrock units can be determined. The shamrocks exclusively locate in the Cu-O ring centers, while corner or edge positions are not populated. The neighboring rings to a shamrock site always remain empty, which gives rise to the characteristic ( √ 3 × √ 3)R30 • superstructure seen in reciprocal space. Typical shamrock features protrude the Cu 2 O surface by about 1.5 Å with slightly smaller values at positive compared to negative bias. A bias-independent shamrock height points to a topographic instead of an electronic contrast regime. The mean height of 1.5 Å hereby matches an extra layer of atoms sitting on top of the ideal (1 × 1) plane. The shamrocks typically exhibit a trifold symmetry that can sometimes be resolved into three discernable maxima (Fig. 3c). The mean distance between the shamrock's symmetry centers has been determined to 3.5 Å, although their exact STM appearance is governed by the tip shape (Fig. 3c). The shamrock units arrange into small domains containing not more than 100 individual species. Neighboring domains are separated by boundaries, along which more than one Cu-O six ring remains unoccupied.

The only other surface irregularities besides the shamrocks are surface O vacancies, whereby only the O cus ions are affected. The density of this defect type is very low, typically smaller than 0.04 nm -2 . In contrast to the shamrock features, the O cus defects exhibit a pronounced bias-dependent contrast (Fig. 4a). They appear bright in topographic images taken above 2.0 V, but develop a negative contrast at low positive and negative bias. Whereas the O cus vacancies are difficult to be identified in topographic data, they leave a clear fingerprint in the STM conductance maps (Fig. 4b). Especially at 2.3 V, they stick out due to their large conductance values, but merge with the background signal of the oxide film at lower bias.

Theoretical Results

To construct an atomistic model consistent with these experimental observations, we have considered a large set of surface configurations compatible with the ( √ 3 × √ 3)R30 • unit cell by introducing both ad-and missing species in surface and sub-surface trilayers of the stoichiometric and Cu-deficient terminations. The surface energies of the key low energy configurations discussed in the following are plotted as a function of oxygen chemical potential in Fig. 5, while a complete set is shown in the SI (Fig. S2).

Before entering into the details of the phase diagram, let us note the substantially different sequences of relative stability predicted by HSE and GGA. Not only is the predicted most stable configuration not the same, but also the relative stability of the simple stoichiometric and Cu-deficient terminations as well as the position of the bulk Cu 2 O stability limit visibly differ. As pointed out in Ref. [START_REF] Nilius | Incorrect DFT-GGA predictions of the stability of non-stoichiometric/polar dielectric surfaces: The case of Cu 2 O (111)[END_REF] hybrid calculations tend to better reproduce the experimental value of the bulk Cu 2 O stability limit and, due to a better account of the Cu 2 O(111) band gap and surface ionization potential, they avoid the serious underestimation of surface energy of the Cu-deficient termination which is otherwise over-stabilized by GGA.

Figure 5 proves that both the ST-O s or CuD-O s vacancy models proposed in the literature are thermodynamically unstable in both GGA and HSE approaches, due to the large cost of oxygen vacancy formation at the Cu 2 O(111) surface. Indeed, filling the surface vacancy with an oxygen atom transforms them into the (1×1) stoichiometric (ST) or Cu-deficient (CuD) terminations, respectively, and brings a very substantial gain in stability. Interestingly, the energy cost to create surface oxygen vacancies is systematically about 0.9 eV/vacancy larger in the HSE approximation than in GGA, in line with larger HSE oxygen vacancy formation energy in bulk Cu 2 O. Moreover, the comparison between X-O s and X-O ss energies in Figure 5 (with X=ST or CuD) shows that moving the oxygen vacancy from the surface to the subsurface site directly below the surface Cu cus cation According to the HSE results, the stability enhancement of the PY configuration with respect to the Cu-deficient termination originates from the reduced number of under-coordinated O cus atoms (from three to one per unit cell) and from the reduced formal surface charge excess at these polar terminations. [START_REF] Goniakowski | Polarity of oxide surfaces and nanostructures[END_REF] While we are not aware of an existing theoretical prediction of such configuration on the Cu 2 O(111) surface, it bears a topological similarity to the (3,3;1,0) reconstruction proposed for the (100) orientation, [START_REF] Soldemo | The surface structure of Cu 2 O (100)[END_REF] and with the Cu 3 O entities seen during the growth of Cu 2 O on Pt(111). [START_REF] Gloystein | Copper Oxidation on Pt (111)-More than a Surface Oxide?[END_REF]

Discussion

Comparison Between Experimental and Simulated STM Images

In order to substantiate the relevance of the nanopyramid model as the likely structure of the ( √ 3 × √ 3)R30 • reconstruction besides stability considerations, we now confront simulated STM images of the pyramidal (PY) and the two vacancy (CuD-O s and ST-O s ) models, Figure 7, with experimental STM images (Figures 2 and 3).

As expected, the shamrock pattern produced by the PY model is correctly centered in the middle of the surface hexagonal rings and, as in the experimental images, its lobes point toward the neighboring O cus atoms. The three Cu adatoms which protrude above the surface O cus atomic plane produce three distinct lobes at negative bias, with a bright spot superimposed at their center, which corresponds to the capping oxygen atom. At positive bias, the shamrock pattern is less distinct. The simulated STM contrast is dominated by the capping oxygen atoms and approaches closely the single-spot image observed experimentally for some tip shapes and/or bias values (Fig. 3c). Thus, the PY model of the ( √ 3 × √ 3)R30 • reconstruction not only agrees with the experimental observation by its position and orientation on the surface but also by a weakly bias dependent STM contrast principally driven by topographic effects.

Conversely, the vacancy models produce a STM contrast centered at one corner of the surface rings and strongly dependent on bias. This is very clear in the ST-O s configuration, where, similarly to the stoichiometric surface, the positive bias contrast is dominated by the contri- bution of empty s state of Cu cus (see SI, Figure S3). The resulting spot distance of 4.7 Å is close to the geometrical distance between the Cu cus ions (4.9 Å). Conversely, negative bias contrast is rather due to the filled d z2 states at the same Cu cus atoms, which dramatically reduces the separation of the three maxima.

Electronic effects are even more striking in the STM image of the CuD-O s configuration, where the overall weak contrast resembles closely that of the pristine Cu-deficient surface. Due to the proximity of Cu csa atoms around the vacancy (2.5 Å), their empty s states (Figure S3) do not produce the three clearly separated spots, characteristic of the observed shamrock reconstruction (Fig. 3c). The O cus vacancy thus rather appears as a single spot in the simulated STM image. At negative bias, the nearly in-plane orientation of the Cu-O bonds makes the occupied Cu d z2 orbitals of the three terminal Cu atoms point towards the vacancy rather than out of the surface plane, making the three atoms hardly visible. The contrast reversal between positive and negative polarity agrees with experimental STM images of the few observed oxygen vacancies (Figure 4).

In summary, the STM signatures of the CuD-O s and ST-O s vacancy models of the ( √ 3 × √

3)R30 • reconstruction are in sound disagreement with the experimental results, regarding the shamrock location, topographic height and orientation with respect to the surface rings, as well as their strongly bias-dependent STM contrast of mainly electronic origin. Conversely, the STM signature of the proposed pyramidal model PY matches very satisfactorily the experiment. Not only is the shamrock structure correctly located and oriented with respect to the Cu 2 O(111) lattice but also its STM contrast is principally driven by topographic effects thus producing a relatively weak bias dependence.

Oxygen Vacancies at the

Cu 2 O(111) Surface with surface oxygen vacancies are unstable with respect to moving the vacancy sub-surface. For example, while the formation of an oxygen vacancy on the top of the Cu 4 O entity in the PY configuration requires a high energy, much less energy (reduction of 0.6/0.7 eV/vacancy in HSE/GGA) is needed if the vacancy is instead created in the subsurface site just below the Cu cus cation. Such strong preference for sub-surface vacancy location at the Cu 2 O(111) surface is shared by no other sub-surface or bulk site. It is quite unusual and, to our knowledge, it has not been previously recognized. Indeed, in a large variety of materials of very different nature, vacancy formation energies grow with the number of broken bonds, which gives a strong bias for vacancies to segregate at surfaces. In the present case, although four O-Cu bonds are broken in the sub-surface site compared to three at the Cu 2 O(111) surface, our results reveal the existence of an electronic mechanism, quasi-independent of the HSE/GGA approximation, which opposes and overrides the usual coordination effect. It is related to the different localization of the two electrons left by the removed oxygen at the surface and sub-surface sites. In the former case, they are evenly distributed over the 3 nearest-neighbor Cu cations, resulting in their partial metallization and fractional oxidation state. In contrast, in the case of a sub-surface vacancy below a Cu cus atom, this Cu cus cation participates to the electron accommodation, to the point of becoming practically neutral, and the charge modification of the three other neighbors is much smaller. As a consequence, despite its charge neutrality and the absence of O-Cu bonds, this Cu cus atom remains strongly trapped in the center of the surface ring. Indeed, when surface Cu cus cations are absent, for example in the CuD-O s configuration, the energy gain of moving the O vacancy to sub-surface is significantly reduced (Figure 5), which clearly highlights the peculiar role of surface Cu cus ions in the unusual behavior of vacancies at the Cu 2 O(111) surface. Such electronic stabilization presents some similarity with the electron localization around subsurface oxygen vacancies found in defective ce-ria films. [START_REF] Jerratsch | Electron localization in defective ceria films: A study with scanning-tunneling microscopy and density-functional theory[END_REF]

Conclusion

By a combination of high-resolution STM measurements and numerical simulations based on a DFT-HSE approach, we challenge the previous interpretations of the Cu 2 O(111) ( √ 3 × √

3)R30 • surface reconstruction. We show that the existing models do not meet the STM observations, that is shamrock features located at the Cu-O ring centers, orientation of the spots towards the surface O cus atoms and topographic nature of the STM contrast. Moreover, they are thermodynamically unstable, due to the preference, not recognized before, of oxygen vacancies to locate sub-surface rather than in the surface plane. We propose a new model of the ( √ 3 × √ 3)R30 • surface reconstruction in which the shamrock protrusions seen in STM are formed by nanopyramids made of three copper adatoms located at the center of the Cu-O rings above the Cu cus atoms and capped by an oxygen atom.

In a broader perspective, we find that the relative stability pattern of the various surface terminations may be rationalized by two main effects: (i) a strong thermodynamic bias for oxygen vacancy filling which might be one cause for the p-type conductance of Cu 2 O and (ii) a clear preference for sub-surface rather than surface oxygen vacancies. While the first effect is much more pronounced in HSE, in line with larger HSE bulk oxygen vacancy formation energy, the second one relies on the enhanced capacity of the surface Cu cus atoms to accommodate an electron excess and is independent of the level of approximation. These effects are responsible for the thermodynamic instability of the vacancy models and result in an enhanced stability of the pyramidal ( √ 3 × √ 3)R30 • reconstruction, being even more stable than the known (1 × 1) configurations.

A Cu 2 O(111) surface model that is not based on low-coordinated Cu cus ions anymore will drastically change common interpretation schemes of the surface reactivity of cuprous oxide.

Acknowledgments

Financial support from the DFG grant Ni 650-5 Towards an atomic-scale understanding of ideal, defective and doped cuprous oxide surfaces and interfaces is gratefully acknowledged.

Supporting Information

The Supporting Information is available free of charge at XXX.

It gives additional information on theoretical results, concerning the stability of several Cu

Figure 1 :

 1 Figure 1: Atomic structure of the stoichiometric (ST) and Cu-deficient (CuD) (1 × 1) Cu 2 O(111) surface and the two vacancy models of ( √ 3 × √ 3)R30 • reconstruction, built by removing 1/3 O cus atoms from either the stoichiometric (ST-O s ) or Cu-deficient (CuD-O s ) surfaces. (1 × 1) and ( √ 3 × √ 3)R30 • unit cells and color codes for atoms are indicated. The three copper atoms Cu vac in the direct vicinity of the oxygen vacancies are represented with a dark blue color.

2 :

 2 2 O(111) termination with 3 dissociated water molecules per ( √ 3 × √ 3)R30 • surface unit cell. Total surface energies of asymmetric slabs as a function of oxygen chemical potential ∆µ O (a) Overview (60 × 60 nm 2 ) and (b) close-up (10.5 × 10.5 nm 2 ) STM topographic image of a high-pressure Cu 2 O film on Au(111) (U B = -1.1 V, I = 50 pA). The insets in (a) and (b) show a LEED pattern (E kin = 37.5 eV) and a Fast-Fourier transformation of the oxide surface, respectively, and depict the fundamental (1 × 1) lattice and ( √ 3 × √ 3)R30 • superstructure spots due to the protruding shamrock units.

Figure 3 :

 3 Figure 3: STM topography of a pristine Cu 2 O (1 × 1) region with a few shamrock protrusions at the edges (3 × 3 nm 2 ). While individual Cu ions are resolved in (a) (U B = -1.1 V), the O ions are detected in (b) (U B = +2.0 V). With this insight, a structure model may be derived for Cu 2 O (1 × 1), in which the shamrock units occupy the centers of Cu-O six rings. (c) STM appearance of the shamrock species for characteristic tip states (3 × 3 nm 2 , U B = +2.0 V, I = 50 pA).

Figure 4 :

 4 Figure 4: STM images (left) and corresponding differential conductance maps (right) of a Cu 2 O(111) surface with shamrock maxima, measured with oxygen contrast (6×6 nm 2 , I = 50 pA). The two encircled O cus vacancies in the lattice undergo a characteristic contrast reversal for bias voltages inside the conduction band (bright) and the oxide band gap (dark).

Figure 5 :

 5 Figure 5: Calculated HSE (left) and GGA(right) surface energies of various ( √ 3 × √ 3)R30 • model reconstructions of the Cu 2 O(111) surface, plotted as a function of oxygen chemical potential. The labels "-O s " and "-O ss " refer to configurations with surface and subsurface oxygen vacancies, respectively (see text).

Figure 6 :

 6 Figure 6: Atomistic representation of the pyramidal reconstruction (PY) of the ( √ 3 × √ 3)R30 • Cu 2 O(111) surface. Left: top view of the ( √ 3 × √ 3)R30 • unit cell. Right: Top and side views of a nanopyramid at the center of a Cu 6 O 6 ring. Same color code as in Fig. 1.

Figure 7 :

 7 Figure 7: From top to bottom, simulated constant current STM images of the pyramidal PY, vacancy CuD-O s and vacancy ST-O s models of the Cu 2 O(111) ( √ 3 × √ 3)R30 • reconstruction, at negative and positive biases (left and right of the atomistic representations, respectively).

Figure 5

 5 Figure 5 reveals an unexpected peculiarity of surface oxygen vacancies. Indeed, whether starting from the pyramidal, the stoichiometric or the Cu-deficient structures, all configurations

3 )

 3 R30 • configurations, their simulated STM images and their Local Densities of States (LDOS) (pdf).

Graphical TOC Entry