Functions of the Laplacian matrix with application to distributed formation control

Keywords: Combinatorial Laplacian, Network systems, Consensus protocol, Matrix functions

In this paper, we study a class of matrix functions of the combinatorial Laplacian that preserve its structure, i.e. that define matrices which are positive semidefinite, and which have zero row-sum and non-positive off-diagonal entries. This formulation has the merit of presenting different incarnations of the Laplacian matrix appeared in the recent literature, in a unified framework. For the first time, we apply this family of Laplacian functions to consensus theory, and we show that they leave the agreement value unchanged and offer distinctive advantages in terms of performance and design flexibility. The theory is illustrated via worked examples and numerical experiments featuring four representative Laplacian functions in a shape-based distributed formation control strategy for single-integrator robots.

I. INTRODUCTION A. Related work

T HE combinatorial Laplacian is ubiquitous in network sci- ence [START_REF] Newman | Networks[END_REF], and over the last decade we have witnessed the emergence of several variants to describe distributed dynamic processes. In particular, the recent advent of graph signal processing [START_REF] Shuman | Vertex-frequency analysis on graphs[END_REF], [START_REF] Ortega | Graph Signal Processing: Overview, Challenges, and Applications[END_REF], graph neural networks [START_REF] Gama | Graphs, Convolutions, and Neural Networks: From graph filters to graph neural networks[END_REF], [START_REF] Wu | A Comprehensive Survey on Graph Neural Networks[END_REF], and the growing popularity of network systems [START_REF] Mesbahi | Graph Theoretic Methods in Multiagent Networks[END_REF], [START_REF] Bullo | Lectures on Network Systems[END_REF], has had a catalytic effect on the research in this field.

An interesting variant of combinatorial Laplacian is the deformed Laplacian, which has found applications in multiagent systems theory [START_REF] Morbidi | The deformed consensus protocol[END_REF], semi-supervised learning [START_REF] Gong | Deformed Graph Laplacian for Semisupervised Learning[END_REF], and in the design of new centrality measures for undirected and directed networks [START_REF] Grindrod | The Deformed Graph Laplacian and Its Applications to Network Centrality Analysis[END_REF]- [START_REF] Arrigo | Beyond non-backtracking: noncycling network centrality measures[END_REF]. Along the same lines as [START_REF] Morbidi | The deformed consensus protocol[END_REF], the parametric Laplacian has been introduced in [START_REF] Morbidi | The Second-order Parametric Consensus Protocol[END_REF]. The Bethe-Hessian matrix, the reversal of the deformed Laplacian in the case of undirected graphs [START_REF] Arrigo | On the exponential generating function for non-backtracking walks[END_REF], has been used for spectral clustering [START_REF] Saade | Spectral Clustering of Graphs with the Bethe Hessian[END_REF] and community detection in sparse heterogeneous networks [START_REF] Dall'amico | Revisiting the Bethe-Hessian: Improved Community Detection in Sparse Heterogeneous Graphs[END_REF] (see also [START_REF] Stoll | A Literature Survey of Matrix Methods for Data Science[END_REF]). A different definition of deformed Laplacian, which encompasses several Laplacian-like matrices available in the literature (connection, magnetic [START_REF] Fanuel | Magnetic Eigenmaps for the visualization of directed networks[END_REF], signed and dilation), was proposed in [START_REF] Fanuel | Deformed Laplacians and spectral ranking in directed networks[END_REF]. Notably, the dilation Laplacian has been shown to be useful for spectral ranking in directed graphs.

The p-Laplacian, a nonlinear generalization of the combinatorial Laplacian, has recently attracted the attention of the machine learning community. The p-Laplacian reduces to the standard Laplacian for p = 2, and it has been successfully

The author is with the MIS Laboratory, Université de Picardie Jules Verne, 33 rue Saint-Leu, 80039 Amiens, France. Email: fabio.morbidi@u-picardie.fr applied to solve two-class [START_REF] Bühler | Spectral clustering based on the graph p-Laplacian[END_REF] and multi-class [START_REF] Luo | On the eigenvectors of p-Laplacian[END_REF] clustering problems.

The exponential of the negated Laplacian is common in the study of classical transport on lattices and networks, and it is often referred to as the "heat kernel" because of its interpretation as a diffusion process related to the heat equation [START_REF] Michelitsch | Fractional Dynamics on Networks and Lattices[END_REF]. It also appears in discrete-time agreement protocols [START_REF] Mesbahi | Graph Theoretic Methods in Multiagent Networks[END_REF].

Laplacian powers have been widely explored in the literature as well. In particular, the bi-Laplacian (or Laplacian squared) has been considered in [START_REF] Li | Image colorization by using graph bi-Laplacian[END_REF] for image colorization. The bi-Laplacian corresponds to the discretization of the biharmonic PDE equation with Neumann boundary condition, which together with its numerical schemes is widely studied and applied in problems of linear elasticity [START_REF] Henrot | Extremum Problems for Eigenvalues of Elliptic Operators[END_REF]Ch. 11], data interpolation, and computer vision (image inpainting). For non-integer powers in the (0, 1] interval, we obtain the socalled fractional Laplacian, which has been used to study random walks on networks, among other applications [START_REF] Michelitsch | Fractional Dynamics on Networks and Lattices[END_REF], [START_REF] Lischke | What is the fractional Laplacian? A comparative review with new results[END_REF]. In [START_REF] Bautista | L γ -PageRank for Semi-Supervised Learning[END_REF], the authors have proposed a generalization of the PageRank algorithm for semi-supervised learning based on a (non-necessary integer) power of the Laplacian matrix. Finally, it has been recently shown that graph filters can be represented by matrix polynomials of the combinatorial Laplacian [START_REF] Segarra | Optimal Graph-Filter Design and Applications to Distributed Linear Network Operators[END_REF], [START_REF] Levie | CayleyNets: Graph Convolutional Neural Networks With Complex Rational Spectral Filters[END_REF].

B. Original contributions, organization and notation

The original contribution of this paper is twofold. First of all, following [START_REF] Michelitsch | Fractional Dynamics on Networks and Lattices[END_REF], we study matrix functions of the Laplacian of an undirected graph, that preserve its main algebraic properties, i.e. that give rise to matrices which are positive semidefinite and which have zero row-sum and non-positive off-diagonal entries. These Laplacian functions promote the emergence of non-local network correlations: in fact, thanks to their ability to reorganize the information exchanges between the nodes of a graph, long-range interactions are enabled. This general formulation allows us to bring together a number of results scattered across different areas and shed new light on them. Second, we revisit the classical continuous-time consensus protocol by replacing the combinatorial Laplacian with our family of Laplacian functions. These functions leave the agreement value unchanged and offer some practical advantages. In fact, in many situations, by simply tuning a scalar parameter, the user can seamlessly modify the dynamic behavior of the multi-agent system, for example, to guarantee faster convergence towards consensus or to adapt to variable external conditions. This generalizes existing work in the literature, where integer powers of the Laplacian have been considered for fast consensus seeking (by introducing a "multihop relay" interpretation [START_REF] Jin | Multi-Hop Relay Protocols for Fast Consensus Seeking[END_REF]), and for dynamic consensus over wireless sensor networks [START_REF] Manfredi | Design of a multi-hop dynamic consensus algorithm over wireless sensor networks[END_REF]. The price to pay for the increased design flexibility, is that the interaction graph associated with our Laplacian functions is generally more dense than that of the original Laplacian. A simple yet effective approximation of the Laplacian functions is thus proposed to ensure a distributed implementation of the new coordination protocols. The theoretical findings are illustrated via worked examples and validated via numerical simulations in which single-integrator robots run a consensus-based formation control algorithm.

The rest of this paper is organized as follows. In Sect. II, we briefly recall the definition of two matrix functions, and some basic notions of algebraic graph theory. In Sect. III, we study a family of functions of the combinatorial Laplacian that retain its special structure, and in Sect. IV we use them in a novel shape-based distributed formation control protocol, which is validated via extensive numerical experiments. Finally, in Sect. V, the main contributions of the paper are summarized and some possible avenues for future research are outlined.

Notation: Throughout this paper, Z >0 and Z ≥0 denote the sets of positive and non-negative integers, respectively, I n the n×n identity matrix, 0 n the n×n matrix of zeros, 1 a column vector of n ones, and ⊗ the Kronecker product. For a complex number z ∈ C, we use Im(z) to denote its imaginary part, and arg(z) its argument.

II. BACKGROUND MATERIAL

In this section, we briefly recall the definitions of matrix logarithm and matrix pth root, and review some elementary notions of graph theory, which serve as the foundation of the next sections.

A. Matrix logarithm

A logarithm of A ∈ C n×n is any matrix X such that e X = A. Any nonsingular matrix A has infinitely many logarithms [START_REF] Higham | Functions of Matrices: Theory and Computation[END_REF]Th. 1.27]. In this paper, we assume that A ∈ C n×n has no eigenvalues on the closed negative real axis, and log(•) always denotes the principal logarithm, which is the unique logarithm all of whose eigenvalues lie in the strip {z ∈ C : -π < Im(z) < π} [START_REF] Higham | Functions of Matrices: Theory and Computation[END_REF]Th. 1.31]. If A is real, then its principal logarithm is real.

B. Matrix pth root

Let p ≥ 2 be an integer. Matrix X is a pth root of A ∈ C n×n if X p = A. We recall the following result on existence of pth roots [START_REF] Higham | Functions of Matrices: Theory and Computation[END_REF]Th. 7.2].

Theorem 1 (Principal pth root) Let A ∈ C n×n have no eigenvalues on the closed negative real axis. There is a unique pth root X of A all of whose eigenvalues lie in the sector {z ∈ C : -π/p < arg(z) < π/p}, and it is a primary matrix function 1 of A. We refer to X as the principal pth root of A and write X = A 1/p . If A is real, then A 1/p is real.

For p = 2, we have the following extension of Theorem 1 which allows A to be singular (see Problem 1.27 in [START_REF] Higham | Functions of Matrices: Theory and Computation[END_REF]).

Proposition 1 (Square root of a singular matrix) Let A ∈ C n×n have no eigenvalues on the closed negative real axis, except possibly for a semisimple 2 zero eigenvalue. There is a unique square root X of A that is a primary matrix function of A and whose nonzero eigenvalues lie in the open right half-plane. If A is real, then X is real.

There exist four main numerical algorithms for the computation of the matrix pth root: Schur, Newton, and Schur-Newton algorithms, and matrix sign method [START_REF] Higham | Functions of Matrices: Theory and Computation[END_REF]Sect. 7]. The Schur algorithm assumes that A ∈ R n×n is nonsingular, whereas the other methods assume that A has no eigenvalues on the closed negative real axis.

C. Graph theory

Let G = (V, E) be a connected graph, where V = {1, . . . , n} is the set of nodes and E is the set of edges. In this paper, all G's are simple, i.e. they are unweighted, undirected graphs containing no self-loops or multiple edges [START_REF] Godsil | Algebraic Graph Theory[END_REF].

Definition 1 (Adjacency and Laplacian matrix)

• The adjacency matrix A of graph G is an n × n matrix defined as

[A] ij = 1 if {i, j} ∈ E and [A] ij = 0 otherwise. • The Laplacian matrix of graph G is an n × n symmetric positive semidefinite matrix defined as L(G) = D -A where D = diag(d 1 , d 2 , . . . , d n ) = diag(A1) is the degree matrix.
It is easy to verify that L1 = 0, i.e. the row-sums of L are equal to zero, and that the off-diagonal entries of L are nonpositive. Since L is a symmetric matrix, we have

Lv i = λ i v i , i ∈ {1,

III. FUNCTIONS OF THE LAPLACIAN MATRIX

In this section, we study matrix functions f (L) of the combinatorial Laplacian L. These functions define matrices which incorporate the network topology information and which promote non-local interactions between the nodes of graph G. We will exploit them in Sect. IV, to design a new distributed formation control strategy for single-integrator robots. Our exposition here follows the outline of [22, Ch. 1]: however, several new examples and results are provided in Sect. III-A and Sect. III-B, and the notation has been modified to conform to that conventionally used in systems theory. For the sake of simplicity, we henceforth use f ij (L) to denote the (i, j) entry of matrix f (L).

Given a well-defined function f (x) : R → R, the matrix f (L) can be obtained using the series expansion f (x) = ∞ =0 c x where c is the real coefficient of 2 An eigenvalue is called semisimple if its algebraic multiplicity equals its geometric multiplicity.

the th term, or using the spectral decomposition of L. The latter option gives

f (L) = n i = 1 f (λ i ) v i v T i . (1)
Equation [START_REF] Newman | Networks[END_REF] shows that to find f (L), one should calculate the spectrum {λ 1 , λ 2 , . . . , λ n } of L and then compute {f (λ 1 ), f(λ 2 ), . . . , f(λ n )}: the eigenvectors f (L) remain the same as those of L. From (1), we also notice that matrix f (L) is symmetric by construction. For further details on the theory of functions of matrices, see [START_REF] Higham | Functions of Matrices: Theory and Computation[END_REF] and [START_REF] Lancaster | The Theory of Matrices: with Applications[END_REF]Ch. 9]. Note that while equation ( 1) allows to calculate general functions of L, we are interested here in functions which preserve the special structure of the Laplacian matrix (as reported in Sect. II-C). In order to retain these desirable properties, the matrix f (L) should satisfy the following three conditions:

• Condition I: f (L) is positive semidefinite, i.e. the eigenvalues of f (L) are restricted to be positive or zero.

• Condition II: the entries f ij (L) satisfy n j=1 f ij (L) = 0 for i ∈ {1, 2, . . . , n}, or equivalently f (L)1 = 0, i.e.
each row-sum is equal to zero.

• Condition III: the off-diagonal entries of f (L) are nonpositive and they are not allowed to be all simultaneously zero. Hence, owing to Condition II, the diagonal elements of f (L) are all strictly positive.

Condition I is guaranteed if f (x) ≥ 0 for x ≥ 0. Then f (λ i ) ≥ 0 for all i, and the eigenvalues of f (L) can be arranged in increasing order as those of L. On the other hand, by using equation ( 1), it is easy to verify that Condition II is fulfilled if f (0) = 0. However, these two conditions on the function f do not guarantee that the off-diagonal entries of f (L) are non-positive, as required by Condition III.

For example, if we take G = S n (the star graph with n nodes), we have that

L(G) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ n -1 -1 -1 • • • -1 -1 1 0 • • • 0 -1 0 1 • • • 0 . . . . . . . . . . . . . . . -1 0 0 • • • 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . The function f (x) = x 2 is non-negative for x ≥ 0 and f (0) = 0. Thereby, f (L) = L 2 , a.k.a. bi-Laplacian [23],
satisfies Conditions I and II, but the structure imposed by Condition III is not preserved: in fact, L 2 has strictly positive (unitary) off-diagonal entries. On the contrary, the function

f (x) = log(x + 1) is non-negative for x ≥ 0, f (0) = 0, and
f (L) = log(L + I n ) has negative off-diagonal entries. Hence, log(L + I n ) is an admissible Laplacian function.
The following notion of completely monotonic function [START_REF] Miller | Completely monotonic functions[END_REF] plays an important role in the characterization of the admissible functions f (L).

Definition 2 (Completely monotonic function

) A function g(x) defined on 0 < x < ∞ is said to be completely monotonic, if it possesses derivatives g (m) (x) = d m dx m g(x) for all m ∈ Z ≥0 and if (-1) m g (m) (x) ≥ 0 for all x > 0.
In [START_REF] Michelitsch | Fractional Dynamics on Networks and Lattices[END_REF]Sect. 1.4.2], the authors have proved that the functions f (L) that satisfy the necessary Condition III, can be constructed via the scalar functions

f (x) = H(0) -H(-x), 0 ≤ x < ∞, (2)
where H(x) denotes a primitive of the auxiliary scalar function h(x), i.e. h(x) = d dx H(x). The function f (x) can be expressed in terms of a function g(x) = h(-x) defined on 0 ≤ x < ∞, which satisfies the following conditions:

d dx f (x) = g(x) > 0, 0 ≤ x < ∞, (3)
and

(-1) m g (m) (x) ≥ 0, 0 ≤ x < ∞, m ∈ Z >0 . ( 4
)
According to Definition 2, a function g(x) that fulfills conditions ( 3) and ( 4) is a completely monotonic function. By (3), g(x) is strictly positive for x ≥ 0, including, in particular, the entire spectral interval 0 ≤ x ≤ λ n of the Laplacian L, and the function f (x) monotonically increases on its interval of definition 0 ≤ x < ∞. This monotonicity property implies that the algebraic multiplicity of the eigenvalues of L is maintained in

f (L) (in fact, f (λ i ) > f(λ j ) for λ i > λ j ). The function f (x)
is also given by the following integral

f (x) = x 0 g(y) dy, 0 ≤ x < ∞, (5)
i.e. it is the primitive of g(x) with f (0) = 0 and with f (x) > 0 for x > 0. In [22, Sect. 1.4.3], the authors have shown that relations (2)-( 5) are indeed sufficient to generate scalar C ∞ functions f (x) which define admissible Laplacian functions f (L) satisfying Conditions I-III. Note that f (x) in (5) has the general structure given in equation ( 2): hence, the admissible Laplacian functions can be represented as

f (L) = H(0)I n -H(-L) = n i = 2 (H(0) -H(-λ i )) v i v T i .
There is no term v 1 v T 1 since f (λ 1 ) = 0 (recall equation ( 1)). The eigenvector of f (L) associated with the zero eigenvalue is 1/ √ n, since L and f (L) have the same set of eigenvectors.

A. Examples of f (L)

There exist several classes of completely monotonic functions which in conjunction with equation ( 5) allow to construct admissible Laplacian functions satisfying Conditions I-III. We focus here on six examples (half of which are not reported in [START_REF] Michelitsch | Fractional Dynamics on Networks and Lattices[END_REF]), that are of special interest.

1) The function g(x) = (λx) m for 0 ≤ x < λ with m ∈ Z ≥0 and 0 < λ n < λ, satisfies conditions (3) and (4). From equation ( 5), we obtain

f (x) = x 0 (λ -y) m dy = 1 m + 1 λ m+1 -(λ -x) m+1 .
This yields the matrix function:

f (L) = 1 m + 1 λ m+1 I n -(λI n -L) m+1 = n i = 1 1 m + 1 λ m+1 -(λ -λ i ) m+1 v i v T i .
Note that in the trivial case of m = 0, the function g(x) = 1, and f (x) is the identity function, which yields f (L) = L. 2) Consider the function

g(x) = a (ax + b) m , a > 0, b ≥ 0, m ∈ Z >0 .
Let us start with the case of m = 1. From equation ( 5), we obtain f (x) = log(ax + b) for x ≥ 0. However, the constraint f (0) = 0, requires f (x) = log(ax+b)-log b. Hence, we have

f (L) = log(aL + bI n ) -log(bI n ).
Let us now consider the general case of m ∈ {2, 3, . . .}. Following the same procedure as above, we obtain:

f (x) = - 1 (m -1)(ax + b) m-1 , which we modify into f (x) = 1 m -1 - 1 (ax + b) m-1 + 1 b m-1 ,
to guarantee that f (0) = 0. Hence, we end up with

f (L) = 1 m -1 -(aL + bI n ) 1-m + b 1-m I n .
3) The exponential g(x) = ae -ax with a > 0, is a completely monotonic function. The corresponding function that retains the Laplacian structure is f (x) = 1e -ax and thus f (L) = I ne -aL . Matrix e -aL is doubly stochastic and it has appeared under different forms in consensus theory [START_REF] Mesbahi | Graph Theoretic Methods in Multiagent Networks[END_REF]. Functions like e -aL and the regularized Laplacian (I n + aL) -1 , have been also used as kernels to compute similarities between the nodes of an undirected graph [START_REF] Fouss | Algorithms and Models for Network Data and Link Analysis[END_REF]. 4) The function [START_REF] Miller | Completely monotonic functions[END_REF]. From (5), we obtain

g(x) = 1 b -ce -ax , a > 0, b > c > 0, defined on 0 < x < ∞, is completely monotonic
f (x) = 1 ab log(b -ce -ax ) + ax .
To guarantee f (0) = 0, we set

f (x) = 1 ab log(b -ce -ax ) + ax -log(b -c) , yielding f (L) = 1 ab log(bI n -ce -aL ) + aL -log(b -c)I n .
5) The function g(x) = (β + 1) x β with β ≤ 0 is a completely monotonic function. From equation ( 5), we obtain f (x) = x β+1 . The condition f (0) = 0, requires that -1 < β ≤ 0. Therefore, the function f (x) = x γ with γ ∈ R such that 0 < γ ≤ 1, preserves the Laplacian structure imposed by Conditions I-III.

We thus obtain the so-called fractional Laplacian, f (L) = L γ , which has been studied in the theory of random walks, and in diffusive and quantum transport on networks [START_REF] Michelitsch | Fractional Dynamics on Networks and Lattices[END_REF]Ch. 2]. Note that integer powers γ ∈ {2, 3, . . .}, fail to satisfy Condition III. 6) Node-invariant graph filters are linear graph-signal operators of the form [START_REF] Segarra | Optimal Graph-Filter Design and Applications to Distributed Linear Network Operators[END_REF], [START_REF] Segarra | Network Topology Inference from Spectral Templates[END_REF]:

f (L) = m-1 j = 0 b j L j , ( 6
)
where b 0 , b 1 , . . . , b m-1 are real coefficients, i.e. graph filters are polynomials of the Laplacian L (or of any other n × n matrix, such as the adjacency matrix A, whose sparsity pattern captures the local structure of graph G). A graph filter can be equivalently defined as

f (L) = a 0 m-1 j = 1 (L -a j I n ),
where a 0 , a 1 , . . . , a m-1 are real coefficients, which also gives rise to a polynomial on L of degree m -1.

From ( 6), we can see that Conditions I and II are fulfilled if b 0 = 0 and b j ≥ 0 for j ∈ {1, 2, . . . , m -1}, but Condition III is not satisfied, in general. Graph filters based on Chebyshev and Cayley polynomials have also been recently introduced to process graph-structured data [START_REF] Wu | A Comprehensive Survey on Graph Neural Networks[END_REF], [START_REF] Levie | CayleyNets: Graph Convolutional Neural Networks With Complex Rational Spectral Filters[END_REF].

B. General properties of f (L)

In the previous section, we have identified the family of functions that preserve the structure of the combinatorial Laplacian. We now briefly discuss some general properties of the matrix f (L), which will be used in Sect. IV.

1) Generalized degree: By construction, the diagonal entries of f (L) are positive, and analogously to the Laplacian L, we can refer to d i = f ii (L), i ∈ {1, 2, . . . , n}, as the generalized degree associated with the function f [22, Sect. 1.5.1]. The average of the generalized degree is defined as

1 n n i = 1 d i = 1 n trace(f (L)) = 1 n n i = 1 f (λ i ).
In the case of the fractional Laplacian, it is called the average fractional degree. Note that d i does not only capture the local information (on its nearest neighbors) of node i, but also incorporates knowledge at the level of the whole network. This "non-locality" property is further explored below.

2) Laplacian functions for regular graphs: To gain some insight into the structure of the Laplacian functions f (L), it is worth focusing on the special case of regular graphs. In a regular graph, each node has the same degree k, and the Laplacian takes the simple form L = kI n -A. By using the series expansion presented at the beginning of Sect. III, we can express f (L) as follows [22, Sect. 1.5.2]:

f (L) = ∞ = 1 c (kI n -A) = ∞ = 1 h = 0 c h k -h (-1) h A h .
(7) Equation ( 7) unveils a connection between f (L) and the integer powers of the adjacency matrix A. Recall that the (i, j) entry of A h for h ∈ Z >0 , is the number of all the possible paths connecting node i to node j with h edges, whereas the diagonal entry (i, i) of A h is the number of closed paths with h edges which start and end at the same node i [START_REF] Godsil | Algebraic Graph Theory[END_REF]. Therefore, equation [START_REF] Bullo | Lectures on Network Systems[END_REF] shows how the function f changes the local character of the Laplacian L and makes it a long-range operator: matrix f (L) is thus well suited to define dynamical processes with non-local interactions on networks.

3) Exponential decay of functions of banded matrices:

A band matrix is a sparse matrix whose nonzero entries are limited to a diagonal band, including the main diagonal and (zero or) more diagonals on either side. Given a matrix B ∈ C n×n , if all entries of B are zero outside a diagonally bordered band whose range is determined by r 1 , r 2 ∈ Z ≥0 (i.e. [B] ij = 0 if j < ir 1 or i + r 2 < j), then r 1 and r 2 are called the lower bandwidth and upper bandwidth of B, respectively. The bandwidth of B is the maximum of r 1 and r 2 , i.e. it is the integer r such that

[B] ij = 0 if |i -j| > r.
For example, a band matrix with r 1 = r 2 = 0 (r 1 = r 2 = 1) is a diagonal (tridiagonal) matrix.

For the entries of functions of Hermitian band matrices, we have the following exponentially decaying bound (further decay results applying to general matrices, are available in [START_REF] Benzi | Decay Bounds and O(n) Algorithms for Approximating Functions of Sparse Matrices[END_REF]):

Theorem 2 (Benzi & Golub [START_REF] Benzi | Bounds for the Entries of Matrix Functions with Applications to Preconditioning[END_REF]) Let B ∈ C n×n be Hermitian and of bandwidth r, and let the function f be analytic in an ellipse containing the spectrum of B. Then,

f (B) satisfies |f ij (B)| ≤ M |i-j| ,
where M is a constant and = μ 1/r , where μ ∈ (0, 1) depends only on f .

Theorem 2 shows that the entries of f (B) are bounded in an exponential decay manner away from the diagonal (the rate of decay depending on f ), with the bound decreasing as the bandwidth r decreases. Note that this does not necessarily mean that "decay to zero" is observed in practice [START_REF] Higham | Functions of Matrices: Theory and Computation[END_REF]Sect. 14.2].

If we turn our attention again to the Laplacian functions f (L), it is easy to verify that the following family of undirected graphs with n nodes, admits a Laplacian L of bandwidth r ∈ {1, 2, . . . , n -2}: B r n = (V, E) where {i, j} ∈ E if and only if j = i + s with i ∈ {1, 2, . . . , n -s} and s ∈ {1, 2, . . . , r}. Graph B 1 n coincides with P n (the path graph with n nodes), whose Laplacian is a tridiagonal matrix. We excluded the case of r = n -1, since it corresponds to the complete graph with n nodes, whose Laplacian has zero bandwidth (see the example in Fig. 1). Note that this family of graphs is not unique. In fact, the entries on the diagonal band of L(B r n ) are all different from zero, but the definition of band matrix does not preclude null values.

IV. APPLICATION OF LAPLACIAN FUNCTIONS TO CONSENSUS-BASED FORMATION CONTROL

In this section, we leverage the Laplacian functions f (L) studied in Sect. III, to design a new class of continuoustime consensus protocols. In particular, we will focus on the following generalization of the shape-based formation control strategy presented in [START_REF] Mesbahi | Graph Theoretic Methods in Multiagent Networks[END_REF]Ch. 6.3]: [START_REF] Morbidi | The deformed consensus protocol[END_REF] where the state vector x(t) = [x 1 (t), y 1 (t), . . . , x n (t), y n (t)] T ∈ R 2n contains the x-and y-coordinates of the positions of n single-integrator robots at time t ≥ 0. This strategy allows to drive the n robots to a rotationally-invariant formation encoded through the formation graph G F = (V, E F ) and the associated constant vector of target locations ξ ∈ R 2n . For the sake of simplicity, in the following, we will assume that G F coincides with the interaction (or communication) graph G = (V, E) of the robots (thus, the inclusion E F ⊆ E is trivially satisfied, cf. [6, Th. 6.12]).

ẋ(t) = (-f (L) ⊗ I 2 )(x(t) -ξ), x(0) = x 0 ,
Remark 1 (Invariance of the agreement value) By virtue of the spectral properties of the admissible Laplacian functions discussed in Sect. III, it is easy to prove that the state vector x(t) of system (8) asymptotically converges, for any initial condition x 0 ∈ R 2n , to a constant vector whose value does not depend on the f (L) chosen (in other words, the agreement value is f -invariant).

Remark 2 (Distributed implementation of (8)) While some admissible functions f conserve, at least in part, the sparsity pattern of the Laplacian L, the emergence of long-range interactions between the nodes of the graph typically translates into dense matrix functions of the Laplacian (i.e. the weighted interaction graph associated with f (L) tends to be fully connected). This means that to implement [START_REF] Morbidi | The deformed consensus protocol[END_REF], the robots should adopt an all-to-all communication pattern, which is undesirable, in practice. However, in many instances (cf. Theorem 2), a large percentage of the off-diagonal entries of f (L) is very close to zero, which correspond to network connections which bring negligible information to the nodes. This calls for an operator which most nearly transforms f (L) into a new matrix f q (L) which is sparser than f (L). A simple solution is to define f q (L) as follows:

f q ij (L) = 0 if |f ij (L)| < q, i = j, f ij (L) otherwise,
where q is a small positive constant (a threshold). This transformation does not break the symmetry of f (L) and the offdiagonal entries f q (L) remain non-positive as required by Condition III. Moreover, the Gershgorin circle theorem ensures that f q (L), as f (L), is positive semidefinite. In conclusion, while f q (L) still remains less sparse than L in general, it makes protocol (8) amenable to a distributed implementation (see Sect. IV-B for more details).

A. Case studies

For later use in Sect. IV-B, we recall here some elementary properties of four of the Laplacian functions examined in Sect. III-A.

1) Logarithmic function: Let f (x) = log(ax + 1) with a > 0. Note that

lim a → 0 + f (L) = log(aL + I n ) = 0 n . Moreover, f (L) grows unbounded, as a → ∞. 2) Exponential function: Let f (x) = 1 -e -ax with a > 0. Matrix f (L) = I n -e -aL ranges between 0 n (as a → 0 + ) and I n -1 n 11 T (as a → ∞). In fact lim a → ∞ e -aL = I n -LL # = 1 n 11 T ,
where L # denotes the group generalized inverse of L [START_REF] Bernstein | Matrix Mathematics: Theory, Facts, and Formulas[END_REF]Prop. 11.8.2]. Finally, if we consider the Loewner ordering (i.e. the partial ordering " " defined by the convex cone of positive semidefinite matrices), the following property holds true for f (L) = I n -e -aL : Eigenvalues of L(P 10 )

0 n f (L) I n - 1 n 11 T . 0 0.5 1 1.5 2 2.5 3 3.5
Eigenvalues of f (L(P 10 ))

f (x) = x f (x) = log(3x + 1) f (x) = 1 -e -3x f (x) = x 1/3 Fig. 2.
Mapping of the spectrum of L(P10) via four admissible functions f (x). The 10 eigenvalues of L(P10) are represented as black dots on the horizontal axis and those of f (L(P10)) as colored dots on the vertical axis.

3) Quadratic function: Let g(x) = λx for 0 ≤ x < λ, which is positive over the spectrum of L, if λ > λ n with 2-2 cos((n-1)π/n) ≤ λ n ≤ n. By choosing such a λ, from (5), we have that f

(x) = 1 2 [λ 2 -(λ -x) 2 ] = 1 2 x(2λ -x), which yields the quadratic matrix function f (L) = -1 2 L 2 + λL. 4) Fractional power: Let f (x) = x γ with 0 < γ ≤ 1.
Matrix f (L) = L γ ranges between I n (as γ → 0 + ) and L (for γ = 1).

For the sake of illustration, Fig. 2 shows how the functions

f (x) = x (red), f (x) = log(3x + 1) (green), f (x) = 1 -e -3x
(blue), and f (x) = x 1/3 (magenta), map the 10 eigenvalues 2 -2 cos(kπ/10), k ∈ {0, 1, . . . , 9}, of the Laplacian L of the path graph P 10 . Note that the exponential function maps the spectrum of L(P 10 ) into the [0, 1) interval: in fact 1-e -3λn = 1 -8 × 10 -6 .

B. Numerical results

In our numerical experiments, we simulated system (8) for 10 robots, using different Laplacian functions, interaction graphs and vectors of target locations. We ran 20-second simulations using Matlab ode45 solver with a variable step (max step size: 0.01 s).

We used Matlab built-in commands logm and expm for calculating the matrix logarithm and matrix exponential, respectively. The computation of the fractional Laplacian being more delicate, we opted for the rootpm real routine of The Matrix Function Toolbox [START_REF] Higham | The Matrix Function Toolbox[END_REF] (see also [START_REF] Higham | Functions of Matrices: Theory and Computation[END_REF]Appendix D]). This routine computes the pth root of a real matrix via the real Schur form: similar results were obtained with rootpm sign which computes the pth root via the matrix sign function. Both routines handle singular matrices (L has a zero eigenvalue), and return non-principal pth roots.

In our first example, the target formation is a regular decagon of unit radius centered at the origin, i.e. ξ = [cos(0), sin(0), cos(2π/10), sin(2π/10), . . . , cos(9π/10), sin(9π/10)] T ∈ R 20 . The interaction graph is G = P 10 , and the Laplacian functions considered are f (L) = L, f (L) = log(3L + I 10 ) and f (L) = I 10e -3L . The initial positions x 0 of the 10 robots have been randomly generated by drawing the x-and y-coordinates from the standard uniform distribution on the open interval (0, 1). The first row of Fig. 3 reports the trajectory of the robots for each of the three functions (the edges of the formation graph are solid black, and the initial positions of the robots are marked with stars). The second row of Fig. 3 shows the corresponding time-evolution of the formation error e e e(t) = x(t) -ξ (x-coordinates, top; y-coordinates, bottom). In the three cases, e e e(t) → [τ x , τ y , . . . , τ x , τ y ] T as t → ∞, where τ x 0.3451 and τ y 0.5841 are constant offsets. The average of the generalized degree associated with the three Laplacian functions is 1.8, 1.54 and 0.7833, respectively (recall Sect. III-B.1), and their second smallest eigenvalues are λ 2 = 2 -2 cos(π/10) 0.0979, log(3λ 2 + 1) 0.2575 and 1e -3λ2 0.2545, respectively. Therefore, the logarithmic and exponential functions accelerate convergence towards the desired formation, and the convergence speed can be adjusted by tuning a single real parameter (the positive scalar a). However, both functions keep the agreement value unchanged.

To study the impact of the approximation introduced in Remark 2, in Fig. 4 we considered the same interaction graph, target formation, and initial conditions as in Figs. 3(b),(e), but we replaced f (L) = log(3L + I 10 ) with f q (L), where the threshold q = 0.01. Note that the exponential decay property of Theorem 2 holds for f (L) = log(3L + I 10 ) with r = 1. Figs. 4(a),(b) report the trajectory of the 10 robots and the time-evolution of the formation error e e e(t), respectively, and they qualitatively show that the proposed approximation has a negligible effect on the dynamic behavior of system [START_REF] Morbidi | The deformed consensus protocol[END_REF]. More specifically, at the end-time, t = 20 s, we have that e e e(t)

2.1453 for f (L) and e e e(t) 1.7851 for f q (L), where • denotes the Euclidean norm. However, the weighted interaction graph associated with f q (L) is sparser than that associated with f (L). This is evident in Fig. 4(c), which shows the sparsity pattern of the adjacency matrices diag(f 1,1 (L), . . . , f 10,10 (L))f (L) (left) and diag(f q 1,1 (L), . . . , f q 10,10 (L))f q (L) (right), where the zero entries are white and the nonzero entries are black.

0.5 1 _ _ x [m] y [m] R1 R2 R3 R4 R5 R6 R7 R8 R9 R10 (b) f (L) = -1 2 L 2 + 5L 1
In our second example, we chose the interaction graph G = C 10 , the cycle graph with 10 nodes, and the desired formation is the pentagram (or five-pointed star). Let where j ∈ {1, 3, 5, . . . , 19}. The Laplacian functions considered in this second example, are f (L) = L, f (L) = -1 2 L 2 + 5L and f (L) = L 1/5 . In the quadratic function, we set λ = 5 > λ 10 = 4, to satisfy the condition discussed in Sect. IV-A. The vector of initial positions x 0 , is not random this time: the 10 robots are initially placed at the vertices of the same regular decagon of unit radius considered as target formation in the first example. As in Fig. 3, the first row of Fig. 5 reports the trajectory of the 10 robots for the three Laplacian functions, and the second row, the corresponding time-evolution of the formation error e e e(t) = x(t) -ξ. Since the regular decagon is centered at the origin, we have that e e e(t) → 0 as t → ∞, in all cases. Finally, note that the average fractional degree of L 1/5 is 1.0115, and that the user can take advantage of the non-integer power γ to modulate the convergence speed of the 10 robots towards the desired formation. In fact, the second smallest eigenvalues of L and L 1/5 , are λ 2 = 2 -2 cos(π/5) 0.3820 and λ 1/5 2 0.8249, respectively, which explains the faster convergence rate observed in Figs. 5(c),(f).

V. CONCLUSION AND FUTURE WORK

In this paper, we have explored a general class of matrix functions of the combinatorial Laplacian, that retain its structural properties. This allowed us to present under a common framework, several variants of the Laplacian scattered across different domains. For the first time, this family of Laplacian functions has been utilized in a consensus-based formation control protocol, revealing some attractive features. In fact, the selected functions do not alter the agreement value, and offer greater design flexibility. For example, our numerical experiments show that in many situations, a single scalar parameter is sufficient to adjust the convergence speed towards consensus without drastically increasing the number of communication exchanges between the agents.

This work provides fertile ground for further research and experimentation. One line of future research will be to relax Condition III and consider a broader class of admissible Laplacian functions. As far as the consensus protocol is concerned, we see potential for extension to weighted directed interaction graphs, and we plan to perform hardware experiments with mobile robots. Finally, in [START_REF] Michelitsch | Fractional Dynamics on Networks and Lattices[END_REF], the authors have shown a connection between the fractional Laplacian of a graph and the operators in fractional calculus. There might then exist a link between the agreement protocol driven by the fractional Laplacian, and the coordination algorithms for fractional-order systems studied in [START_REF] Cao | Distributed Coordination of Networked Fractional-Order Systems[END_REF].

4 Fig. 1 .

 41 Fig. 1. Graphs B r 5 for r ∈ {1, 2, 3}. Note that B 4 5 (light gray) is the complete graph with 5 nodes.

Fig. 3 .Fig. 4 .

 34 Fig. 3. First example, ξ is the regular decagon and G = P10: (1st row) Trajectory of the 10 robots. The initial and final positions are marked with a star " * " and a hollow circle "•", respectively, and the edges of the formation graph are solid black; (2nd row) Time evolution of the formation error e e e(t) = x(t) -ξ of the 10 robots (for color coding, see the legends in the 1st row). (a),(d) f (L) = L, (b),(e) f (L) = log(3L + I10), and (c),(f) f (L) = I10 -e -3L .

5 Fig. 5 . 2 L 2 +

 5522 Fig. 5. Second example, ξ is the pentagram and G = C10: (1st row) Trajectory of the 10 robots. The initial and final positions are marked with a star " * " and a hollow circle "•", respectively, and the edges of the formation graph are solid black; (2nd row) Time evolution of the formation error e e e(t) = x(t) -ξ of the 10 robots (for color coding, see the legends in the 1st row). (a),(d) f (L) = L, (b),(e) f (L) = -1 2 L 2 + 5L, and (c),(f) f (L) = L 1/5 .

  be the circumradius of the pentagram and the circumradius of its inner pentagon, respectively. Then, the target locations of the 10 robots can be expressed in polar coordinates as jπ 10 ), sin( jπ 10 ) T if j ∈ {1, 5, 9, . . .}, R cos( jπ 10 ), sin( jπ 10 ) T if j ∈ {3, 7, 11, . . .},

For more details on primary and nonprimary matrix functions, the reader is referred to[START_REF] Higham | Functions of Matrices: Theory and Computation[END_REF] Sect. 1.4].


Thi B H Tran

Patrice Huguet

Arnaud Morin

Mike Robitzer

Stefano Deabate
email: stefano.deabate@umontpellier.fr

Operando µ-Raman Measurement of the Water Distribution Along and Across the Membrane in the Fuel Cell

To obtain a fundamental understanding of the transport regimes governing the electrochemical behavior and ageing of the polymer electrolyte membrane fuel cell, the measurement of the membrane in-and through-plane water distribution is essential. In this paper, operando -Raman spectroscopy is used to probe the water content of the Nafion ® membrane in a model fuel cell working at constant stoichiometry and relative humidity. Water concentration profiles crossing the membrane thickness are obtained with m resolution at different locations of the active surface: at the middle and close to the inlets/outlets of the reactant gases, at the feed gas channel and under-lands areas. The influence of the operating temperature and of the current density delivered by the cell are investigated. Of particular interest for perspective fuel cells, the membrane inner water partition appears highly heterogeneous when operating at ambient conditions, and dependent on the delivered current. The increase of temperature decreases the membrane water content at all the probed positions and induces a less uneven water distribution. The membrane hydration appears to be a key parameter for understanding the water redistribution between cathode and anode with the change of the cell operating conditions.

Introduction

Fuel cells (FC) represent a highly efficiency energy conversion technology. Combined with dihydrogen fuel, a clean energy carrier, the polymer electrolyte membrane fuel cell (PEMFC) technology has a great potential to reduce both carbon dioxide emissions and dependence on fossil fuels. Thus, PEMFCs are currently regarded as strategic products by major industrial countries. 1,2 Nevertheless, two main interrelated technical challenges delay PEMFC large-scale marketing: cost and durability. [3][4][START_REF]Fuel Cell System Cost -2017[END_REF] Consequently, improving water and heat management arises as a major factor. [6][7][8] The water distribution within the electrochemical system results from the interplay between in-and through-plane fluxes. The in-plane water transfer is mainly governed by the flow of the reactant gases and liquid water accumulation along the distribution channel. Simultaneously, water is redistributed by the membrane between the opposite electrodes, via coupled mass and charge transport processes involving both H2O and H + species: mainly backdiffusion and electroosmosis. The significance of other contributions, e.g. convective or thermoosmotic, remains under debate. 6 The occurring of coupled transport phenomena and mechanisms, together with the complexity of the FC design (consisting in a layered structure of media with different porosity and permeability), results in heterogeneous through-and in-plane conditions concurrently affecting temperature, current, concentration and pressure of reactants and water. Accordingly, the way the FC performs and degrades is also inhomogeneous, difficult to predict and manage.

The development of suitable models of the PEMFC electrochemical behavior, performances and ageing requires the fundamental understanding of the different transport mechanisms and, therefore, the detailed knowledge of the through-plane water distribution across the different components of the membrane-electrode assembly (MEA). The polymer membrane, placed at the middle of the electrochemical system and separating the electrodes, appears as the keystone of the 3D water distribution and represents an example of the typical back-coupling of transport phenomena which confers particular complexity to the water distribution issue. Indeed, the membrane permeability to water depends on the membrane water content and the last depends, in turn, on the external water concentration as it is determined by the in-and through-plane fluxes. Further, the water sorption ability of the membrane presumably depends on the different, local environmental conditions inside the FC: typically, the uneven constriction occurring at the scale of the alternating gas distribution channels and lands for the current collection, hindering the homogeneous swelling of the membrane and locally altering the porosity of the active (ALs) and gas diffusion layers (GDLs). The information of the local water content in the membrane during operation is, then, a critical issue for the validation of mass transfer models, the definition of optimized operating conditions and the development of efficient systems with innovative designs for efficient water management. A consequent research effort has been carried out on the development of operando techniques able to probe the water content and partition in the membrane. 7,9,10 Among these, operando Raman microspectroscopy has been developed more recently by a limited number of research groups. [11][12][13][14][15][16][17][18][19][20] Operando -Raman is a laboratory-scale, sensitive and selective technique based on the accumulation of Raman spectra across the membrane thickness. Once the raw signal is properly corrected for the different optical and spectroscopic aberrations, the concentration profile of the water crossing the membrane is obtained with both in-and through-plane m spatial resolution. More details about the principle of operando -Raman and depth-profiling measurements through perfluorosulfonic polymers can be found in Deabate et al.. 17 In our previous article, 20 we reported the evolution of the membrane hydration at the middle of the active surface as a function of the FC operating temperature and delivered current. The work here reported completes and goes further. Water concentration profiles are collected close to the inlets and outlets of the reactant gases, at locations corresponding to the gas feed channel and to the lands for the current collection, at operating temperatures ranging from ambient conditions to 80°C. In order to progress in the understanding of limiting phenomena occurring under the dry operating conditions currently envisaged and relationships with water transfers across the MEA, the cell is operated at low relative humidity (20% RH) and high flows of the reactant gases. Results allow to map the uneven water partition both in-and troughplane and attest the crosswise relationships between the water distribution, the current delivered by the cell and the FC operating temperature.

Experimental

Cell. The model cell used for operando -Raman measurements is the same used in Tran et al. 20 and consists of stainless steel monopolar plates with a single serpentine gas distribution channel (active area 4 cm 2 , channel width and depth 1 and 0.5 mm respectively, lands width 1.2 mm) coated by a gold layer (10 m thick) improving electrical contact and preventing steel corrosion. To provide the access of the Raman exciting irradiation to the membrane, the upper plate (cathode side) is drilled at the center (1 mm hole). A round cover glass and a compressible gasket squeezed by a tightening ring ensure the tightness. Four distinct cathode BPs are used, differing for the position of the hole which is placed either at the gas distribution channel or at the land for the current collection, at the middle of the active surface or close to the inlets/outlets of the gaseous reactives (Figure 1). Materials and Assembly. Nafion ® N115 membrane (thickness ~ 125 m at 50% RH, EW = 1100 g mol -1 ) from DuPont De Nemours is used as received. Electrodes are Johnson Matthey on Toray paper TGP-H-060 (thickness 190 m, 0.4 mg cm -2 Pt). To provide the access of the Raman exciting irradiation to the membrane, GDL and AL are removed over a circular surface (∅400 m) as previously reported. 20 The cell is finally clamped with a tightening torque of 5 Nm and placed at 120°C during 1 hour. FC Operating Conditions. The PEM fuel cell is operated with pure oxygen and hydrogen at atmospheric pressure. The input and output conditions (gases temperature, pressure, flow rate and humidity) are managed and monitored by a home-made test bench. RH is measured by a Rotronic HygroPalm capacitive sensor. The cell temperature is managed by a Pt 100 probe placed at the cathode plate. The reactants are supplied at the same temperature as the cell, in counter-flow arrangement (anode and cathode channels face each other across the membrane and respective flows have opposite directions), at constant RH (20% at both inlets) and stoichiometry (25/100 for H2/O2 respectively). The high flow rates are needed to avoid water condensation in the channel and at the optical window (the last representing a cold point).

The electrochemical management of the cell is carried out by a Biologic SP-150 potentiostat equipped with a VMP3 booster. The FC polarization curves are established point by point using constant current chronopotentiometry. The impedance of the cell is measured at each point by electrochemical impedance spectroscopy (EIS) (current perturbation 10 mA, frequency range 2 × 10 -1 to 5 × 105 Hz).

Previous to Raman measurements, the PEMFC is operated during at least 30min, until to attain steady conditions of voltage, ohmic resistance and outlet gas RH.

Mass Balance Measurements. Parameters regarding the water mass balance measurement, i.e. gas flow rates supplied at the electrodes, RH, pressure at the inlet and outlet and cell temperature, are recorded continuously during experiments. The amount of water (D / mol s -1 ) incoming and outgoing each electrode is calculated for each operating condition. The extra water flow (J / mol s -1 cm -2 ), which corresponds to the molar flow rate of water from the electrochemical reaction, is:

𝐽 = (𝐷 -𝐷 )/𝑆 (1)
where 𝐷 and 𝐷 are the molar flow rates of water at the outlet and inlet, respectively, and S is the cell active area (S / cm 2 ). The ratios between the extra water measured at each electrode by the RH probes and the overall flow related to the whole water produced by the electrochemical reaction, i.e. calculated from the current delivered by the FC, are:

𝑅 = 𝐽 𝐽 ⁄ (2) 𝑅 = 𝐽 𝐽 ⁄ (3)
It follows that:

𝑅 = 𝑅 + 𝑅 = (𝐽 + 𝐽 ) 𝐽 ⁄ (4)
Since 𝐽 and 𝐽 correspond to experimental values from the RH measurements, while 𝐽 is calculated from the current delivered by the cell, R tot is expected to be equal to 1 when the FC operates at the steady state.

Raman Measurements. The PEMFC is placed on an automatic sliding plate and depth-resolved measurements are carried out across the membrane thickness with a step of 1 m (Figure S1). The experimental in-depth resolution is 7±1 μm. 17 The experimental volume probed at a given position is ≤ 14 µm 3 . The collection time for a single spectrum is limited to 1 s to avoid local warming and ~ 220 spectra are cumulated at different positions across the membrane thickness for a single cartography. To ensure high statistics, 15 cartographies are carried out for each experiment and their average performed.

Raman spectra are excited with a 647.1 nm irradiation generated by an Ar-Kr laser. The laser power at the membrane surface is 20 mW. The back-scattered radiation is recovered by a 50X Fluotar objective (maintained at 30°C to avoid damage) and transmitted to a LABRAM1B confocal Raman spectrometer (Jobin-Yvon S.A., Horiba, France) equipped with a charge coupled device detector with a 300 grooves mm-1 grating and cooled by a double Peltier effect.

The water volume fraction V and hydration number  of the membrane are calculated at each probed position from the ratio between the intensities of the (O-H) mode and some selected Raman bands belonging to the ionomer chemical groups, as detailed elsewhere. 17,20,21 The soobtained sorption isotherms of the Nafion ® 115 membrane, attesting the reliability and sensitivity of the method, are reported in Figure S2.

The water concentration-profiles obtained by depth-resolved Raman measurements are affected by artificial gradients and poor resolution. Actual water profiles are then recovered by the geometrical correction of spherical aberrations and by a convolution procedure using the instrumental spreading function, according to the procedure detailed in Deabate et al.. 17 After treatment of the raw data, the depth-resolution of the -Raman measurement is assumed to correspond to ~ 1 µm (probed volume  1 m 3 ). 20 and attest the reproducibility of measurements obtained with the different BPs. Since operating conditions, cell design and materials are not optimized for efficiency but for operando microspectroscopy and for enhancing limiting phenomena related to low RH operation, the electrochemical performances of the model FC are lower than the state-of-the-art. It should therefore be pointed out that performances do not affect neither the nature of the water transport phenomena involved nor their relationships with the FC electrochemical behavior, so that results obtained here can be used for the better understanding of state-of-the-art cells.

Results and Discussion

The FC performance increases with T (Figure 2a), which corresponds to the general decrease of the integral ohmic resistance (Figure 2b), the last attributable (as a first approximation) to the simultaneous enhancing of proton mobility in the ionomer phase (of both the electrodes and the membrane), reactive gases diffusion and catalyst activity. Concurrently, ohmic losses are also affected by the current density j delivered by the FC. At T  50°C, the cell resistance appears to increases with the increase of j. At 80°C, resistance initially drops and then stay unaffected. This last behavior has already been reported under similar operating conditions and attributed either to lack of temperature management (i.e. resistance values are recorded before the temperature stabilization) 22 or to lacking membrane hydration when using low hydrated dioxygen at low current densities. 23 Neither hypotheses can be applied here since steady conditions have been carefully checked and -Raman measurements attest no significant changes of the membrane water content with j at this temperature (vide infra). As discussed in the following, the particular behavior of the ohmic resistance observed at 80°C (and, to a lesser extent at 50°C) can be attributed to the evolution with j of the water amount transferred from the cathode across the membrane and hydrating the anode ionomer phase.

Figures 3 and4 show the evolution with current of the water concentration gradients crossing the membrane in the PEMFC operating at Troom and 80°C, respectively. The behavior observed at the intermediate T is in-between these two extreme cases (see Figures S3 andS4). All water concentration profiles are practically linear, attesting that the value of the water diffusion coefficient does not vary significantly along the membrane thickness for a given T and j. Confirming our previous results, 18,20 the membrane average water content decreases with the increase of T (Figure 5). When the cell is operated at T  50°C, the membrane hydration also depends on the applied current density: it first increases at the FC start up and then decreases with j. The last results allow to explain the general increase of the FC resistance with j (Figure 2b). On the whole, the most hydrated position corresponds to the middle of the active area and the lower the cathode inlet/anode outlet (see also Figure S5). The highest water content is observed at the anode channel (V = 0.133 at 25 mA cm -2 and Troom) but, overall, the membrane appears more hydrated under the lands for the current collection than at the channel. It should be noted that the hydration difference at the lands/channel scale was yet observed by -Raman in a FC with parallel channels but, in that case, the phenomenon was largely more enhanced. 20 The discrepancy between the water contents at the under-lands and channel positions increases with the increase of j, which strongly suggests the accumulation of water from the electrochemical reaction at the lands because of simultaneous less effective drainage by gases and more effective heat evacuation by lands. A further, concomitant reason could be the fact that the channel has been identified as the area hosting the larger local currents, 24,25 therefore subjected to higher local warming inducing lower actual RH.

The membrane water loss with T was previously attributed to the relative amount of water produced by the electrochemically reaction becoming negligible compared to that injected via the reactants hydration when the cell T increases: the water activity of the gases at the inlets staying invariant, the actual RH inside the cell decreases. 20 Indeed, the maximum values measured at Troom for the RH at the cathode and anode outlets are 65% and 45% respectively, while both stay very close to 20 % at 80°C. The membrane dehydration with j was also previously reported 18,20 and mainly ascribed to the spontaneous increase of the FC inner T due to the increase of electrochemical heat and ohmic losses, with a lesser contribution from to the increase of pressure drops due to the increase of the gas reactants flows with current (inducing the decreasing of the actual amount of water injected in the cell). The observation that the membrane water content becomes unaffected by j at T = 80°C agrees with the consistent decrease of the ohmic losses due to the enhanced proton mobility at this temperature (which implies the lower spontaneous warming of the cell for a given increase of j) and to the concomitant lower impact of any T increase on the actual RH (because of the higher water saturation pressure). Results here reported agree with EIS measurements and definitively confirm what was previously hypothesized, 18,20 namely that the membrane dehydration phenomena observed with the increase of T and j do not concern only the middle of the active surface but equally involve the whole membrane volume.

Membrane inner water gradients (V/d where V = V / cathode -V / anode and d is the membrane thickness as measured by -Raman, see Figures S6 andS7) are observed whatever the probed position and the operating conditions, even if strongly affected by all the investigated parameters. Data show that: (i) at T  50°C, the inlet and outlet positions (referring to the cathode) show opposite gradients (particularly at the lower j), which should correspond to the occurrence of diffusion water fluxes in opposite direction at these locations as expected for the counter-flow configuration; (ii) gradients measured at the middle of the active surface and, to a lesser extent, at the under-lands outlet position exhibit intermediate behavior i.e. initially directed from the anode to the cathode, the gradient progressively reverses with the increase of j (which strongly suggests the inversion with the current increase of the external water gradient and, consequently, water fluxes crossing the membrane); (iii) on the whole, gradients are largest at the inlets/outlets positions than at the center (but this difference becomes negligible with the increase of both T and j); (iv) on the whole, gradients are slightly less marked at the under-lands than at the channel area (for a given location at the active surface); (v) but for the under-lands position at the cathode outlet at Troom and T = 32°C, gradients decrease with the increase of both T and j. At T = 80°C, the behavior of the weak gradients remaining becomes completely erratic.

As far as self-or low-humidified PEMFC are concerned, experimental and theoretical studies have so far considered the water concentration gradient at the membrane always oriented from the cathode to the anode, independently from the direction of the net water flux (see Feindel et al. 26 and references therein). Then, a main result here is to show that inner gradients (and so, diffusion fluxes which most often represent the main contribution) can exhibit opposite direction depending on the position at the active surface and can reverse as the current delivered by the FC varies.

The decrease of the membrane inner water gradient with T reflects the decrease of the external gradient, i.e. the difference between the gaseous water activities at the anode and cathode sides of the cell. This is due to both (i) the fact that the amount of the water produced by the electrochemical reaction at the cathode becomes progressively less weighty compared to that injected into the cell by hydration of the reactants and (ii) the increase of the MEA (particularly the membrane) permeability with T, 27 enhancing water flows due to backdiffusion. Water mass balance measurements reported in Figure 6 clearly show that, even if the gaseous water activity measured at the FC cathode outlet stays always larger (as expected) for all the temperature and current conditions, the transfer of water from the electrochemical reaction to the anode increases with T.

The decrease of the membrane inner water gradient with j (observed for the most of the probed positions) is related to the larger dehydration of the membrane anode side (initially more hydrated) compared to the cathode for a given current increase (see Figure 3, S3 and S4). In principle, assuming that the spontaneous warming of the cell with j should equally affect both sides of the membrane (at least at the macroscopic scale), increasing j at a given T should correspond to the increase of the external water gradient. In spite of this, the concomitant decrease of the membrane permeability due to dehydration (as attested by -Raman measurements), hindering water transfer to the cathode, and enhancing of the eletroosmosis flows with j (also reinforced by the increase of the electroosmotic drag coefficient with the membrane water loss 28 ) induce the membrane anodic side to dehydrate more rapidly. The different evolutions with j of the water transfer across the MEA observed at different T (Figure 6) result then from the interplay between the antagonistic effects of temperature and current on the net through-plane water flux and their relationships with the membrane water content. At the lower T (Troom and 32°C), when the membrane hydration changes considerably with j (i.e. decreases), water transfer from the cathode to the anode decreases continuously with the current increase, which attest the global loss of the membrane permeability, overwhelming any effect due to the spontaneous warming of the cell. Differently, at the operating T = 50°C and more particularly 80°C, when the membrane hydration stays almost unaffected by current, water transfer from the cathode to the anode first increases with j, presumably because of the increase of the external gradient, and then decreases with the prevailing of the electroosmosis flow over back-diffusion. Since, at the maximum current, R cathode = 1 at the lower T while remaining < 0.8 at the higher T, these results suggest that a better water removal from the cathode by transfer through the membrane to the anode should be included among the different phenomena improving the FC performance with the increase of the operating T. UiR-corrected reported in Figure 2c confirms that, at 80°C, losses at low j mainly originates from the electrodes, presumably because of lower actual local RH.

The m resolution of the depth-resolved measurements allows access to the water content of the membrane surface i.e. to compare the membrane hydration behavior at the opposite sides. Figure 7 represents the water distribution at the anode and cathode membrane surfaces as a function of the FC operating T and j (according to data reported in Figure S5, the evolution of the water content between two next probed positions being assumed as continuous). Both surfaces exhibit unevenly water content at the channel-land scale as well as along the inlet to outlet direction but, as long as T  50°C, the heterogeneity of the water distributions at the opposite sides is completely different. This attest the complexity of the interplay between the different water transport mechanisms (diffusion, electroosmosis and convection), FC design (serpentine configuration of the channels) and operation conditions (counter-flow feeding). At Troom, the membrane cathode side, directly exposed to water from the electrochemical reaction, exhibits gradual water accumulation from the inlet to the outlet, with the driest area corresponding to the channel inlet (where water evacuation by gas is more effective) while identical water concentrations are observed at the outlet channel and under-land positions. Somehow unexpected, the cathode side appears globally less hydrated than the anode. Differently from the cathode, but as usually expected for the FC fed counter-flow, the anode side shows higher water accumulation at the middle of the active surface, particularly at the channel. The inlet and outlet areas exhibit very similar water contents, with the channel as the less hydrated position. The increase of both T and j dehydrates both sides of the membrane to very similar water contents and decreases all differences between channel and ribs, inlets and outlets until an almost homogeneous water distribution is attained at 80°C, close to the membrane hydration state observed at the equilibrium at 20% RH.

Of interest for PEMFC intended to operate at ambient conditions, the coupling of macroand microscopic measurements allows to characterize the evolution of the actual feed gases RH along the flow path as a function of the cell operating conditions. Figure 8 reports the example of the RH evolution at the FC cathode and anode sides as a function of the delivered j, at T = 32°C. Values at the ends of the channels (full symbols) correspond to the direct measurements of the water activity by the RH probe while the others (open symbols) are deduced from the water content of the membrane surface according to -Raman i.e. water contents are converted into RH values using the sorption isotherm reported in Figure S2 and assuming the membrane surface at the equilibrium with the gas phase. At the cathode side, the water activity always increases from the inlet to the outlet, almost linearly at low j but with a larger water accumulation close to the end of the flow path when the current increase (a larger spread between macro-and microscopic measurements is observed at this area). The water vapor concentration distribution along the anode is more strongly affected by j. At low current, the water activity first increases from the inlet to exceed values measured at the cathode and reach the maximum at the midstream area, presumably owing to the strong back-diffusion from the highly wetted cathode at this position. Then, moving to the downstream region, water concentration decreases due to the forward diffusion to the drier cathode and larger electroosmosis (the last increasing because of the lower average water content of the membrane at this area 28 ) but always remaining larger than at the cathode. Thus, at the lower j, the anode side appears to have larger RH along most part of the flow path. The opposite is observed at larger currents. With the increasing of j, the position where anode and cathode show the same water activity moves progressively from the cathode outlet/anode inlet region towards the opposite side. The global value of the water activity decreases at both sides but more rapidly at the anode channel. At the largest j, the RH stays homogeneous and very close to 20% along the anode flow path (despite the larger water transfer from the cathode, Figure 6) while increases suddenly at the cathode downstream area.

Conclusions

Raman microspectroscopy has been applied to improve the knowledge of the steady-state water distributions occurring in the PEMFC under a range of operating temperatures. This work, even if mainly representing a phenomenological description of the Nafion hydration behavior under specific operating conditions, highlights aspects of the water content distribution which could hardly have been anticipated by modeling. This points out the concern of operando investigation to address the actual behavior of perfluorosulfonic materials and produce useful inputs for forthcoming numerical studies.

As far as low T are concerned (which, together with low humidity represents the operating conditions targeted for automotive applications), results attest the highly heterogeneous in-and through-plane water distribution inside the membrane, and this despite the high flow rates used in this work. The membrane hydration, strongly dependent on both the operating temperature and the current delivered by the FC, in turn affects (together with T) the flows redistributing gaseous water between the cathode and anode compartments of the cell. This implies that the effective modeling of through-plane water fluxes requires the accurate knowledge of the local water content of the membrane.

At low T, two different water distribution are observed at the opposite surfaces of the membrane. Hydration at the cathode side is governed by the accumulation of the electrochemically produced water along the flow path. Water distribution at the anode results from the counter-flow management of gaseous reactants and exhibits larger water concentration at the midstream area. Under the specific operating conditions concerning this work, i.e. constant stoichiometry and RH, the increase of both FC temperature and current dehydrates both sides of the membrane to almost identical and homogeneous water contents. This implies that the effect of current on the water distribution becomes negligible when the cell is operated at the standard T = 80°C.

Globally, the FC performances appear to strongly depend on two main properties of the membrane: ohmic resistance (strictly related to the water content) and water permeability, the last governing the adequate removal of water from the cathode and, thus, probably ensuring the better hydration of the electrodes ionomer phase. The evolution of both resistance and permeability with j affects concomitantly the cell electrochemical behavior up to the operating T = 50°C. At T = 80°C, the FC performances are no longer dependent from the membrane hydration, the last staying unaffected by current, and the membrane ability to transfer water from the cathode to the anode becomes the main parameter affecting performances at low j. Forthcoming operando -Raman studies will be devoted, on the one hand, to extend our investigation to further operating conditions affecting mass transport losses (e.g. cell compression, land/channel size and design, weak stoichiometry…) and, on the other, to improve experimental conditions and cell design as to investigate state-of-art MEA.

Figure 2

 2 Figure 2 reports representative polarization curves of the model PEMFC operated at different T (Troom corresponding to around 25°C) and related ohmic resistances. Results are very similar to those previously reported20 and attest the reproducibility of measurements obtained with the different BPs. Since operating conditions, cell design and materials are not optimized for efficiency but for operando microspectroscopy and for enhancing limiting phenomena related to low RH operation, the electrochemical performances of the model FC are lower than the state-of-the-art. It should therefore be pointed out that performances do not affect neither the nature of the water transport phenomena involved nor their relationships with the FC electrochemical behavior, so that results obtained here can be used for the better

Figures

  Figures

Figure 1 .

 1 Figure 1. In scale schematic representation of the configuration of the gas distribution channel and positions probed by the Raman exciting irradiation.

  Figure 1. In scale schematic representation of the configuration of the gas distribution channel and positions probed by the Raman exciting irradiation.

Figure 2 .

 2 Figure 2. (a) FC polarization curves at different operating T (reported in the legend) and corresponding evolutions of (b) the ohmic resistance and (c) the iR corrected potential as a function of the delivered current.

Figure 3 .

 3 Figure 3. Comparison of the evolution of the membrane inner water distribution with the current density (reported in the legend) at the different probed position (top: channel, bottom: lands, left: cathode inlet, middle: center, right: cathode outlet) in the cell operated at Troom.

Figure 4 .

 4 Figure 4. Comparison of the evolution of the membrane inner water distribution with the current density (reported in the legend) at the different probed position (top: channel, bottom: lands, left: cathode inlet, middle: center, right: cathode outlet) in the cell operated at 80°C.

Figure 5 .

 5 Figure 5. Evolution with the current density of the membrane average water volume fraction at the channel (open symbols) and under-lands (full symbols) positions in the cell operated at different T. The minimum and maximum water contents measured for a given T are also reported as hydration number  = [H2O]/[SO3].

Figure 6 .

 6 Figure 6. (a) Evolution with current of the extra water flows measured at the cathode (open symbols) and anode (full symbols) sides (the evolution of the overall flow of the electrochemically produced water is also reported in full black symbol; errors are inside the symbol size) and of (b) the ratios between the extra water flow measured at each electrode and the overall extra water, at the different FC operating T (reported in the legend).

Figure 7 .

 7 Figure 7. Schematic representation of the water distribution at the anode and cathode membrane surfaces and the effects of the FC operating T and delivered j.

Figure 8 .

 8 Figure 8. Evolution of the RH distribution along the anode (red) and cathode (blue) gas flow paths with the current density (reported in the legend) in the FC operated at T = 32°C. Full symbols correspond to direct measurements at the FC inlets and outlets whereas open symbols are RH values deduced from the membrane surface water content as given by -Raman.


Isabelle Trenque

Stéphane Mornet

Etienne Duguet

Manuel Gaudon

New Insights into Crystallite Size and Cell Parameters Correlation for ZnO Nanoparticles Obtained from Polyol-Mediated Synthesis

ZnO nanocrystals were prepared from polyol-mediated synthesis. Two key parameters, that is, the zinc precursor concentration and the alcohol mixture chosen as synthesis medium, were varied. The increase of the precursor concentration and the decrease of the permittivity of the alcohol mixture were shown to favor the crystallite growth, leading to crystallite sizes ranging from 5 to 35 nm. The aggregation behavior of the nanocrystal units to form or not polycrystalline spheres, depending on the van der Walls interaction density, was shown and explained.

Every sample was accurately characterized by X-ray diffraction; cell parameters were extracted from full pattern matching refinements.

A clear correlation between crystallite size and cell parameters, that is, an asymptotic decrease of a and c cell parameters versus the crystallite size, was established. A simple model was also successfully developed to interpret the as-established correlation.

Introduction

Nanocrystals (NCs) are of growing interest in material science as well as for numerous industrial applications. A type of synthesis that turned out to be well-suited for the preparation of spherical polycrystalline particles (30-500 nm in size) made of nanosized building units (5-50 nm) is the so-called polyol route. It has been initiated by the Fievet and colleagues (1)(2)(3) in the 1990s and then largely developed by Feldmann et al. (4)(5)(6)(7)(8) The polyol route was especially used for the elaboration of ZnO NCs with high purity/homogeneity and tunable particle morphology. (9) In a first step, a forced hydrolysis of zinc precursor (typically hydrated zinc acetate) is performed at temperatures of about 150-160 °C in a multivalent and high-boiling alcohol, for example, diethylene glycol (DEG, bp 246 °C). A second condensation step leads to the formation of ZnO NCs aggregated to form polycrystalline spheres. With different metallic precursors, Tseng et al. have recently obtained various shapes of the polycrystalline aggregates, for example, spheres and rods. (10) Several synthesis parameters can be varied: the composition of the polyol medium (ethylene glycol, glycerol, tetraethylene glycol), the zinc precursor source, its concentration, and its hydrolysis rate. (11)(12)(13)(14) Most of the literature discussions dealt with the polycrystalline size, but not with the crystallite unit size. In one of our previous studies, it has been shown that polyol-mediated synthesis is also adequate to obatine A 3+ -doped ZnO (A 3+ = Al 3+ , Ga 3+ , or In 3+ ) particles for applications as TCO (transparent conductive oxides) thin films or a thermal insulation barrier. (15) On the other hand, correlation between the crystallite size of inorganic oxides and the cell parameters, even for textbook structures, was found to be poorly reported in literature. The unit cells of most NCs are significantly different from those of the coarse-grain or single-crystal counterparts. Even if the phenomenon is well-known, it has been sparingly reported as the consequence of grain-surface relaxation (GSR effect). Stoichiometric effects related to vacancies or interstitial atoms are generally cited. (16)(17)(18)(19)(20)(21)(22) connected to the displacement of atoms next to the surface of a grain from the equilibrium position is not usually considered even if this effect seems to be important. (23) In this study, ZnO NCs were prepared from polyol-mediated synthesis using zinc acetate as precursor. The zinc precursor concentration and the composition of the DEG/hexanol mixtures were independently varied. The particle morphologies were investigated using TEM microscopy, and the crystallite size and cell parameters were extracted from X-ray diffraction analyses. Considering a GSR effect connected to the displacement of atoms next to the grain surface due to the decrease of the surface atom coordination number (bond cut), a model based on the bond valence model [START_REF] Brown | The Chemical Bond in Inorganic Chemistry: The Bond Valence Model[END_REF][START_REF] Brown | [END_REF] was developed to explain the observed correlation between the unit cell parameters and the crystallite dimension.

Experimental Section

Suspensions of Highly Dispersed ZnO NCs

ZnO colloids were obtained through polyol-mediated precipitation, that is, forced hydrolysis at 170 °C of zinc acetate in DEG. In a typical procedure, a known amount of zinc acetate dihydrate (Sigma, ≥99%) was first dissolved in a known volume of DEG (Sigma, 99%). After heating up the solution to 160 °C for 1 h under stirring in a Pyrex round-bottom flask fitted with a reflux column, the precipitation of ZnO crystallites occurred. The solid phase was recovered by centrifugation and washed with ethanol.

The effects on the morphology of the metal cation concentration and of the chemical composition of the alcoholic medium were investigated. Thus, syntheses were performed, on one side, with various Zn 2+ concentrations: 0.05, 0.1, 0.5, and 0.9 M in pure DEG, and, on the other side, with a fixed Zn 2+ concentration (0.1 M) in various DEG/hexanol mixtures, with 1-hexanol (Sigma, 99%) fractions of 0, 50, 75, 90, and 100 vol %.

Characterization of the ZnO Particles/Suspensions Electron Microscopy

The size and the morphology of the so-obtained particles were evaluated by scanning electron microscopy (SEM) and transmission electron microscopy (TEM). The SEM micrographs were recorded with a HITACHI 4500-I apparatus fit out with a field emission gun (SEM-FEG) working at 3.0 kV. Such equipment allows a high spatial resolution (roughly equal to 5 nm) while operating at such low accelerating voltages. TEM experiments were performed on a HITACHI H7650 equipment with a field emissive gun operating at 80 kV and a point resolution of 1 nm.

X-ray Diffraction

X-ray diffraction patterns were collected using a PANalytical X'Pert MPD PW 3040 with Cu Kα1 radiation (λ = 1.540562 Å), and the measurements were performed in continuous scan mode with a scan step of 0.028 in the 2θ range from 20 to 80° (overnight measurements). The setup of the instrument is of Bragg-Bentano geometry with a primary monochromator Ge(111) equipped with fixed divergence and antiscattering slits. The Fullprof program was used for both full pattern matching and structural refinements by the Rietveld method. (26,27) Nevertheless, the atomic position refinements bringing no supplementary information, only the results and the conventional reliability factors (cRp and cRwp) corresponding to the full profile matching mode are herein reported. X-ray diffraction refinements using the peak profile Fullprof function n°7 (Thomson, Cox, and Hastings peak profile function) were also used to evaluate the size of the crystallites. IRF (instrumental resolution function) was defined using very similar collecting parameters than for the samples on the LaB6 "SRM660a" standard material. The shape of the crystallites deduced the application of an anisotropic size-strain model (platelet-type:

-1 or needle-type: +1 model), which was performed empirically, that is, testing one by one different anisotropic axes, and finally selecting the shape-type and the deformation axis, which lead to the smallest reliability factors (the best fit). With Fullprof software, the parameter relative to the anisotropy intensity is the "LorSiz". Size-Model (xxx) ± 1 means: the deformation axis is the (xxx) axis (a "flattening" axis for platelet-type model: +1 or an elongation axis for needle-type shape:

-1.)

Results and Discussion

Effect of the Zinc Cation Concentration TEM and SEM images on the primary particles obtained using different Zn 2+ concentrations in DEG show a drastic evolution of morphology (Figure 1). Spherical polycrystalline aggregates were obtained in the ethanol suspensions for low Zn 2+ concentrations, for example, 0.05 and 0.1 M (Figure 1a,d). Statistical analyses performed on SEM images showed an average diameter of about 310 nm with a quite narrow size distribution (the particles ranging between 150 nm for the 10% smallest ones and 480 nm for the 10% biggest ones). Increasing the ZnO precursor concentration up to 0.5 M induced two important changes: (i) the average diameter of the spherical aggregates increased up to 800 nm, and (ii) "softly agglomerated" crystallites appeared besides the aggregates (Figure 1b,e). The latter phenomenon was drastically enhanced when increasing once again the precursor concentration up to 0.9 M, giving agglomerate-free dispersions (Figure 1c,f). The increase of the precursor concentration induced a significant decrease of the spherical agglomerates' cohesion strength until the complete "dispersion" of the constituting crystallites in the ethanol suspensions. The obtaining of less coherent agglomerates with the metal cation concentration increase was actually due to an increase of the constituting crystallite size (Table 1). First, it has to be noted that the crystallites can exhibit anisotropic shapes in resonance with the würtzite structure, needing hence to make pattern refinements with anisotropic peak profile function. As illustration, the pattern refinement of the 0.1 M sample using an isotropic or anisotropic model along the (110) axis is reported in Figure 2 (experimental, calculated and differential spectra are both represented). The zoom-in on some characteristic 2θ angle range of the differential spectra illustrates the importance of getting correct refinement quality and low reliability factors to perform the refinements with an anisotropic peak profile model. The crystallite sizes as well as the crystalline anisotropy are both deeply impacted by the precursor concentration. As shown in Table 1, the anisotropic shape of the crystallite changed dramatically from (110) needle to (110) platelet; that is, the crystallite exhibited an opposite morphology, when the precursor concentration varied from 0.1 to 0.5 M. Additionally, the crystallite size increased from 7 to 21 nm along the (110) axis and continued to expand while a higher concentration (0.9 M) was used. The progressive increase of the crystallite size with the precursor concentration was also evidenced from the progressive reduction of the peak width (Figure 3). The concentration of metal cation during the precipitation controls the thermodynamic balance between the nucleation of the crystallites, the growth of the pre-existing crystallites, and the aggregation of the crystallites. Here, the precipitation performed in concentrated medium favored the crystallite growth. It can be supposed that, even if the increase of the concentration, and so the supersaturation reached before the germination step, tends to increase the nucleation velocity (La Mer model (28,29) ), the ratio between the number of nuclei and the available Zn 2+ cations in solution decreases with increasing Zn 2+ precursor concentration, leading thus to larger crystallites at the end of the nucleation-growth process. The increase of the crystallite size induces a collapse of the van der Walls interactions between crystallites responsible of the overall cohesion of the agglomerates, the interface (contact) area (between two crystallites in the spherical aggregates) per volume unit becoming smaller and smaller. Thus, suspensions of highly dispersed NCs can be readily prepared from a high Zn 2+ concentration in DEG.

Effect of the Alcohol Nature

This study was conducted in a symmetric way to that related to the influence of the precursor concentration. TEM-extracted morphologies obtained for various hexanol concentrations in DEG are shown in Figure 4. In parallel, the crystallite size was calculated from X-ray pattern refinements; diffractograms were superimposed (Figure 5). The refinement results are reported in Table 2. Briefly, the replacement of DEG by hexanol had the same effect as the increase of the precursor concentration in DEG, that is, an increase of the crystallite sizes leading to spherical aggregates less and less cohesive until the achievement of dispersion at the NCs level. Up to a hexanol concentration equal to 50 vol %, the net increase of the crystallite size remained low: from about 7 nm for the synthesis performed in pure DEG to 11 nm for the 50-50 vol % DEG-hexanol medium.

The spherical aggregates' morphology was also weakly impacted; a slight growth of the aggregates, positively correlated to the constituting crystallite size, can be nevertheless evidenced between panels (a) and (b) in Figure 4. The X-ray diffraction refinements showed that, with 75 or 90 vol % of hexanol, the crystallite size increased to about 15 nm or about 25 nm, respectively. The agglomerates had disappeared to give way to crystallites that were found isolated or grouped to form small agglomerates of few tenths of a nanometer (Figure 4c). It can be noted that the crystallite shape tended to be changed as a function of the polyol medium used, as it was previously shown as a function of the metal precursor concentration.

Here, it is obvious that the parameter directly governing the nucleation-growth thermodynamics is the medium permittivity (εr): 31.7 and 13.3 for DEG and hexanol, respectively. Whatever, the obtaining of less cohesive agglomerates with the higher hexanol concentrations is due, like in the previous experiments, to a decrease of the van der Walls interactions related to an increase of the average crystallite size.

Correlation between Crystallite Size and Cell Parameter

By varying the key parameters, that is, the composition of the alcohol medium or concentration of the metal precursor, the polyol-mediated syntheses allowed controlling the crystallite size and shape in the nanometer range. The pattern matching of the X-ray diffraction spectra using the (110) anisotropic shape allowed independently and accurately the characterization of both a and c cell parameters. Bringing together the results obtained in the two previous studies, the evolution of the c and a cell parameters versus the crystallite size along (002) and (100) axes are plotted in Figure 6a,b, respectively. Both a and c cell parameters decreased versus the crystallite size with an evolution that can be described as a clear and monotonic asymptotic curve. A model from solid-state chemistry as the bond-valence model of Brown et al. [START_REF] Brown | The Chemical Bond in Inorganic Chemistry: The Bond Valence Model[END_REF][START_REF] Brown | [END_REF] can help to interpret A Pauling-Brown bond relaxation approach was attempted to try to explain the variation of the cell parameters versus the crystallite size. In this approach, in a first step, we tried to identify "natural" surface planes considering that the most easy and only acceptable way (on a thermodynamic point of view) to cut the tetrahedral sites is to cut only one bond on its four ones, that is, to keep three neighbors for the O 2-or Zn 2+ surface cations. From this hypothesis, it is easy to identify the (010) or (100) planes, as represented in Figure 7, as natural cutting planes. A second step is to attribute the bond valence of the successive bond layers from the surface to the bulk, layer by layer, in order to respect the ions' V valence value (V = 2).

Hence, with the cutting plane considered, for the "surface" bonds represented in bold black in From eq 1, the bvi values are inversely proportional to the bond length. This equation was used to calculate the bond length of the successive half-layers from the surface (Figure 8a). In accordance with the calculated bond lengths, some shortenings or some increase of the interlayer distances have to occur; the relaxation displacement of each layer is represented by arrows in Figure 7. Actually, the half-layers are alternatively thicker or thinner. In a next step, it is possible to calculate each layer thickness, from the sum of the two corresponding half-layers, adding up the projection along the axis perpendicular to the surface plane (Figure 8b,c). With the consideration of an almost regular tetrahedral shape for each site (with bond angles of 109°), that is, not taking into account the relaxation atomic displacements, the layer thickness, for each layer, can be derived from eq 2 where bl1 and bl2 are the bond length of the first half-layer (bonds represented in Figure 7 in deep colors) and the bond length of the second half-layer (in light colors), respectively.

From this mathematical treatment, the layer thickness versus their number is represented in Figure 8b, showing the sum for each layer thickness of both half-layer components. A zoom-in on the top of the histogram curve clearly shows the asymptotic decrease of the layer thickness from the surface to the core of the crystallites (Figure 8c). Considering that the cutting plane perpendicular to each cell axis can be modeled like the cutting plane represented in Figure 7, a last mathematical treatment allowed calculating the c or a average parameters versus the crystallite size from the thickness of the layers depending on their position from the surface. The evolution of the elongation factor described as cexp/cbulk versus the crystallite size brought to the fore a similar evolution: an asymptotic decrease extended on the same crystallite size range as the one extracted experimentally from the X-ray pattern refinements (Figure 9). It seems that the "Pauling-Brown" approach used to simulate the atomic displacements issued from the bond missing on the crystallite surface is appropriate enough to predict the dependence of cell parameters versus crystallite size in the würtzite structure.

Actually, the curve shape is perfectly predicted as shown in the Figure 9 curve; nevertheless, the a and c numerical values are slightly overestimated by our model. An explanation could be that no pending bonds are taken into account in our model, whereas the zinc oxide can interact with the gaseous atmosphere, forming pending bonds (with hydroxyl or carbonyl groups). Consequently, a decrease of some surface bond valences in comparison with our model can occur due to the zinc oxide and atmosphere interaction.

Conclusions

We reported on a simple and cheap method for the elaboration of ZnO NCs with crystallite size as well as the aggregation behavior that can be controlled. The method already known as polyol-mediated synthesis was here optimized by the systematic study of the zinc precursor concentration and the composition of the DEG/hexanol mixture chosen for the forced P a g e | 10

hydrocondensation step. The lattice parameter values showed a significant size effect: an exponential increase of both a and c cell parameters while the crystallite size decreases was observed. This was interpreted as a grain-surface relaxation effect in nanocrystals. On the basis of the bond valence law, a simple model was successfully proposed, leading to a satisfying fit of the experimental cell parameters curve versus crystallite size. It was shown that, even if the evolution of cell parameters evidenced in this paper is obviously not universal and depended on the nature of the surface plane, the compound crystallographic structure, the occurrence or not of hydrated water or hydroxyl group on the surface, etc., the Pauling-Brown approach developed here can be applied in a universal way to predict the evolution of cell parameters versus crystallite size.

P a g e | 3

 3

Figure 1 .

 1 Figure 1. Typical SEM and TEM micrographs of ZnO particles showing the spherical aggregates and dispersed nanocrystallites, as a function of the zinc precursor concentration used in the polyol-mediated synthesis in pure DEG: 0.1 M (a, d), 0.5 M (b, e), and 0.9 M (c, f).

P a g e | 4

 4

Figure 2 .

 2 Figure 2. Experimental, calculated and differential XRD patterns illustrating the pattern matching quality for ZnO compounds obtained using a Zn precursor concentration in DEG equal to 0.1 M: isotropic peak profile model (a) and anisotropic (110 vector) peak profile model (b).

Figure 3 .

 3 Figure 3. XRD patterns of ZnO particles as a function of Zn precursor concentration in pure DEG.

Figure 4 .

 4 Figure 4. Typical TEM micrographs of ZnO particles showing the spherical aggregates and dispersed nanocrystallites, as a function of the fraction of hexanol in the DEG/hexanol mixture used as synthesis medium for a Zn precursor concentration of 0.1 M: 0 (a), 50 (b), 75 (c), and 90 vol % (d).

Figure 5 .

 5 Figure 5. XRD patterns of the ZnO particles as a function of the fraction of hexanol in the DEG/hexanol mixture used as synthesis medium for a Zn precursor concentration of 0.1 M.

P a g e | 7

 7 this behavior. It can be noted that, from knowledge, such a clear observation was never reported. The modification of the ZnO average cell parameters depending on the crystallite size obviously comes from drastic modifications of the interatomic distances between the planes near the crystallite surface. One surface can be described as an area of missing bonds (with the subsequent excess of energy), and intuitively, the missing bonds around the surface atoms should cause a shortening of the bonds between these surface atoms and the atoms of the first underlayment. Indeed, the bond valence law expresses that the distance and the coordination number of an ion are linked by the equation where V is the valence of the ion (equal to the oxidation degree modulus: here, V = 2 for Zn 2+ and O 2-), bvi is the bond valence of each bond implying the central ion, r0 is an empirically determined parameter characteristic of the A-X bond (r0 is tabulated equal to 1.704 Å for Zn 2+ -O 2-bond), and Li is the length of each i bond of the first coordination sphere. The valence is then a function of the number of bond lengths; when the coordination number i decreases, each bond length should also decrease to maintain the V value to 2 for each ion of the ZnO würtzite structure.

Figure 6 .

 6 Figure 6. Cell parameters extracted from XRD pattern refinements for different ZnO compounds versus their grain diameter.

Figure 7 ,

 7 due to the missing of one Zn 2+ for each surface O 2-, and one O2-for each surface Zn 2+ , all bond valences (bvi) have to equal 2/3. The second half-layer of ions gets four neighbors: two of them are surface ions, that is, exchanging with the considered ion a flux of anP a g e | 8electronic charge of 2/3; then the two other bonds have to get a valence of 1/3. This second sublayer of bonds with bv = 1/3 is represented in light gray. The reasoning was extended to the successive sublayers from the surface to the core of the particles and led to the bond valence values reported in Figure7. For example, for the second layer with bonds in blue, the first half-layer (bonds in deep blue) gets a bond with bv = 0.555 and the second half-layer (bonds in light blue) gets a bond with bv = 0.444. Very quickly, the calculated bv values are converging at the bulk bv value of 2, demonstrating the validity of the applied model.

Figure 7 .

 7 Figure 7. Calculation of the bond valences (bv) for bonds of successive layers from the surface to the bulk considering, as an illustration, the section plane perpendicular to the b axis.

Figure 8 .

 8 Figure 8. Calculation of the bond length and layer thickness (from bond projection on the vertical axis to the surface plane) for successive layers: the layers are numbered in growing order (from 0 to n) from the surface layer to the bulk. In (a) can be distinguished the variation of each half-layer, whereas, in (b) and (c), the half-layers are cumulated in order to draw the evolution of the successive layers.

Figure 9 .

 9 Figure 9. Expansion coefficient of the cell parameter perpendicular to the section plane calculated considering the bond valences forsuccessive planes from the surface to the bulk versus grain diameter.

Acknowledgment

The Conseil Régional d'Aquitaine and the Agence Nationale pour la Recherche (project ANR-A2010 JCJC 0805 01) are gratefully acknowledged.

Author Information

Notes

The authors declare no competing financial interest.