Guillaume Chapuis
email: gchapuis@lanl.gov

Hristo Djidjev
email: hristodjidjevdjidjev@lanl.gov

Georg Hahn
email: g.hahn@lancaster.ac.uk

Guillaume Rizk
email: guillaume.rizk@inria.fr

Finding Maximum Cliques on the D-Wave Quantum Annealer

Keywords: Maximum clique, Quantum annealing, D-Wave 2X, Optimization, Gurobi

This paper assesses the performance of the D-Wave 2X (DW) quantum annealer for finding a maximum clique in a graph, one of the most fundamental and important NP-hard problems. Because the size of the largest graphs DW can directly solve is quite small (usually around 45 vertices), we also consider decomposition algorithms intended for larger graphs and analyze their performance. For smaller graphs that fit DW, we provide formulations of the maximum clique problem as a quadratic unconstrained binary optimization (QUBO) problem, which is one of the two input types (together with the Ising model) acceptable by the machine, and compare several quantum implementations to current classical algorithms such as simulated annealing, Gurobi, and third-party clique finding heuristics. We further estimate the contributions of the quantum phase of the quantum annealer and the classical post-processing phase typically used to enhance each solution returned by DW. We demonstrate that on random graphs that fit DW, no quantum speedup can be observed compared with the classical algorithms. On the other hand, for instances specifically designed to fit well the DW qubit interconnection network, we observe substantial speed-ups in computing time over classical approaches.

Introduction

The emergence of the first commercially available quantum computers by D-Wave Systems, Inc. [START_REF] Denchev | What is the computational value of finite-range tunneling?[END_REF] has provided researchers with a new tool to tackle NP-hard problems for which presently, no classical polynomial-time algorithms are known to exist and which can hence only be solved approximately (with the exception of very small instances which can be solved exactly).

One such computer is D-Wave 2X, which we denote here as DW. It has roughly 1000 units storing quantum information, called qubits, which are implemented via a series of superconducting loops on the DW chip. Each loop encodes both a 0 and 1 (or, alternatively, -1 and +1) value at the same time through two superimposed currents in both clockwise and counter-clockwise directions until the annealing process has been completed and the system turns classical [START_REF] Bunyk | Architectural considerations in the design of a superconducting quantum annealing processor[END_REF][START_REF] Johnson | Quantum annealing with manufactured spins[END_REF].

The device is designed to minimize an unconstrained objective function consisting of a sum of linear and quadratic binary contributions, weighted by given constants. Specifically, it aims at minimizing the Hamiltonian H = H (x 1 , . . . , x N ) = i∈V a i x i + (i,j )∈E a ij x i x j [START_REF] Balas | Finding a maximum clique in an arbitrary graph[END_REF] with variables x i ∈ {0, 1} and coefficients a i , a ij ∈ R, where V = {1, . . . , N} and E = V × V [START_REF] King | Benchmarking a quantum annealing processor with the time-totarget metric[END_REF]. This type of problem is known as a quadratic unconstrained binary optimization (QUBO) problem. When the coefficients a i and a ij are encoded as capacities of the couplers (links) connecting the qubits, H describes the quantum energy of the system: During annealing, the quantum system consisting of the qubits and couplers tries to settle in its stable state, which is one of a minimum energy, i.e., of a minimum value of H . In order to solve a given optimization problem, one has to encode it as a minimization problem of a Hamiltonian of type [START_REF] Balas | Finding a maximum clique in an arbitrary graph[END_REF].

Similarly to the random moves considered in a simulated annealing classical algorithm, a quantum annealer uses quantum tunneling to escape local minima and to find a lowenergy configuration of a physical system (e.g., constructed from an optimization problem). Its use of quantum superposition of 0 and 1 qubit values enables a quantum computer to consider and manipulate all combinations of variable values simultaneously, while its use of quantum tunneling allows it to avoid hill climbing, thus giving it a potential advantage over a classical computer. However, it is unclear if this potential is realized by the current quantum computing technology, and by the DW computer in particular, and whether DW provides any quantum advantage over the best available classical algorithms [11,[START_REF] Rønnow | Defining and detecting quantum speedup[END_REF].

This article tries to answer these questions for the problem of finding a maximum clique (MC) in a graph, an important NP-hard problem with multiple applications including network analysis, bioinformatics, and computational chemistry. Given an undirected graph G = (V , E), a clique is a subset S of the vertices forming a complete subgraph, meaning that any two vertices of S are connected by an edge in G. The clique size is the number of vertices in S, and the maximum clique problem is to find a clique with a maximum number of vertices in G [START_REF] Balas | Finding a maximum clique in an arbitrary graph[END_REF].

We will consider formulations of MC as a QUBO problem and study its implementations on DW using different tools and strategies. We will compare these implementations to several classical algorithms on different graphs and try to determine whether DW offers any quantum advantage, observed as a speedup over classical approaches.

The article is organized as follows. Section 2 starts with an overview of related work that aims to solve graph and combinatorial problems with DW, in particular previous work on the maximum clique problem. Section 3 proceeds by introducing the qubit architecture on the DW chip as well as available software tools. We also describe a QUBO formulation of MC together with its implementations on DW and present methods for dealing with graphs of sizes too large to fit onto the DW chip. Section 4 presents an experimental analysis of the quantum software tools and a comparison with several classical algorithms, both for graphs small enough to fit DW directly as well as for larger graphs for which decomposition approaches are needed. We conclude with a discussion of our results in Section 5.

In the rest of the paper, we denote a graph as G = (V , E), where V = {1, . . . , N} is a set of N ∈ N vertices and E is a set of undirected edges.

Related work

Several publications available in the literature aim at searching for a quantum advantage within a variety of problem classes. Existing publications often target a particular (NP-complete) problem and compare the performance of a quantum annealer (by D-Wave Systems, Inc. [START_REF] Denchev | What is the computational value of finite-range tunneling?[END_REF]) to state-of-the-art classical or heuristic solvers. Early examples include multiple query optimization in databases, analyzed by [START_REF] Trummer | Multiple Query Optimization on the D-Wave 2X adiabatic Quantum Computer[END_REF], who investigate scaling behavior and show a speed-up of several orders of magnitude over classical optimization algorithms, and [START_REF] Cao | Solving Set Cover with Pairs Problem using Quantum Annealing[END_REF], who in contrast do not detect any quantum speedup for the set cover with pairs problem, one of Karp's 21 NP-complete problems. Other work include graph partitioning via quantum annealing on DW in the context of QMD (quantum molecular dynamics) applications [START_REF] Mniszewski | Graph Partitioning using the D-Wave for Electronic Structure Problems[END_REF][START_REF] Ushijima-Mwesigwa | Graph Partitioning using Quantum Annealing on the D-Wave System[END_REF], in which graph partitioning is shown to reduce the computational complexity of QMD. Test sets for integer optimization are investigated in [START_REF] Coffrin | Challenges and successes of solving binary quadratic programming benchmarks on the dw2x qpu[END_REF], who observe an advantage of DW over Gurobi [START_REF]Gurobi optimizer reference manual[END_REF] both in terms of speed and quality of solution.

In [START_REF] Stollenwerk | Discrete optimisation problems on an adiabatic quantum computer[END_REF], a general introduction to the DW architecture and the representation of problem instances in Ising and QUBO format is given as well as a QUBO formulation for the maximum clique problem. However, the authors do not actually report any computation results for finding maximum cliques on DW, nor do they compare DW to state-of-the-art heuristic solvers. In contrast to [START_REF] Stollenwerk | Discrete optimisation problems on an adiabatic quantum computer[END_REF], we solve the maximum clique problem on DW for a variety of test graphs and compare its solutions to the ones of state-of-the-art classical solvers. Moreover, we present a graph splitting algorithm allowing to solve problem sizes larger than those embeddable on DW, analyze its scaling behavior, and investigate the influence of alternative QUBO formulations on the solution quality.

In [START_REF] Boothby | Fast clique minor generation in Chimera qubit connectivity graphs[END_REF], the authors consider finding large clique minors in the DW hardware Chimera graph C(m, n, l), defined as the m × n grid of K l,l complete bipartite graphs (also called unit cells). The authors present a polynomial time algorithm for finding clique minors in the special case of the Chimera graph only. Such clique minors are needed to embed problem instances of arbitrary connectivity onto the current and future Chimera architectures, given the problem size is not larger than the clique minor. In contrast, in the present article we consider finding maximum cliques in arbitrary graphs with DW by minimizing a QUBO for the maximum clique problem (applied to the user-specified arbitrary input graph). This step requires an embedding of our QUBO onto DW's Chimera graph, for which the algorithm of [START_REF] Boothby | Fast clique minor generation in Chimera qubit connectivity graphs[END_REF] can be beneficial. However, the work of [START_REF] Boothby | Fast clique minor generation in Chimera qubit connectivity graphs[END_REF] does not substitute for the embedding and minimization of a QUBO when finding cliques in arbitrary graphs as considered in our work.

Instead of attempting a full solution via DW, other publications propose using a quantum annealer to assist in finding a solution of certain problem classes, which often are of a practical and thus more complex nature. For instance, [START_REF] Dridi | Homology computation of large point clouds using quantum annealing[END_REF] consider computing the (algebraic) homology of a data point cloud and propose to reduce this computation to a minimum clique covering, which is then suggested to be solved using DW. No empirical results are presented. A real-world application (the network scheduling problem) is considered in [START_REF] Wang | Quantum versus simulated annealing in wireless interference network optimization[END_REF]. The authors demonstrate an advantage of quantum over simulated annealing; moreover, they show how to obtain more admissible solutions with DW by introducing an additional weight into the QUBO that incrases the gap between linear and quadratic QUBO terms. In [START_REF] Nguyen | Solving sparse representations for object classification using the quantum d-wave 2x machine[END_REF], a sparse coding model is trained using samples obtained via DW from a Hamiltonian with L p sparseness penalty. A graph flow problem in real-world traffic network analysis is considered in [START_REF] Neukart | Traffic flow optimization using a quantum annealer[END_REF], who employ gps coordinates of cars in Beijing as training data.

Another class of publications is concerned with the theory of quantum annealing and the problem of benchmarking quantum computations. For instance, [START_REF] Rogers | Ising Simulations on the D-Wave QPU[END_REF] empirically verify the known phase transition in magnetization for the 2D Ising model with DW. In [START_REF] Thulasidasan | Generative Modeling for Machine Learning on the D-Wave[END_REF], the author proposes to run Markov Chain Monte Carlo using samples generated by DW from a suitable Boltzmann distribution. The question whether random spin-glass problems are a suitable type of problem to detect a quantum advantage over classical approaches is considered in [START_REF] Perdomo-Ortiz | On the readiness of quantum optimization machines for industrial applications[END_REF], who also study the problem of benchmarking quantum annealing vs. classical CMOS computation.

Solving MC on D-Wave

This section introduces the DW chip architecture and briefly presents three tools provided by D-Wave Inc. to submit quadratic programs to the quantum computer.

We also introduce the QUBO formulation of MC needed to submit an MC instance to DW. The section concludes with an algorithmic framework designed to solve instances of MC which are not embeddable on DW.

DW hardware and software

The Qubit architecture

DW operates on roughly 1000 qubits. The precise number of available qubits varies from machine to machine (even of the same type) due to manufacturing errors which render some of the qubits inoperative. The qubits are connected using a specific type of network called Chimera graph, C 12,12,4 (see Fig. 1), comprised of a lattice of 12 × 12 cells, where each cell is a 4 × 4 complete bipartite graph. DW can naturally solve Ising and QUBO problems where non-zero quadratic terms are represented by an edge in the Chimera graph.

The particular architecture of the qubits implies two important consequences: First, the chip design actually only allows for direct pairwise interactions between two qubits which are physically adjacent on the chip. For pairwise interactions between qubits not physically connected, a minor embedding of the graph describing the non-zero structure of the Hamiltonian matrix into the Chimera type graph is needed, which maps a logical variable into one or several physical qubits on the chip. Minor embeddings are hence necessary to ensure arbitrary connectivity of the logical variables in the QUBO. The largest complete graph that the DW can embed in theory has 1+4•12 = 49 vertices. In practice, the largest embeddable graph is slightly smaller (n ≈ 45) due to missing qubits arising in the manufacturing stage.

When more than one qubit is used to represent a variable, that set of qubits is called a chain. The existence of chains has two vital consequences, which will play an important role in the analyses of Section 4. On the one hand, the need for chains uses up qubits, which would otherwise be available to represent more variables in the quadratic program, thereby reducing the maximum problem sizes that can directly be solved on DW. This is the reason for the relatively small sizes of N = 45 for QUBO problems (1) that fit onto DW when the corresponding Hamiltonians are dense (contain nearly all quadratic terms), despite the fact that more than 1000 qubits are available in DW.

On the other hand, due to the imperfections of the quantum annealing process caused by environmental noise, limited precision, and other shortcomings, solutions returned by D-Wave do not always correspond to the minimum energy configuration. In the case of chains, all qubits in a chain encode the same variable in (1) and hence should have the same value, but for the reasons outlined above this may not be the case. This phenomenon is called a broken chain, and it is not clear which value should be assigned to a variable if its chain is broken. Clearly, chains can be ensured to not break by increasing their coupler weights, but as we will see in the next section this may significantly reduce the accuracy of the solver.

D-Wave solvers

D-Wave Inc. provides several tools that help users submit their QUBO problems to the quantum processor, perform the annealing, apply necessary pre-and post-processing steps, and format the output. This section briefly describes several such tools used in this article.

Sapi Sapi stands for Solver API and provides the highest level of control one can have over the quantum annealer. It allows the user to compute minor embeddings for a given Ising or QUBO problem, to choose the number of annealing cycles, or to specify the type of post-processing. Sapi interfaces for the programming languages C and Python are available.

One can also use a pre-computed embedding of a complete 45-vertex graph, thus avoiding the need to run the slow embedding algorithm.

QBsolv QBsolv is a tool that can solve problems in QUBO format which are of a size that cannot natively fit onto DW. Larger problems (with more variables or more connections than can be mapped onto the corresponding Chimera graph) are analyzed by a hybrid algorithm, which identifies a small number of significant rows and columns of the Hamiltonian. It then defines a QUBO on that subset of variables which fits DW, solves it, and extends the found solution to a solution of the original problem.

QSage In contrast to Sapi or QBsolv, QSage is a blackbox hybrid solver which does not require a QUBO or Ising formulation as input. Instead, QSage is able to minimize any function operating on a binary input string of arbitrary size. For this it uses a tabu search algorithm enhanced with DW-generated low-energy samples near the current local minimum. To ensure that also input sizes larger than the DW architecture can be processed, QSage optimizes over random substrings of the input bits.

QUBO formulations of MC

Recall that a QUBO problem can be written as minimize

x i ∈{0,1} H = 1≤i<j ≤N a ij x i x j , ( 2
)
where the weights a ij , i = j , are the quadratic terms and a ii are the linear terms (since x 2 i = x i for x i ∈ {0, 1}). There are multiple ways to formulate the MC problem as a QUBO. One of the simplest is based on the equivalence between MC and the maximum independent set problem.

An independent set S of a graph H is a set of vertices with the property that for any two vertices v, w ∈ S, v and w are not connected by an edge in H . It is easy to see that an independent set of H = (V , E) defines a clique in graph G = (V , E), where E is the complement of set E. Therefore, looking for the maximum clique in G is equivalent to finding the maximum independent set in H . The corresponding constraint formulation for MC is maximize

x i ∈{0,1} N i=1 x i subject to (i,j )∈E x i x j = 0, (3)
where G = (V , E) is the input graph and E is the complement of E. The equivalent unconstrained (QUBO) minimization of (3), written in the form [START_REF] Batagelj | An o(m) algorithm for cores decomposition of networks[END_REF], is

H = -A N i=1 x i + B (i,j )∈E x i x j , (4)
where one can determine that the coefficients/penalties A and B can be chosen as A = 1, B = 2 (see [START_REF] Lucas | Ising formulations of many np problems[END_REF]). A disadvantage of the formulation ( 4) is that H contains an order of N 2 quadratic terms even for sparse graphs G, which limits the size problems for which MC can be directly solved on DW.

Solving larger MC instances

To solve the MC problem on an arbitrary graph, we develop several algorithms that reduce the size of the input graph by removing vertices and edges that do not belong to a maximum clique and/or split the input graph into smaller subgraphs of at most 45 vertices, the maximal size of a complete graph embeddable on DW. Let G(V , E) be a connected graph of n vertices.

Extracting the k -core

The

k-core of a graph G = (V , E) is the maximal subgraph of G whose vertices have degrees at least k. It is easy to see that if G has a clique C of size k + 1, then C is also a clique of the k-core of G (since all vertices in a k-clique have degrees k-1)
. Therefore, finding a maximum clique of size no more than k + 1 in the original graph G is equivalent to finding such a clique in the k-core of G (which might be a graph of much smaller size).

One can compute the k-core iteratively by picking a vertex v of degree less than k, removing v and its adjacent edges, updating the degrees of the remaining vertices, and repeating while such a vertex v exists. The algorithm can be implemented in optimal O(|E|) time [START_REF] Batagelj | An o(m) algorithm for cores decomposition of networks[END_REF].

We also apply another reduction approach, which we refer to as edge k-core, to reduce the size of an input graph g using a known lower bound lower bound on the clique size. This approach combining k-core and edge k-core is given in pseudo-code notation as Algorithm 1. In Algorithm 1, we first aim to reduce the size of g by simply extracting its k-core, where k is set to the currently known lower bound. It is easily shown that for two vertices v, w in a clique of size c, the intersection of the two neighbor lists of v and w has size at least c -1. We therefore choose a random vertex v in sg and remove all edges (v, e) satisfying |N(v) ∩ N(e)| < lower bound -2 (here N(v) denotes the set of neighbor vertices of v), as such edges cannot be part of a clique with size larger than lower bound. Since this changes the graph structure, we attempt to extract the lower bound-core at the end again before returning the reduced graph.

Algorithm 1 Graph reducing -core based algorithm

Graph partitioning

This divide-and-conquer approach aims at dividing G into smaller subgraphs, solves the MC problem in each of these subgraphs, and combines the subproblem solutions into a solution of the original problem. If one uses standard (edgecut) graph partitioning, which divides the vertices of the graph into a number of roughly equal parts so that the number of cut edges, or edges with endpoints in different parts, is minimized, then the third step, combining the subgraph solutions, will be computationally very expensive. Instead, we will use CH-partitioning, recently introduced in [START_REF] Djidjev | Graph partitioning methods for fast parallel quantum molecular dynamics[END_REF].

In CH-partitioning, there are two levels of dividing the vertices of G into subsets. In the core partitioning, the set V of vertices is divided into nonempty core sets C 1 , . . . , C s such that i C i = V and C i ∩ C j = ∅ for i = j . There is one halo set H i of vertices for each core set C i , defined as the set of neighbor vertices of C i that are not from C i . Recall that a vertex w is a neighbor of a vertex v iff there is an edge between v and w. We define the cost of the CH-partitioning

P = ({C i }, {H i }) as cost(P) = max 1≤i≤s (|C i | + |H i |). ( 5 )
The

i C i = V
, where V is the set of the vertices of G, then v belongs to some core C j .

We will next show that for any vertex w = v from K, w ∈ C j ∪ H j , which will imply that all vertices of K are in C j ∪ H j , implying the correctness of the proposition.

If w ∈ C j then the claim follows.

Assume that w ∈ C j . We will show that w ∈ H j . Since K is a clique, there is an edge between any two vertices from it, and hence there is an edge between v and w. Since, by definition, H j consists of all neighbors of vertices from C j that are not in C j , v ∈ C j , w ∈ C j , and w is a neighbor of v, then w ∈ H j .

Using Proposition 1, the solutions to all subproblems of a CH-partitioning can be combined into a solution of the original problem at an additional cost of only O(s) = O(n), where s is the number of the sets of the partition.

One may conjecture that increasing s in (5) will always reduce the cost, but this is not always the case. If the minimum cost is achieved for s = 1, or if some of the parts of the partition are still too large, then the method in the next subsection might be applied.

Vertex splitting

This method is similar to a special case of the previous one, obtained by choosing s = 2, letting C 1 contain only a single vertex v, and letting C 2 contain all other vertices V \ {v}. Moreover, while the halo H 1 of C 1 is defined as above, we set C 1 = ∅ and H 2 = ∅. As a result, G is divided into two subgraphs, G 1 containing all neighbors of v without v itself, and G 2 containing all vertices of G except v, see Fig. 2. Because this partitioning is uniquely determined by a single vertex, we call it a vertex-splitting partitioning. The cost of such a partitioning is again given by [START_REF] Bunyk | Architectural considerations in the design of a superconducting quantum annealing processor[END_REF].

Proposition 2 Given a vertex-splitting partitioning of G, ({C 1 , C 2 }, {H 1 , H 2 = ∅}), the size of the maximum clique of G is equal to max{k 1 + 1, k 2 }, where k i , i = 1, 2, is the size of a maximum clique of the subgraph of G induced by C i ∪ H i .
Proof Let v be the vertex that defines the partition. If there is a maximum clique of G that contains v let K be such a clique, otherwise let K be any maximum cliques of G. Consider the following two cases.

Case 1: v belongs to K. Then the set of the vertices of K consists of v and a subset V 1 of vertices from G 1 . Moreover, since there is an edge between any two vertices of K, there is an edge between any two vertices of

V 1 , which means that V 1 defines a clique K 1 in G 1 .
Assume that K 1 is not a maximum clique of G 1 , i.e., there exists a clique K 1 in G 1 with more vertices than K. Then adding v to the vertices of K 1 will result in a clique in G of size larger than K, which is a contradiction to the choice of K. Hence K 1 is a maximum clique in G 1 , whose size was denoted by k 1 . Since K consists of v and the vertices of K Since H 2 = ∅, vertex splitting can be used in cases where CH-partitioning fails. Moreover, if there is a vertex of degree less than n -1, this method will always create subproblems of size smaller than the original one. However, the total number of subproblems resulting from the repeated use of this method can be too large. A more efficient algorithm can be obtained if all the above methods are combined.

Combining the three methods

We use the following algorithm to decompose a given input graph G into smaller MC instances fitting the DW size limit. We assume that the size k + 1 of the maximum clique is known. (Otherwise, use the procedure of this section in a binary-tree search fashion to determine the size of the maximum clique. This increases the running time by a factor O(log k) = O(log n) only.) We also have an implementation that, instead of "guessing" the exact value of k, uses lower bounds on k determined by the size of the largest clique found so far.

The algorithm works in two phases. First, we apply the kcore algorithm on the input graph and then CH-partitioning on the resulting k-core.

Consequently, we keep a list L of subgraphs (ordered by their number of vertices), which is initialized with the output of the CH-partitioning step. In each iteration and until all produced subgraphs fit the (DW) size limit, we choose a vertex v from the largest subgraph sg, extract the subgraph ssg induced by v and its neighbors and remove v from sg. The k-cores of the two subgraphs produced at this iteration are then inserted into L. Second, we compute the maximum clique on DW for any subgraph in L of size small enough.

Algorithm 2 gives the pseudo-code of this approach. It returns a list of subgraphs of an input graph g sorted in increasing order of their number of vertices as well as an updated lower bound on the maximum clique size. The parameters of Algorithm 2 are the input graph g, a maximal number of vertices vertex limit for which the maximum clique problem is solved directly on a subgraph, and a lower bound on the clique size found so far (lower bound). All returned subgraphs have the property that their size is at most vertex limit. Since the algorithm attemps to solve MC exactly on graphs not larger than vertex limit, the parameter vertex limit in our case can be set to the maximal number of vertices embeddable on DW. Algorithm 2 works as follows. First, a sorted list of graphs called subgraphs (sorted in descending order of the degree of the subgraphs) is created and initialized with g. As long as the largest subgraph (denoted as subgraphs[-1] in Python notation) has at least vertex limit nodes, the current largest graph sg in the list (command pop()) is returned, sg is removed from list subgraphs, and a vertex v is chosen according to some rule specified in the function choose vertex (see the end of Section 3.3 for possible approaches). Then, the induced subgraph ssg by vertex v is extracted and deleted from sg.

A graph reduction step via the function reduce graph is then applied to sg which reduces the size of the graph using the currently known best lower bound lower bound on the clique size. The graph reduction is given separately as Algorithm 1.

Suppose sg still contains vertices after reduction. If the degree of sg after reduction is less than vertex limit, we attempt to solve the MC problem exactly on sg using some function solve() (for instance via DW) and update lower bound. Otherwise, sg is inserted again into the list subgraphs.

The same step is repeated for the subgraph ssg with the exception that ssg does not have to be re-inserted into list subgraphs at the end. This is because the subgraph induced by a single vertex v either contains a clique or can be removed.

Removing a vertex v in line 5 of Algorithm 2 decreases the size of sg by one in each iteration, thus the algorithm terminates in finite time once all generated subgraphs have size at most vertex limit.

Lastly, we describe our procedure choose vertex(sg) for choosing the next vertex to be removed from sg. A vertex with high degree will potentially greatly reduce the size of sg, however at the expense of also producing a large subgraph ssg.

In order to maximize the impact of removing a vertex, we successively try out three choices: a vertex of highest degree, a vertex of median degree and, if necessary, a vertex of lowest degree in sg. If the vertex of lowest degree has degree |V | -1, then sg is a clique: In this case, solving MC on sg can be omitted and lower bound can be updated immediately.

Experimental analysis

The aim of this section is to investigate if a quantum advantage for the MC problem can be detected for certain classes of input graphs. To this end, we compare the DW solvers of Section 3.1.2 to classical ones on various graph instances -from random small graphs that fit the DW chip to (larger) graphs tailored to perfectly fit DW's Chimera architecture. We also evaluate our graph splitting routine of Section 3.3 on large MC instances. First we briefly describe classical solvers that will be used in the comparison.

Classical solvers

Apart from the tools provided by D-Wave Inc., we employ classical solvers in our comparison, consisting of: A simulated annealing algorithm working on the Ising problem (SA-Ising), a simulated annealing algorithm specifically designed to solve the clique problem (SA-clique, see [START_REF] Geng | A simple simulated annealing algorithm for the maximum clique problem[END_REF]), softwares designed to find cliques in heuristic or exact mode (the Fast Max-Clique Finder fmc, see [START_REF] Pattabiraman | Fast algorithms for the maximum clique problem on massive sparse graphs[END_REF]), the software tool pmc (see [START_REF] Rossi | A fast parallel maximum clique algorithm for large sparse graphs and temporal strong components[END_REF]), and the Gurobi solver [START_REF]Gurobi optimizer reference manual[END_REF].

SA-Ising

This is a simulated annealing algorithm working on an Ising problem formulation. The initial solution is a random solution, and a single move in the simulated algorithm is the flip of one random bit.

SA-clique

We implemented a simulated annealing algorithm specifically designed to find cliques, as described in [START_REF] Geng | A simple simulated annealing algorithm for the maximum clique problem[END_REF]. As SA-clique only finds cliques of a user-given size m, we need to apply a binary search on top of it to find the maximum clique size. Its main parameter is a value α controlling the geometric temperature update of the annealing in each step (that is, T n+1 = αT n ). A default choice is α = 0.9996. A value closer to 1 will yield a better solution but will increase the computation time.

Fast Max-Clique Finder (fmc, pmc) These two algorithms are designed to efficiently find a maximum clique for a large sparse graph. They provide exact and heuristic search modes. We use version 1.1 of software fmc [START_REF] Pattabiraman | Fast algorithms for the maximum clique problem on massive sparse graphs[END_REF] and pmc (github commit 751e095) [START_REF] Rossi | A fast parallel maximum clique algorithm for large sparse graphs and temporal strong components[END_REF].

Post-processing heuristics alone (PPHa)

The DW pipeline includes a post-processing step: First, if chains exist, a majority vote is applied to fix any broken chains. Then a local search is performed to ensure that any solution is indeed a local minimum (the raw solutions coming from DW might not be in a local minimum, see [START_REF]D-Wave post-processing guide[END_REF]). For a given solution coming out of the pipeline, one might wonder what the relative contributions of DW and of the post-processing step are. For some small and simple problems, the post processing step alone might be able to find a good solution.

We try to answer this issue by solely applying the post-processing step, and by comparing the result with the one obtained by quantum annealing. However, postprocessing by DW runs on the DW server and is not available separately.

To enable us to still use the DW post-processing alone, we employ the following procedure. We set a very high absolute chain strength (e.g., 1000 times greater than the largest weight in our Ising problem), and turn on the auto-scale feature mapping QUBO weights to the interval [-1, 1]. Because of the limited precision of the DW hardware (DW maps all QUBO weigths to 16 discrete values within [-1, 1]), chain weights will be set to the minimum value -1 while all other weights will be scaled down to 0. In this way, the quantum annealer will only satisfy the chains rather than the actual QUBO we are interested in. As chains will not be connected to other chains, and as all linear terms will be zero, each chain will be assigned a random value -1 or +1. Applying the DW post-processing step to such a QUBO with large chain weights will therefore result in the post-processing step being called with a random initial solution. We hence expect to obtain results stemming from the post-processing step only (with random starting point). This method will be referred to as PPHa, post-processing heuristic alone.

Gurobi Gurobi [START_REF]Gurobi optimizer reference manual[END_REF] is a mathematical programming solver for linear programs, mixed-integer linear and quadratic programs, as well as certain quadratic programs. We employ Gurobi to solve given QUBO problems (Ising problems can be solved as well, nevertheless Gurobi explicitly allows to restrict the range of variables to binary inputs, making it particularly suitable for QUBO instances). Instead of solving MC directly with Gurobi, we solve the dual The edge probability used to generate those graphs is given in the first column. Since for such small graphs, every software returned the correct solution, we only report the running times. Gurobi solves the dual problem, leading to reversed graph densities and timings problem, that is we computed a maximum independent set on the complement graph.

Small graphs with no special structure

We generate four random graphs with increasing edge densities for our experiments. We considered edge probabilities ranging from 0.3 to 0.9 in steps of 0.05. We compare the execution times of DW using the Sapi interface and the different solvers listed in Section 3.1.2 to the classical solvers of Section 4.1.

Results are shown in Table 1. For small graphs, every solver returns a maximum clique, therefore the table shows execution times only. We can see that (a) software solvers are much faster than DW, with pmc being the fastest by several order of magnitudes; (b) DW and PPHa exhibit equal results and execution times. This shows that for these small graphs, even the simple software heuristic included in the DW pipeline is capable of solving the MC problem. The similar performance of DW and PPHa therefore makes it impossible to distinguish between the contributions from the post-processing heuristic and the actual quantum annealer; (c) Gurobi finds the best solution as well (for the dual of MC, the maximum independent set problem, thus timings decrease in the last column of Table 1), but since Gurobi is an exact solver, its running time is higher than the one of the other methods. We note that the timings for Gurobi are for finding the best solution -letting Gurobi run further to subsequently prove that a found solution is optimal requires a far longer runtime. Moreover, we have observed in Eq. 4 in Section 3.2 that the QUBO for MC leads to an order of N 2 quadratic terms even for sparse graphs. This in turn typically causes the QUBO matrix to be very dense, making it difficult to embed the QUBO onto the Chimera graph. If embedding the QUBO is indeed possible, then usually at the cost of incurring long chains. This is due to the fact that a dense QUBO necessarily contains a large number of couplers between qubits not adjacent on the DW chip, thus requiring re-routing through chain qubits.

In our experiments we observe that this is a delicate case for the quantum annealer: Using high coupler strengths for the chains results in consistent chains after annealing, but comes at the cost of downscaling the actual QUBO weights, thus leading to meaningless solutions. Lower chain strengths often cause many of the chains to be broken, i.e. the physical qubits constituting the chains have different values. Therefore some processing needs to be applied to obtain valid solutions. The most simple one is a majority vote, however all postprocessing rules offered by DW are merely heuristic ways of assigning final values to the qubits. It is not guaranteed that a weighting scheme exists which preserves the QUBO and prevents chains from breaking at the same time.

As an example, Fig. 3 shows the first four broken chains in a typical DW execution of the MC problem on a 45 vertex graph. The chain for x 2 has more zeros and less ones than the one for x 7 , yet after the DW postprocessing algorithm was applied, the variables got correct values x 2 = 1, x 7 = 0 (with apparently PPHa overwriting the inferior DW solution). Our experiments with randomly assigned values to broken chains (see the discussion for PPHa in Section 4.1) similarly show that accurate solutions obtained for small graphs are often mostly due to the post-processing algorithm rather than the quantum annealing by DW.

Graphs of sizes that fit DW

In Section 4.2, we performed experiments with random graphs that can be embedded onto DW. We observed that Figure 3 The first four broken chains (out of 16) produced by DW on a test 45-vertex graph. The first column shows the name of the variable the chain corresponds to and the third column gives the correct value for that variable.

highly optimized software solvers outperformed DW in of speed. This is due to the fact that the largest random graphs we are sure to embed on DW (around 45 vertices) are still comparably small and can hence be solved efficiently with an optimized heuristic. In order to detect a difference between DW and classical solvers, we need to consider larger graphs. In this section we will analyze the behavior of the quantum annealer on subgraphs of DW's chimera graph, i.e. the largest graph we can embed on the DW architecture.

Chimera-like graphs

Since on small graphs we did not observe any speedup of DW compared to the classical algorithms, we now consider graphs that fit nicely the DW architecture. The largest graph that fits DW is the Chimera graph C, and since formulation (4) uses the complement edges, the largest graph that we can solve MC on is the complement of C. Let G denote the complement of any graph G. Note that the graphs C and C are not interesting for the MC problem since C is bipartite and hence C consists of two disconnected cliques, which makes MC trivial on this graph.

Consider now the graph C 1 obtained by contracting one random edge from C. An edge contraction consists of deleting an edge (v 1 , v 2 ) and merging its endpoints v 1 and v 2 into a new vertex v * . With N 1 and N 2 the set of neighboring vertices of v 1 and v 2 , the neighbors of v * are N 1 ∪ N 2 \ {v 1 , v 2 }. Solving the MC problem on C 1 requires the embedding of the complement of C 1 onto DW, which is C 1 . The natural embedding of C 1 onto C maps v * onto a chain of two vertices and all other vertices of C 1 onto single vertices of C. Moreover, if we add any edge to C 1 , the resulting graph will not be embeddable onto C any more since C 1 already uses all available qubits and edges of C. We can thus say that C 1 is one of the densest graphs of size |V | -1 than can be embedded onto C.

We can generalize the aforementioned construction to m random edge contractions; the resulting graph C m will have |V |m vertices, will be one of the densest graphs of size |V |-m that fits onto C, and the chains of such an embedding will be the paths of contracted edges. This family of graphs C m with 0 < m < 1100 is therefore a good candidate for the best-case scenario for the MC problem: The C m family are large graphs whose QUBOs can be embedded onto C and whose solutions of MC are not trivial.

Experiments

We solve the MC problem on the C m family of graphs using DW's Sapi, PPHa and the SA-Ising software, SA-clique, and fmc.

Figure 4 shows the result. We observe that for graph sizes up to 400, PPHa finds the same result as DW. For these small graphs the problem is likely simple enough to be solved by the post-processing step alone. As expected, the simulated annealing algorithms designed specifically for MC (fmc, pmc) are behaving better than the general SA-Ising algorithm. The fmc software is run in its heuristic mode. The comparatively lower quality results we obtain with fmc could be due to the fact that fmc is designed for large sparse graphs but run here on very dense graphs.

For large graphs (≥ 800 vertices), DW gives the best solution. (Note we do not know if that solution is optimal.)

Speedup

Since SA-clique seems to be the best candidate to compete against DW, and moreover since it is considered the classical analogue of quantum annealing, we choose to compute the DW speedup relatively to SA-clique on the C m graph family.

We employ the following procedure: For each graph size, we run DW with 500 anneals and report the best solution. The DW runtime is the total qpu runtime for 500 anneals (approximately 0.15s). For SA-clique, we start with a low α parameter (i.e., a fast cooling schedule), and gradually increase α until SA-clique finds the same solution as DW. The value of α for which SA-clique finds the same solution as DW gives us the best execution time for SA-clique given the required accuracy. The SA-clique algorithm is run on one CPU core of an Intel E8400 @ 3.00GHz.

Figure 5 shows the speedup for different graph sizes of the C m family. We observe that DW is slower than SAclique for graphs with less than 200 vertices. For larger graphs, DW gets exponentially faster, reaching a speedup of the order of a million for graphs with 1000 vertices. This behavior is not unexpected: For small graphs, optimized software solvers can terminate with runtimes far less than the constant anneal time of DW (see Table 1). The larger the graphs, the more pronounced the advantage of DW is due to the fact that the C m graph instances investigated in this experiment are similar to the topology of DW's native Chimera graph. Further detail is given in the following section.

Overall, our experiments show that for large graphs whose QUBOs can be embedded onto C, DW is able to find very quickly a solution that is very difficult to obtain with classical solvers.

Topology

In summary, the results of Sections 4.2 and 4.3 demonstrate that the closer the topology of a problem is to the native Chimera graph (Fig. 1) of the DW chip, the more pronounced the advantage of DW over classical solvers. Moreover, with an increasing problem size, the problem becomes exponentially more difficult for classic solvers, while it takes the same time to run on DW (as long as it can be embedded onto the hardware). Note however, that the larger problem we can fit on DW (with a fixed number of qubits), the smaller average chain length we get. This means that these experiments benefit DW in the comparison with classical solvers twice: on the one hand, the problem becomes much more difficult for the classical solvers due to larger graphs involved; on the other hand, it becomes somewhat easier for DW because the shorter chains improve the accuracy, thereby biasing the results in favor of DW.

Using decomposition for large graphs

We investigate some properties of the graph splitting routine of Section 3.3 which enables us to solve MC instances larger than the size that fits onto the DW chip. In this section, we always use our graph splitting routine to divide up the input graphs into subgraphs of 45 vertices, the largest (complete) graphs that can be embedded on the DW chip. First, we test our graph splitting routine on random graphs with 500 vertices and an edge probability (edge density) ranging from 0.1 to 0.4 in steps of 0.05. Figure 6 shows the number of generated subgraphs (or equivalently, the number of solver calls) against the edge probability. Each data point is the median value of ten runs, the standard deviation is given as error bars. The number of solver calls seems to follow an exponential trend with respect to the edge probability. Second, we investigate the scaling of our graph splitting routine with an increasing size |V |. Since, with a fixed edge probability, graphs become denser (their vertex degrees increase) as their size goes to infinity, we take an alternative approach and fix the average degree d of each vertex: We then generate graphs of size 3000 to 20,000 (in steps of 500) using edge probability p = d/(|V | -1). This ensures that the average vertex degree stays constant as |V | goes to infinity.

We measure both the time t (in seconds) of the graph splitting alone as well as the number n of problems/subgraphs being solved by DW. According to Table 1 (column for DW's interface Sapi), the time to solve each subgraph on the DW chip is 0.15 seconds, thus leading to an overall time for computing MC of t + 0.15 • n seconds.

Figure 7 shows average timings from 100 runs for three fixed average degrees d ∈ {50, 100, 200}. We observe that if d is relatively large in comparison to |V | (which, in particular, appears to hold for |V | ≈ 5000 and d = 200), the k-core and CH-partitioning algorithms are less effective, while the vertex-splitting routine alone produces too many subgraphs, causing the computing time to get disproportionately high. With increasing the number of vertices, we observe a roughly linear increase of the runtime. As expected, higher average degrees d result in denser graphs and thus higher runtimes.

To demonstrate the applicability of our graph splitting routine outside of random graphs, we apply the graph splitting to families of graphs from the 1993 DIMACS Challenge on Cliques, Coloring and Satisfiabilty [START_REF]Clique, Coloring, and Satisfiability: Second DIMACS Implementation Challenge[END_REF], also used in [START_REF] Boros | Preprocessing of Unconstrained Quadratic Binary Optimization[END_REF]. These are Hamming and c-fat graphs. Both graph families depend on two parameters: the number of vertices n and an additional internal parameter, the Hamming distance d for Hamming graphs and the partition parameter c for c-fat graphs. We use the generation algorithms of [START_REF] Hasselberg | Test Case Generators and Computational Results for the Maximum Clique Problem[END_REF] for both graph families. We also employed g and U graphs, defined in [START_REF] Kim | A Hybrid Genetic Algorithm for the MAX CUT Problem[END_REF] (including their generation mechanism), which have previously been used for graph assessments in [START_REF] Boros | Preprocessing of Unconstrained Quadratic Binary Optimization[END_REF][START_REF] Kim | A Hybrid Genetic Algorithm for the MAX CUT Problem[END_REF].

Table 2 shows results for all four graph families. We see that for the graph parameters used in the aforementioned studies, our graph splitting algorithm finds a maximum clique (mostly) within a fraction of a second. The number of generated subgraphs along the way varies widely, from none or one subgraph for g graphs to almost two hundred for Hamming graphs. Lastly, we aim to assess the performance of future generations of DW systems on our clique finding approach for arbitrary large graphs. Essentially, we turn the previous question around: Instead of assessing the graph splitting for a variety of graphs and a fixed DW system, we now look at the evolution of possible future DW machines with an increasing number of qubits and investigate the number of solver calls needed by the graph splitting algorithm (applied to a fixed realization of a random graph with 500 vertices and edge presence probability 0.3).

First, assuming a similar Chimera topology for future generations of DW systems, doubling the number of available qubits will increase the size of the maximal complete subgraph that can be embedded by a factor of √ 2. The maximal size of an arbitrary graph embeddable on DW is shown in Fig. 8 in red (right y-axis). If we assume that the number of qubits doubles with each new generation, seven generations of DW machines are required in order to be able to directly embed and solve an arbitrary 500 vertex graph.

Second, Fig. 8 (blue data line; left y-axis) shows the evolution of the number of solver calls for future DW systems with an increasing number of qubits. We use the envisaged size of the maximal complete subgraph embeddable on future DW machines to set the lower bound parameter of the graph splitting algorithm. In this experiment we applied the graph splitting algorithm to the fixed graph generated with 500 vertices and edge presence probability of 0.3. Each data point is the median of ten runs. The standard deviation of those ten runs is given with error bars. The number of required solver calls of our graph splitting algorithm rapidly decreases in what seems like an exponential trend.

Conclusion

This article evaluates the performance of the DW quantum annealer on maximum clique, an important NP-hard graph problem. We compared DW's solvers to common classical solvers with the aim of determining if current technology already allows us to observe a quantum advantage for our particular problem. We summarize our findings as follows.

1. The present DW chip capacity of around 1000 qubits poses a significant limitation on the MC problem instances of general form that can be solved directly with DW. For random graphs with no special structure that are small enough to fit onto DW, the returned solution is of comparable quality to the one obtained by classical methods. Nevertheless the highly optimized classical solvers available are usually faster for such small instances. 2. Special instances of large graphs designed to fit DW's chimera architecture can be solved orders of magnitude faster with DW than with any classical solvers. 3. For MC instances that do not fit DW, the proposed decomposition methods offer a way to divide the MC problem into subproblems that fit DW. The solutions of all subproblems can be combined afterwards into an optimal solution of the original problem (assuming DW solves the subproblems optimally, which is usually true, but cannot be guaranteed). Our decomposition methods are highly effective for relatively sparse graphs; however the number of subproblems generated grows exponentially with increasing density. We demonstrate that this issue can be alleviated when/if larger D-Wave machines become available (Fig. 8).

Overall, we conclude that general problem instances that allow to be mapped onto the DW architecture are typically still too small to show a quantum advantage. But quantum annealing may offer a significant speedup for solving the MC problem, if the problem size is at least several hundred, roughly an order of magnitude larger than what it typically is for general problems that fit D-Wave 2X. (2015,2016,2017). His research interests include algorithmic/computational mathematics and statistics, mathematical optimisation, Monte Carlo methods, multiple testing, and information theory.

Guillaume Rizk holds a Ph.D. from the Université de Rennes 1 and a master degree from the École Nationale Supérieure d'Informatique et de Mathématiques Appliquées de Grenoble (ENSIMAG). His research interests include DNA/RNA sequencing, graph data structures for nextgeneration sequencing (NGS), GPU accelerated RNA folding, and global alignment of short DNA sequences. He is one of the founders of "Algorizk" which develops interactive physics simula tion apps for mobile devices.

1 3 5 6 7

 3567 def reduce graph(Graph g, int lower bound): 2 extract k core(g, lower bound) Vertex v = choose random vertex(g) 4 for each vertex n in neighbors(g, v) do Set nv = neighbors(g,v) Set nn = neighbors(g,n) List common neighbors = intersection(nv, nn) 8 if length(common neighbors) lower bound-2 then 9 remove edge(g, v, n) (g, lower bound)

Figure 2

 2 Figure 2 Illustration of the vertex splitting algorithm.

Algorithm 2 2 List

 22 Graph splitting algorithm 1 def split(Graph g, int vertex limit, int lower bound):

Figure 4

 4 Figure 4 Best clique size found by the different solvers, relatively to the DW result, on the C m family of graphs.

Figure 5

 5 Figure 5 Speedup on artificial graphs designed to fit the Chimera topology.

Figure 6

 6 Figure 6 Number of solver calls against edge probability. Log scale on the y-axis.

Figure 7

 7 Figure 7 Time of the graph splitting routine as a function of the graph size.

  CH-partitioning problem is finding a CH-partitioning of G of minimum cost. The next statement shows how CH-partitions can be used for solving MC in larger graphs. Given a CH-partitioning ({C i }, {H i }) of a graph G, the size of the maximum clique of G is equal to max

	Proposition 1

i {k i }, where k i is the size of a maximum clique of the subgraph of G induced by C i ∪ H i .

Proof Let K be a maximum clique of G and let v be any vertex of K. Since, by definition of CH-partitioning,

  1 , its size is k 1 + 1. Moreover, G 2 cannot have a clique larger than K since any clique in G 2 is also a clique in G. Hence, k 2 ≤ |K| = k 1 + 1 and |K| = max{k 1 + 1, k 2 }. Case 2: v does not belong to K. Then, K is entirely contained in G 2 and hence |K| ≤ k 2 . On the other hand, G 2 cannot have a larger clique than |K| since any clique in G 2 is also a clique in G, hence |K| = k 2 . Moreover, by the choice of K, any clique containing v is of size less than K, so |K| > k 1 + 1, and therefore |K| = max{k 1 + 1, k 2 }.

Table 1

 1 Running on 45 vertex random graphs.

	Graph	Max. clique size	Runtime [s]
			Sapi	PPHa	QBsolv	fmc	pmc	SA	Gurobi
	p=0.3	5	0.15	0.15	0.05	8 • 10 -6	3 • 10 -5	0.15	102
	p=0.5	8	0.15	0.15	0.06	3 • 10 -4	5 • 10 -5	0.37	38
	p=0.7	13	0.15	0.15	0.04	0.002	8 • 10 -5	0.19	33
	p=0.9	20	0.15	0.15	0.04	0.135	8 • 10 -5	0.28	2

Table 2

 2

	Graph splitting
	algorithm applied to a variety	Graph family	Vertices	Parameter	Largest clique	No. subgraphs	Runtime [s]
	of graph families (first column)
	including their graph	Hamming	128	1	64	0	0.09
	parameters (number of vertices			2	32	196	6.4
	in second column, internal parameter in third column).			4 6	4 2	20 1	0.1 0 . 0 9
		c-fat	200	1	12	3	0.02
				5	58	1	0.53
			500	1	14	3	0.08
				5	6 4	2	2 . 1
				10	126	0	25.6
		g graph	100	10	1	1	0.01
			200	10	1	1	0.01
			500	10	1	1	0.07
			1000	10	1	1	0.26
			2000	10	1	1	1.1
			5000	10	1	1	6.8
			10000	10	1	1	28.3
		U graph	1000	5	7	6	0.49
				10	10	9	0.49
				20	14	11	0.57
			2000	5	7	7	1.9
				10	11	10	2.0
				20	17	14	2.2
		Largest clique found, number of generated subgraphs and overall runtime in seconds for the splitting is
		reported

  Hristo Djidjev received his Ph.D. from the University of Sofia, Bulgaria. He is currently a Scientist at Los Alamos National Laboratory, USA. Before Los Alamos, he has been with the Bulgarian Academy of Sciences, Rice University, and University of Warwick. His interests are in graph algorithms, discrete optimization, software/hardware codesign, bioinformatics, and quantum computing. Georg Hahn received his Ph.D. from Imperial College London in 2015. Following a one-year EPSRC doctoral prize fellowship at Imperial's statistics department (2016-17), he is currently a senior research associate affiliated with the StatScale project, a joint EPSRC programme grant to develop next generation statistical methods for streaming data. Previously he was a post-doc at the statistics department of Columbia University in New York (2015-16) and a visiting researcher at Los Alamos National Laboratory

Acknowledgments

The authors acknowledge and appreciate the

support provided for this work by the Los Alamos National Laboratory Directed Research and Development Program (LDRD). They would also like to thank Dr Denny Dahl for his help while working on the D-Wave 2X machine.


domain: chim phys

text:
Alexandru Boţan
email: alexandru.botan@gmail.com

Virginie Marry

Benjamin Rotenberg

Pierre Turq

Benoit Noetinger

How Electrostatics Influences Hydrodynamic Boundary Conditions: Poiseuille and Electro-osmostic Flows in Clay Nanopores

Keywords: molecular simulation, water/solid interface, charged surface, Poisson-Boltzmann theory, dynamics

come   L'archive ouverte pluridisciplinaire

Introduction

The transport properties of fluid in clay minerals are of great practical importance. On the one hand, the low permeability of clays makes them attractive barriers for the storage of carbon dioxide or toxic chemical and radioactive waste. On the other hand, the slow fluid migration hinders the recovery of shale gas and oil. The theoretical description of fluid transport in clays is difficult because of the multiscale structure of these porous material with a wide range of pore sizes, usually classified as interlayer (< 2 nm), interparticle (2-50 nm) and intergranular (>50 nm) pores. Moreover, adsorption/desorption processes between clay layers can lead to the swelling or shrinkage of clay particles and thereby affect the permeability of the material. Multiscale approaches, both from the theoretical and experimental points of view, are required to understand the transport of fluids through clays.

Quasi-elastic neutron scattering [1][2][3] and molecular dynamics simulations [4][5][6] have shown that the transport of water confined in interlayer pores is purely diffusive. This motion between clay layers strongly depends on many factors such as the clay type, the nature of the counter-ions, the presence of other molecules, temperature, etc. [7][8][9][10] Transport on the macroscopic scale is usually described by Darcy's law, arising from convective flow in the larger pores coupled to the complex geometry of the pore network; the influence of microscopic details such as the chemical composition is less pronounced.

For interparticle pores of a few nanometers, convective flows are more sensitive to the details of the water/clay interactions. The latter strongly depend on the nature of the clay, in particular on whether the mineral layer is charged or not. From a static point of view, most uncharged clays behave as hydrophobic (e.g. large contact angles) while with charged clays the presence of counterions induces hydrophilic behaviours (e.g. swelling in the presence of water). Studies focussing on the dynamic aspects of this interaction (non-slip or slip) in clays are scarce and one thus usually assumes a correspondence between the hydrophilic/phobic static properties and non-slip/slip in the presence of a flow. Contact angle measurements may depend on the experimental technique and other factors such as surface roughness, fluid viscosity, electrical properties, gas bubbles, pressure, etc. [START_REF] Lauga | Chapter Microfluidics: The No-Slip Boundary Condition[END_REF][START_REF] Cao | [END_REF] Nevertheless, Huang et al. found a quasi-universal relationship between water slippage and the contact angle on model hydrophobic surfaces. 13 The hydrodynamic boundary conditions are determined by a number of factors, including electrostatic ones. Cho et al. have shown experimentally that the slip length is primarily influenced by the dipole moment of the fluid rather than the wettability of the liquid at the surfaces, and reported a decrease in the slip length with increasing dipole moment. 14 As an explanation, the authors suggested that the dipole-dipole interactions result in a "surface lattice structure" in the liquid near the surface. Stein et al. studied the surface charge dependence on the ionic transport in nanochannels via electro-osmotic flow and found that the ionic liquid conductivity increases with increasing surface charge of the channel. 15 Dukhin et al. measured ζ -potentials for alumina-KCl interfaces and showed that an electric double layer can even be observed at uncharged surfaces. 16 Joseph and Aluru went into details of this effect using molecular simulations. 17 For this purpose, two neutral silica slit channels with and without partial charges filled by 1 M KCl were modeled. The presence of an electrical double layer and electro-osmotic flow were observed for both channels. The presence of the ions and the wall partial charges change the interfacial orientation of water from that of water near an uncharged surface affecting the properties of the electric double layer near the interface.

Clays provide an interesting class of systems to investigate the microscopic origin of the nonslip/slip behaviour at a mineral surface. Indeed, they cover ranges of nearly identical systems which may vary by the presence or absence of charge, by the charge distribution inside the solid and by the nature of the counter-ions. The non-slip or slip nature of hydrodynamic boundary conditions in montmorillonite clay nanopores has been studied by molecular simulations for both Poiseuille and electro-osmotic flows. In order to investigate the non-slip/slip behaviour at such surfaces, Marry et al. 18 and Dufrêche et al. 19 have shown using non equilibrium molecular dynamics (NEMD) that the Stokes equation provides a reasonable description for pores with sizes larger than 3 nm but that neglecting slip results in large errors in the overall flow (as high as 27% for a 4.5 nm pore and 15% for 8 nm). 20 The purpose of the present paper is to address the question of how electrostatics influences the non-slip/slip behaviour of clay surfaces. To that end, we investigate by molecular simulations Poiseuille and electro-osmotic flows in 6 nm nanopores between surfaces of two clay minerals, montmorillonite and pyrophyllite, which differ only by the presence in the former of substitutions in the mineral lattice, resulting in a net negative charge and the presence of counter-ions. We consider two force fields with different partial charges and analyze their predictions for the hydrodynamic boundary conditions at the surface of both types of clays. In addition, this analysis allows to discuss the relevance of continuous theories to predict the hydrodynamic and electro-osmotic flows in clays nanopores. Molecular simulations finally provide insight into the microscopic origin of the non-slip/slip behaviour in terms of molecular orientation of surface water molecules and the strength of hydrogen bonds they donate to the surface. 4 . Counterions Na + in the interlayer region compensate the negative charge of clay layers. The simulation box contains two clay layers, forming a 4.5 nm slit pore centered at z = 0, water molecules distributed in the pore and Na + ions, the latter only for montmorillonite. Each clay layer contains 8×4×1 unit cells (totalling 1280 atoms) with dimensions 41.44×35.88×6.54 Å 3 . The number of water molecules was chosen as a result of a prior grand canonical Monte Carlo simulation as presented in our previous study: 20 2437 (resp. 2474) for montmorillonite with CLAYFF (resp. SSFF) and 2437 (resp. 2478) for pyrophyllite with CLAYFF (resp. SSFF). Figure 1 illustrates the simulation box. Periodic boundary conditions are applied only in the x and y directions and empty space was added in the z direction. Several force fields have been proposed to describe clay materials. [21][22][23][24][25][26][27][28][29] They generally (but not always) describe the interactions between atoms by pairwise additive Lennard-Jones and Coulomb potentials:

Systems and methods

Microscopic model

V i j = q i q j 4πε 0 r i j + 4ε i j [ ( σ i j r i j ) 12 - ( σ i j r i j ) 6 ] (1)
where r i j is the distance between sites i and j of different molecules, q i is the partial charge of the site, and σ i j and ε i j are LJ parameters deduced from the conventional Lorentz-Berthelot mixing rules. [START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF] In the present work, we consider two of the most commonly used force field to simulate clays, which are of the form described above but differ by the value of the corresponding parameters: (a) the one proposed by Skipper [START_REF] Skipper | [END_REF] and subsequently adapted to the SPC/E water model by Smith, 25 denoted as SSFF in the following, and (b) CLAYFF, introduced by Cygan et al. 27 The LJ parameters and charges are presented in Table 1. Both force fields predict thermodynamic properties of clays in fair agreement with experiments; 27,[32][33][34][35] for diffusion coefficients better predictions are obtained with CLAYFF. 9,36 In order to study the influence of electrostatics on boundary conditions the clay substrate, we treated the clay as rigid for both force fields. While CLAYFF allows to describe the flexibility of the clay framework, this is not possible with SSFF. For a discussion of the effect of flexibility on boundary conditions, using the CLAYFF force field, we refer the reader to Ref. 37. Water molecules are described using the SPC/E (Simple Point Charge/Extended) model, 38 while Na + cations are modeled using the parameters proposed by Smith and Dang. 39 Table 1: Charges q and Lennard-Jones parameters σ and ε of montmorillonite and pyrophyllite for SSFF and CLAYFF

SSFF CLAYFF specie q i (e) σ i (Å) ε i (kcal/mol) q i (e) σ i (Å) ε i (kcal/

Non equilibrium molecular dynamics simulations

In non equilibrium molecular dynamics, the force acting on a particle results from interparticle interactions and from an external force, e.g. from a gravitational or electrical field. In our simulations, we apply either (a) a constant force F = 0.75 cal/mol/Å, or (b) an electric field E = 0.01 V/Å in the x direction (parallel to the clay layers) to simulate a steady hydrodynamic or electroosmotic flow, respectively. We verified the linear response of the system for all reported results.

All simulations were carried out with the LAMMPS molecular dynamics package 40 in the NV T ensemble. The temperature was fixed at 300 K using a Nose-Hoover thermostat, coupled only to the degrees of freedom of the fluid in the y-direction, along the surface and perpendicular to the fluid flow. The velocities for H 2 O and Na + are collected every 0.1 ps for 6 ns, in bins along z of width 0.1 Å, after an equilibration period of 2 ns.

Continuous theory

The traditional description of hydrodynamic flows and electrokinetic phenomena relies on a continuous description of the system, whereby the solvent is a continuum of viscosity η and dielectric constant ε r and the charged surfaces are infinite walls with a uniform charge density. At steady state, the balance between viscous stresses and other forces is described by Stokes equation

η∆v + F v = 0,
v(z) = v(L hyd /2) - F v 2η [ z 2 - ( L hyd 2 ) 2 ] (2)
The flow profile is then known when v(L hyd /2) is set by the corresponding boundary condition, as will be discussed below.

In the case of electro-osmosis, F v arises from an applied electric field

E = Ee x , i.e. F v = c ion (z)
v(z) = v(L hyd /2) - eE η ∫ z L hyd /2 ∫ z ′ 0 c ion (z ′′ )dz ′′ dz ′ (3)
where v(L hyd /2) is again determined by the boundary conditions. The ionic density profile can be determined either from molecular simulations or from an analytical theory such as Poisson-Boltzmann. The predictions of the latter model, presented in appendix, will be discussed below.

Hydrodynamic boundary conditions from molecular dynamics

In this section, we assess the validity of continuous theories by comparing analytical velocity profiles with those obtained in molecular dynamics simulations. When such a description is appropriate, this comparison further allows to determine the relevant hydrodynamic boundary conditions (non-slip or slip) and provides the parameters associated with the continuous model, such as the viscosity of the confined fluid or the slip length. We first analyze the hydrodynamic flow and determine the corresponding boundary conditions. We then use the same conditions and parameters to model the electro-osmotic flow, as a validation of these features. We finally discuss the relevance of Poisson-Boltzmann theory for the prediction of electro-osmotic flows in the present case.

Hydrodynamics

The velocity profiles in pyrophyllite and montmorillonite are given on Figure 2 for the two force fields CLAYFF and SSFF. In the central region of the pore (between -15 Å and 15 Å where the fluid is far enough from the surfaces not to be disturbed 41 ), the profile is parabolic, as predicted by Eq. 2. The viscosity can be deduced from the curvature, while the value of the velocity allows to characterize the boundary conditions.

When the fluid slips on the surface, the velocity does not vanish at the interface and a slip length b can be defined as: 42,43 b

= - v(L hyd /2) v ′ (L hyd /2) (4)
and determined from the molecular simulation velocity profile if L hyd is known. In this slip case, a natural choice to apply the boundary condition is to indentify ±L hyd /2 with the positions ±z GDS of the Gibbs Dividing Surfaces (GDS), which partition the space between a region of homogeneous fluid and a region without fluid: where ρ bulk and ρ(z) are respectively the fluid density in the bulk, i.e. in the central region of the pore, and that at z. Thus z GDS can be calculated from the density profile of the fluid ρ(z) obtained by simulation. From the profiles of Figure 2, we can conlude that this slip condition applies in the cases of pyrophyllite simulated with CLAYFF and of montmorillonite simulated with SSFF (in agreement with our previous work 20 ). The values of η, L hyd and b obtained from the fits are reported in Table 2. When the fluid sticks to the surfaces, a layer of thickness d is immobilized at the surface. The positions of the shear planes are then further from the surface than the GDS:

∫ ±z GDS 0 [ρ bulk -ρ(z)]dz = ∫ +∞ ±z GDS ρ(z)dz (5)
± L hyd 2 = ±z GDS ∓ d. ( 6
)
In this stick case v(L hyd /2) = 0, so that L hyd and consequently d can be determined from the molecular simulation velocity profiles. As can be seen on Figure 2, the non-slip boundary condition applies in the cases of montmorillonite simulated with CLAYFF and pyrophyllite simulated with SSFF. The corresponding parameters L hyd and d can be found in Table 2. In both cases d is of the order of one molecular diameter. The location of the GDS is also correlated with the hydrodynamic boundary conditions: As can be seen in Table 2, the slippage is more pronounced when the GDS is further from the surface oxygen atoms. This result is in agreement of the work of Wang et al., who computed z GDS for muscovite (hydrophilic surface) and talc (hydrophobic). 44 One can first observe that the fluid viscosities deduced from the velocity profiles are all close to the bulk value for SPC/E water (0.67 cP 45 ). The viscosity measured for montmorillonite is only slightly larger than with pyrophyllite. Moreover, slight variations in the location of the GDS are observed depending on the considered case. The most surprising result is that the two force fields lead to opposite conclusions. In the case of CLAYFF, the fluid slips along the surfaces of pyrophyllite and sticks to those of montmorillonite, in agreement with the hydrophobic behaviour of the former which is uncharged and with the presence of hydrated counter-ion at the surface of the latter. [START_REF] Chesworth | Encyclopedia of Soil Science[END_REF] Simulations with SSFF predict the opposite behaviours. In this case, the value of b found with montmorillonite is in good agreement with the one reported in a smaller pore. 20 Even if the values for b and d are small, the prediction of less slippage with pyrophyllite than with montmorillonite is unphysical. The present study thus supports the idea that CLAYFF should be preferred to SSFF when considering dynamical quantities. However, as will be shown below, the comparison between both force fields provides useful insights into the microscopic origin of the hydrodynamic behaviour of water at clay surfaces.

Electro-osmosis

In order to check that the same boundary conditions apply in a different situation, One can first notice a good agreement between the solution of Stokes equation and the MD results when the ionic density profile is taken from MD and when the boundary conditions for the hydrodynamic flow are used. This further supports the above conclusions on the slip/non-slip behaviour of montmorillonite as predicted by the two force fields and underlines the relevance of the corresponding parameters. Indeed, neglecting the slip or stagnation length deteriorates the prediction of continuous description. Even if these lengths are small, their influence extends far from the surface: As an example, we had shown in 20 for montmorillonite simulated with SSFF that the error on the overall hydrodynamic flux is as high as 15 % even for a large pore of 8 nm. The influence of the boundary conditions is even more important for electro-osmosis because the driving force in that case is very localized near the charged surfaces. 18,19 Finally, as can be seen in Figure 3, PB theory is insufficient to predict the electro-osmotic flow, even if the correct hydrodynamic boundary conditions are included, as a result of the poor description of the ionic density profile, and hence of the force acting on the fluid. Similar conclusions had been reached

with model surfaces and a simple LJ fluid by Joly al. [START_REF] Joly | [END_REF] or water by Huang et al. 48 The discrepancy is expected to increase with surface charge density and cation valency, as in Qiao and Aluru, 49 who even observed a charge reversal at a highly charged surface.

Microscopic origin of the boundary conditions

In this section, we elucidate the microscopic origin of the non-slip or slip behaviour of the clay surfaces by analyzing the change in boundary conditions b and d under different simulation conditions. The only difference between pyrophyllite and montmorillonite, which are expected to be hydrophobic and hydrophilic, respectively, is the presence in montmorillonite of negative charges in the middle of the clay layer, arising from the substitution of Al 3+ by Mg 2+ , and of the counterions (Na + in our case). One can envision two main reasons for the differences of behaviour between the two clays, which both result from the presence of this charge. The net negative charge of the surface and the positive layer of counterions at the surface result in an electric field, experienced by the water molecules in between, which influences their orien-tation, and hence the hydrogen bond network, both between the surface water molecules and the surface and among surface molecules. 6,41,50,51 Both force fields predict more hydrogen bonds with the surface in the case of montmorillonite, compared to pyrophyllite, and a similar number of such bonds for a given clay. However nonslip or slip is observed depending on the force field. Thus the number of hydrogen bonds is not sufficient to determine the hydrodynamic behaviour. In the following, we first analyze the role of the net charge and the resulting electric field by comparing montmorillonite and pyrophyllite. Then we discuss the role of the polarity of the charge distribution inside the clay layer by comparing the two force fields in the case of montmorillonite.

Montmorillonite vs pyrophyllite: role of the average electric field

The net charge of the layers and the counterions in montmorillonite generate an electric field experienced by surface water molecules. In order to investigate how this field influences the boundary conditions independently of the presence of the cations (which may modify the water structure in their solvation shell), we simulate hydrodynamic flows between pyrophyllite surfaces in the presence of an applied electric field perpendicular to the surfaces E = Ee z . The profiles are shown on Figure 4 for both force fields for various values of the electric field E = 0.5 and 1.5 V/Å, together with the profiles obtained previously in the absence of electric field. The presence of E induces a dissymmetry of the system which is visible on the velocity profiles. The latter can still be described by parabolas. From the curvature, we deduce the viscosity of the fluid and observe an increase with increasing E: η = 0.80 cP (resp. 0.72 cP) for E = 0.5 V/Å with SSFF (resp.

CLAYFF) and η = 0.85 cP (resp. 0.88 cP) for E = 1.5 V/Å with SSFF (resp. CLAYFF). Indeed, applying E through the whole fluid changes the nature of the fluid also in the middle of the pore, by modifying the orientation of the water molecules and hence the hydrogen bond network.

However, the most striking result is that the boundary conditions change: for CLAYFF, the fluid sticks to the left surface although it slipped in the absence of electric field; for SSFF, the fluid which stuck to the surfaces now slips on the right surface. These effects also coincide with the appearance (resp. disappearance) of a small peak in the water oxygen atoms density profiles with CLAYFF (resp. SSFF), which corresponds to water molecules above hexagonal cavities, as illustrated on Figure 4. The molecules in this region are approximately horizontal with hydrogen atoms pointing slightly away from the surface. In the absence of an applied field, such molecules are present with SSFF but not with CLAYFF.

When applying E, the water molecules tend to orientate their dipole in the direction of the field. On the left surface, this causes water molecules to turn their dipole away from the surface.

In the case of CLAYFF, this new orientation seems to facilitate a closer approach to the surface and the entrance of molecules above the cavities. In the case of SSFF, this increases the occupancy of these sites. On the contrary, on the right surface, water molecules orientate with their hydrogen atoms closer to the surface than their oxygen. With SSFF, this reversal of orientation of the water molecules above hexagonal cavities causes their exit from this position.

All these observations point to a strong correlation between the presence of water molecules on these sites in the case of SSFF and the unexpected behaviour of the fluid at the surface of pyrophyllite. These molecules pinned down to preferential sites on the surface (by stronger interactions than with CLAYFF) may slow down the fluid by donating hydrogen bonds to other water molecules.

This would also explains why the fluid slips more in montmorillonite when simulated with SSFF.

Indeed, the first layer of water molecules is located between the negatively charged clay surfaces and the first peak of the cation density (see Figure 5). Thus these water molecules experience an electric field directed toward the surface, which corresponds to the situation at the right surface with the applied field. As a consequence, water molecules can orientate with their hydrogen atoms toward the surface and form hydrogen bonds with it, as already observed in montmorillonite. 41 In addition, at this surface the electric field causes the molecules above the cavities to exit and one observes slip for montmorillonite with SSFF. This further suggests that these molecules are at least partly responsible for the non-slip boundary conditions in pyrophyllite with this force field.

The change of boundary conditions on the right surface for CLAYFF is less obvious to explain:

an increase in E results in a wider immobilized layer. However, when looking in details at the den- sity profiles of water oxygen atoms given on Figure 6, one can notice a strong similarity between the densities in pyrophyllite with E=1.5 V/Å and that in montmorillonite, for both force fields, at least close to the surfaces. This can be understood by considering the electric field experienced by surface water molecules in montmorillonite. Assuming first that they are located in a plane capacitor of surface charge density σ = ±0.008071 e/Å 2 (between the negative plane inside the clay layer containing the substitutions and the positive plane containing the counter-ions), one can estimate that the field is approximately σ /ε 0 ∼ 1.45 V/Å, which corresponds well to the above value.

A more accurate estimate can be obtained by integrating Poisson equation with the cation density profiles from molecular simulations. In the considered region (corresponding to the first peak of water oxygen densities), the electrostatic fields are on average 1.47 and 1.42 V/Å for CLAYFF and SSFF, respectively. Thus the molecules at the surface of montmorillonite experience an average electric field similar to those at the right surface of pyrophyllite with E=1.5 V/Å, resulting in similar hydrodynamic velocity profiles. The main conclusion of this section is thus that the average electric field plays a key role in determining the hydrodynamic boundary conditions. The resulting reorientation of surface water molecules largely explains the differences observed between montmorillonite and pyrophyllite.

Although the average electric field created by the charge separation seems to be a sufficient explanation without invoking the disturbance of surface water molecules in the hydration shell of cations, the role of the latter is not trivial. Indeed, their occurrence as fully hydrated ions (outer sphere complexes) in the Na + case results from the competition between solvation and attraction to the surface. For other ions (e.g. Cs + ) and other surfaces (e.g. more highly charged or with substitution in the tetrahedral layer, closer to the surface), inner sphere complexes can be observed.

This would in turn modify the electric field experienced by the surface water molecules, hence the hydrodynamic behaviour.

Comparing force fields: role of local interactions

While the previous analysis allows to understand the difference between pyrophyllite and montmorillonite, it remains to clarify the origin of the differences between CLAYFF and SSFF in montmorillonite. To that end, we finally compare the results obtained with both force fields in that case.

Although the field experienced by the first layer of water is almost the same in both cases, the fluid sticks with CLAYFF and slips with SSFF. In fact, the comparison of the densities represented on Figure 5 reveals that the structure of the fluid is slightly different for the two force fields: The peaks are larger and closer to the surface with CLAYFF. This suggests that the fluid is more tighlty bound to the surface in that case, as a result of the differences existing between the two force fields (partial charges and Lennard-Jones parameters).

One of the reasons for a slowing down of the fluid with CLAYFF compared to SSFF can be found in differences in the strength of the hydrogen bonds formed with the surface. Indeed, the number of hydrogen bonds formed with the surface is almost the same for both force fields: 2.6 and 2.5 per unit cell for CLAYFF and SSFF, respectively. In order to roughly estimate their strength, ) and larger with SSFF (-0.8e). This explanation was already invoked in our previous work on fluoro-hectorite. 9 The fact that water molecules at the surface interact more strongly with the surface than in the bulk with CLAYFF explains at least partially why the fluid is slowed down in that case. It is also illustrates from the dynamic point of view the role of the competition between adhesive (surface-water) and cohesive (water-water) interactions in determining the hydrophilic/hydrophobic behaviour of a surface, as recently analyzed on pyrophyllite, talc and fluorotalc. 52

Conclusion

We investigated hydrodynamic and electro-osmotic flows in clay nanopores using NEMD simulations, to understand the role of electrostatic interactions on the hydrodynamic boundary conditions.

For this purpose, we simulated charged (montmorillonite) and uncharged (pyrophyllite) clays using two commonly used force fields: SSFF and CLAYFF. For the 6 nm pore investigated here, the molecular velocity profile can be well described by the solution of the Stokes equation only if anisms at play at the interface between water and charged surfaces in general, and oxides such as silica in particular. In the latter case, hydrodynamic boundary conditions will also be affected by the rugosity arising from the amourphous nature of this material. 53 where A is a constant which can be determined from the electroneutrality condition. The electrostatic potential ψ(z) and ion concentration c ion (z) can be calculated analytically in the present case. 54 For that purpose, one has to define a distance L el where the boundary conditions ψ ′ (±L el /2) = ±σ /ε 0 ε r apply. This distance coincides with L for point ions, but otherwise differs due to ther finite size. 18,19 Moreover, L el and L hyd defined in the main text are a priori different. The solution then reads:

c ion (z) = 1 2πL B α 2 cos 2 (αz) ( 8
)
where α is defined by α tan(αL el /2) = -2πL B σ /e, with L B = e 2 /4πε 0 ε r k B T the Bjerrum length.

For the comparison with molecular simulations, we use σ = -0.008071 e/Å 2 and the value ε r = 71.5, which is the dielectric constant of bulk SPC/E water. 55 The simulated normalized cation density profiles for the two force fields are represented on Figure 7, together with calculated PB profiles for several L el . The profiles are very sensitive to the molecular structure of the solvent close to the surfaces. For an easier comparison with between MD and PB results, we consider the integral density, which corresponds to the fraction of cations found between -L and z, also reported in Figure 7. One can observe that no value of L el allows to reach a satisfactory agreement. Increasing L el improves the description in the center of the pore but deteriorates that in the vicinity of the surfaces. As a compromise, for the calculations of electro-osmotic flows we use L el = 48 Å for both force fields. pressure gradient electric eld

Graphical TOC Entry

  have used equilibrium molecular dynamics simulations and Kubo relations to compute electro-osmotic flows. The main conclusion was that the Stokes equation correctly predicts these flows only if slip boundary conditions are taken into account. In a recent study of hydrodynamic (Poiseuille) water flow through Na-montmorillonite nanopores, Boţan et al.

Figure 1 :

 1 Figure 1: Snapshot of the simulation box for montmorillonite: pink, Al; green Mg; yellow, Si; red, O; white, H; blue, Na.

  where v is the velocity and F v the volumic force acting on the fluid. In the case of a hydrodynamic (Poiseuille) flow, the external force is a uniform pressure gradient along x. The profile v = v(z)e x for the present case of a slit pore is parabolic. Defining z = 0 as the z position of the center of the pore and ±L hyd /2 the position of the hydrodynamic interfaces, v(z) reads:

  eE, where c ion (z) is the local counter-ion concentration. The solution of Stokes equation now satisfies:

Figure 2 :

 2 Figure 2: Hydrodynamic velocities (F hydro =0.75 cal/mol/Å) as a function of z, for pyrophyllite (black) and montmorillonite (red) simulated with CLAYFF (left) and SSFF (right). The simulation results (dashed lines) are reported together with fits by parabolas (solid lines). The vertical dashed lines indicate the location of the clay surface oxygen atoms.

  we now analyze the electro-osmotic flows in montmorillonite. The ionic density profiles c ion (z) from molecular simulation are introduced in Eq. 3, together with the slip length b or stagnation length d determined in the previous section. The resulting velocity profiles are compared to that obtained in NEMD simulations in Figure 3. For comparison, we also consider the solution of the same equation assuming a priori boundary conditions (b = 0 and d = 0) or using the Poisson-Boltzmann ionic density profile (see the appendix for the determination of the latter).

Figure 3 :

 3 Figure 3: Electro-osmotic velocities as a function of z in montmorillonite. The simulation results (solid line) are compared to the predictions of Eq. 3 for different situations: Using the ionic density profile from molecular dynamics (MD) or from Poisson-Boltzmann theory (PB), including or not the boundary conditions determined from the hydrodynamic flows (slip length b or stagnation length d).

  (a) The hydrated cations close to the surface change the structure and dynamics of water in their hydration shell. Moreover, the presence of the cation may disturb the hydrogen bond network of water in the surface layer. (b)

Figure 4 :

 4 Figure 4: Hydrodynamic flow in pyrophyllite under a force F hydro =0.75 cal/mol/Å in the presence of an applied field E perpendicular to the surfaces. The velocity profiles (left) and water oxygen densities (right) are reported for E=0 (black) 0.5 (red) and 1.5 V/Å (green). The velocity profiles for montmorillonite in the absence of electric field are also indicated (blue). The vertical dashed lines indicate the location of the clay surface oxygen atoms.

Figure 5 :

 5 Figure 5: Density profiles in montmorillonite for the CLAYFF (black) and SSFF (red) force fields: Na + (top), water oxygen atoms (bottom, solid lines) and water hydrogen atoms (bottom, dashed lines). The vertical dashed lines indicate the location of the clay surface oxygen atoms.

Figure 6 :

 6 Figure 6: Density profile for water oxygen atoms in pyrophyllite with CLAYFF (left) and SSFF (right), for various values of applied field E in the direction perpendicular to the surface: 0 (black), 0.5 (red), 1.5 (green) and 2 V/Å (blue). Results for montmorillonite in the absence of electric field are also indicated (dotted lines). The vertical dashed lines indicate the location of the clay surface oxygen atoms.

  (a) proper boundary conditions, with a slip or stagnation length determined from molecular simulation, are taken into account and (b) the ionic density profiles from MD simulations are used in the case of electro-osmotic flow, since the Poisson-Boltzmann equation fails to reproduce the ionic profiles, hence the force acting on the fluid. While the CLAYFF force field predicts a hydrophobic pyrophyllite and hydrophilic montmorillonite as expected from experimental behaviour, the SSFF force field leads to the opposite conclusions.The non-slip or slip boundary conditions at clay surfaces strongly depend on electrostatic interactions of surface water molecules with the surface, both on the presence/absence of a net charge and on the details of the charge distribution inside the clay layer. The main effect of the net charge is the resulting average electric field experienced by surface water molecules between the charged surface and the condensed layer counter-ions, which influences their orientation. The charge distribution inside the clay layer determine the polarity of the surface and hence the strength of hydrogen bonds donated by water molecules to surface oxygen atoms. The competition between hydration and attraction to the surface determines the location of counter-ions at the surface (inner or outer sphere complexes) which should in turn impact the field experienced by surface water and hence hydrodynamic boundary condition. While a systematic study is out of the scope of the present work, we are currently investigating other ions and different clays, e.g. with a larger surface charge or with a charge located in the tetrahedral layer, closer to the fluid. These conclusions should also hold for other types of surfaces and the present study provides insight into the microscopic mech-

Figure 7 :

 7 Figure 7: Cation concentration profiles and their integrals for the two force fields. The MD results (bold black line) are compared to the Poisson-Bolztmann result for L el = 46 Å (red), 48 Å (thin black line) and 51 Å (blue line). The latter case corresponds to L el /2 = z GDS . The vertical dashed lines indicate the location of the clay surface oxygen atoms.

Two different types of clay are considered: pyrophyllite, a dioctahedral 2:1 neutral clay with unit formula

  Al 4 Si 8 O 20 (OH) 4 , and montmorillonite, a cationic clay, which can be obtained from pyrophyllite by partial isomorphous substitution of Al 3+ ions by Mg 2+ ions in the octahedral layers to give a layer composition of Na 0.75 [Si 8 ](Al 3.25 Mg 0.75 )O 20 (OH)

Table 2 :

 2 Hydrodynamic parameters deduced from molecular dynamics simulations

		η (cP)	z GDS (Å) L hyd /2 (Å)	b (Å)	d (Å)
	pyrophyllite CLAYFF	0.65 ± 0.02 24.8 ± 0.2 24.8 ± 0.2 4.1 ± 0.2
	pyrophyllite SSFF	0.63 ± 0.04 25.8 ± 0.2 23.3 ± 0.2		2.5 ± 0.4
	montmorillonite CLAYFF 0.72 ± 0.02 25.5 ± 0.2 23.4 ± 0.2		2.1 ± 0.4
	montmorillonite SSFF	0.68 ± 0.03 25.5 ± 0.2 25.5 ± 0.2 1.5 ± 0.3

  S and 4.69 kcal/mol/Å with O W , indicating that the interactions with the surface are stronger than with the surrounding water molecules. For SSFF, the average forces are 4.26 kcal/mol/Å with O S and 4.69 kcal/mol/Å with O W , i.e. an opposite situation. This is likely a consequence of the fact that the charge of SPC/E oxygen (-0.8476e) is somewhat smaller than that for O S in CLAYFF (-1.05e

we calculate the average force exerted on a water hydrogen atom, either bonded to another water molecule O W or to a surface oxygen atom O S . For CLAYFF, the average forces are 4.75 kcal/mol/Å when bonded to an O

Acknowledgement

The authors acknowledge financial support from the Agence Nationale de la Recherche under grant ANR-09-SYSC-012 and from the Groupement National de Recherches PARIS. The authors thank Jean-François Dufrêche for useful discussions.

Appendix

In this appendix, we discuss the relevance of Poisson-Boltzmann (PB) theory for the description of ionic profiles at the surface of Na-montmorillonite in the absence of added salt. Combining Poisson equation with a Boltzmann distribution of point ions results in the PB equation: