Boly Seck
email: seck.boly@ugb.edu.sn

Pierre-Louis Cayrel
email: pierre.louis.cayrel@univ-st-etienne.fr

Idy Diop

Vlad-Florin Dragoi
email: vlad.dragoi@uav.ro

Kalen Couzon
email: kalen.couzon@ens.uvsq.fr

Brice Colombier

Vincent Grosso

Key-Recovery by Side-Channel Information on the Matrix-Vector Product in Code-Based Cryptosystems

Keywords: Code-Based Cryptography, Side-Channel Attack, Matrix-Vector Product Problem, NIST PQC Standardization

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In recent years, a lot of research has been done on quantum computers [START_REF] Larsen | Deterministic multi-mode gates on a scalable photonic quantum computing platform[END_REF][START_REF] Takeda | Toward large-scale fault-tolerant universal photonic quantum computing[END_REF][START_REF] Gyongyosi | A survey on quantum computing technology[END_REF]. These are computers that exploit the phenomena of quantum mechanics to solve difficult mathematical problems in number theory, such as the Integer Factorization Problem or the Discrete Logarithm Problem. Shor [START_REF] Shor | Algorithms for quantum computation: discrete logarithms and factoring[END_REF] proved that if large-scale quantum computers are built, they will be able to break most of the current asymmetric cryptography like RSA, ECDSA, and ECDH schemes. This would seriously compromise the confidentiality and integrity of all digital communications.

Since then, cryptographic community proposed alternative solutions which remain safe in the quantum era. These schemes are called post-quantum resistant. In 2016, the National Institute of Standards and Technology (NIST) made a call to the community to propose post-quantum secure solutions for standardization. The process consists of several rounds, and only some of the candidates in each round are chosen to enter the next round. The most popular approaches are those based on the search for low-weight words for lattice, the problem of decoding random codes, solving multivariate polynomial systems, isogenies, and hash functions [START_REF] Buchmann | Post-quantum cryptography: state of the art[END_REF][START_REF] Bernstein | Post-quantum cryptography[END_REF]. Lattice-based cryptography has the reputation of being very efficient. Code-based cryptography using some codes is often considered to be already more mature and reliable such as McEliece [START_REF] Mceliece | A public-key cryptosystem based on algebraic[END_REF] and Niederreiter [START_REF] Niederreiter | Knapsack-type cryptosystems and algebraic coding theory[END_REF] cryptosystems.

The majority of code-based post-quantum cryptosystems base their security on the classic hard problem in code-based cryptography: the binary Syndrome Decoding Problem (SDP). Informally, for a binary linear code C of length n and dimension k, having a parity-check matrix H, the SDP is defined as follows: given s ∈ F n-k 2 , find a binary vector x having less than t values equal to one, such that Hx = s. The best algorithm to solve this problem in this original version is the Information Set Decoding (ISD) proposed by Prange [START_REF] Prange | The use of information sets in decoding cyclic codes[END_REF]. The ISD techniques are considered the best strategy for message recovery. It consists, in randomly permuting the matrix H (denote P such a permutation) until the support of the permuted x is included in the set {0, . . . , n -k -1}, i.e., the set where the HP is in upper triangular form. It has been considerably improved since then [START_REF] Becker | Decoding Random Binary Linear Codes in 2 n/20 : How 1 + 1 = 0 Improves Information Set Decoding[END_REF][START_REF] Lee | An Observation on the Security of McEliece's Public-Key Cryptosystem[END_REF][START_REF] Leon | A probabilistic algorithm for computing minimum weights of large error-correcting codes[END_REF][START_REF] May | Decoding Random Linear Codes in ˜O(2 0.054n )[END_REF][START_REF] May | On Computing Nearest Neighbors with Applications to Decoding of Binary Linear Codes[END_REF][START_REF] Stern | A method for finding codewords of small weight[END_REF], although the complexity remains exponential in t.

A recent possible solution to solve the syndrome decoding problem is to use Integer Linear Programming (ILP). The idea was first proposed by Feldman [START_REF] Feldman | Decoding error-correcting codes via linear programming[END_REF] and later improved by Feldman et al. [START_REF] Feldman | Using linear programming to Decode Binary linear codes[END_REF]. Since the initial problem is nonlinear, some relaxation was proposed in order to decrease the complexity. Most recently, Cayrel et al. [START_REF] Cayrel | Messagerecovery laser fault injection attack on the classic McEliece cryptosystem[END_REF] showed that the SDP becomes considerably easier to solve if the syndrome is computed in N. They perform a laser fault injection attack on the matrix-vector product when computing the syndrome in the encapsulation of Classic McEliece. This allows them by corrupting some specific instructions during this operation to have a syndrome in N. To solve the syndrome decoding problem in N, they propose to define the SDP as an ILP inspired by the ideas of Tanatmis et al. [START_REF] Tanatmis | A separation algorithm for improved LP-decoding of linear block codes[END_REF]. The complexity of recovering the secret message from the faulty syndrome is polynomial O(n 2 ) with an optimized version of their algorithm.

Afterwards, Colombier et al. [START_REF] Colombier | Profiled Sidechannel Attack on Cryptosystems based on the Binary Syndrome Decoding Problem[END_REF] proposed to perform a message-recovery attack in Classic McEliece that relies on side-channel information only instead of laser fault injection in the previous work [START_REF] Cayrel | Messagerecovery laser fault injection attack on the classic McEliece cryptosystem[END_REF]. The latter depends on the very strong attacker model and does not apply to optimized implementations of the algorithm that make optimal usage of the machine words capacity. Improvements include the power consumption analysis that is sufficient to obtain an integer syndrome using machine learning techniques. To recover the secret message they use the computationally-efficient score function and known information-set decoding methods.

Contribution:

In this work, a key-recovery chosen-ciphertext attack against codebased cryptosystems is performed. We analyze in particular the secret operation of matrix-vector multiplication in Niederreiter decryption using a physical attack. First, we will introduce a new formalism in code-based cryptography. Informally, for z in F n-k 2 of any weight, the Matrix-Vector Product Problem (MVPP) is defined as follows: given

z * in N n-k , find S ∈ F (n-k)×(n-k) 2
such that Sz T = z * . To get z * in N n-k , we will use the same method of the power analysis attack in [START_REF] Colombier | Profiled Sidechannel Attack on Cryptosystems based on the Binary Syndrome Decoding Problem[END_REF]. This method is based on side-channel analysis using random forests to recover z * from the Hamming weight information obtained from the matrix-vector product in the first step of Niederreiter decryption. Second, we show that if we can construct the matrix

Z * = (z * 1 , • • • , z * n-k
) correctly, we can directly find the secret of the cryptosystem. We obtain directly the secret without solving the syndrome decoding problem unlike in [START_REF] Cayrel | Messagerecovery laser fault injection attack on the classic McEliece cryptosystem[END_REF][START_REF] Colombier | Profiled Sidechannel Attack on Cryptosystems based on the Binary Syndrome Decoding Problem[END_REF] and this is applicable for most of the code-based cryptosystems.

Organization: The paper is organized as follows. In Section 2, we focus on codebased cryptosystems, and in particular on the results of the NIST PQC competition. Section 3 defines the new formalism in code-based cryptography, the Matrix-Vector Product Problem (MVPP). In Section 4, we present our attack on the matrix-vector product in Niederreiter decryption using a side-channel attack. Finally, we conclude this paper in Section 5.

Code-Based Cryptosystems

Encoding Theory

Notations The following conventions and notations are used. A finite field is denoted by F, and the ring of integers by N. Vectors (column vectors) and matrices are written in bold, e.g., a binary vector of length n is

x ∈ {0, 1} n , an m × n integer matrix is A = (a i,j ) 0≤i≤m-1 0≤j≤n-1 ∈ M m,n (N). A[i] denotes the i-th
line of A and a row sub-matrix of A indexed by a set I ⊆ {0, . . . , m -1} is denoted by A I, = (a i,j ) i∈I 0≤j≤n-1

. The same applies to column vectors, i.e., x I is the sub-vector induced by the set I on x.

Error-Correcting Codes An [n, k] linear code C over F q is a vector subspace of F n q , where k, n are positive integers with k < n. The elements of C are called codewords. The support of a codeword Supp(c) is the set of non-zero positions of c. We will represent a code either by its generator matrix G ∈ F k×n q such that its lines form a basis of the vector space C, or by its parity-check matrix,

H ∈ F (n-k)×n q
, where HG T = 0 holds. One of the key elements of decoding is the use of metrics. In the Hamming metric, we consider codes with coefficients in F q (generally, F 2 ).

Definition 1 (Hamming metric). Let x ∈ F 2 , the Hamming weight wt(x) is the number of non null coordinates in x, and the distance between two vectors x and y is the number of non null coordinates in wt(xy).

The hardness of general decoding for a linear code is an N P-complete problem in coding theory [START_REF] Berlekamp | On the inherent intractability of certain coding problems (Corresp.)[END_REF]. This is the Syndrome Decoding Problem (SDP), which is the hard problem in code-based cryptography.

Definition 2 (Binary-SDP). Let H ∈ F (n-k)×n 2 , a vector s ∈ F n-k 2 and t ∈ N.
The syndrome decoding problem is to find x ∈ F n 2 such that Hx T = s and wt(x) ≤ t.

The best-known algorithms for solving this problem are all exponential in t. Except if the syndrome is computed in N instead of F 2 [START_REF] Cayrel | Messagerecovery laser fault injection attack on the classic McEliece cryptosystem[END_REF].

Definition 3 (N-SDP). Let H ∈ M n-k,n (N), with h i,j ∈ {0, 1} for all i, j, a vector s ∈ N n-k and t ∈ N * . The syndrome decoding problem in N is to find x ∈ N n with x i ∈ {0, 1} such that Hx T = s and wt(x) ≤ t.
H and x are sampled in the same way as for the binary SDP, only the matrixvector multiplication operation changes, and thus its result s.

Thus we define the new problem on the matrix-vector product as follows,

Definition 4 (Binay-Matrix-Vector Product Problem (MVPP)). Let z ∈ F n-k 2
of any weight, a vector z * ∈ N n-k . The matrix-vector product problem is to find S ∈ F (n-k)×(n-k) such that Sz T = z * .

We can find z * for side-channel information with power consumption analysis and then with a chosen-ciphertext attack we find S.

NIST PQC Standardization and Result

On July 5, 2022, NIST released the first four winning algorithms from a campaign launched in 2016 to standardize post-quantum cryptographic algorithms. These future standards will be default options for selecting post-quantum algorithms in the majority of security products. Provided that these post-quantum algorithms are also combined with proven classical algorithms through hybrid mechanisms. The main goal of the process started by NIST is to replace three standards that are considered the most vulnerable to quantum attacks, i.e., FIPS 186-4 5 (for digital signatures), NIST SP 800-56A 6 , and NIST SP 800-56B 7 (both for keys establishment in public-key cryptography). For the first round of this competition, 69 candidates met the minimum criteria and the requirements imposed by NIST. 26 out of 69 were announced on January 30, 2019, for moving to the second round. Of these, 17 are public-key encryption and/or key-establishment schemes and 9 are digital signature schemes. In July 2020, NIST started the third round of this process where only seven finalists were admitted (four PKE/KEM and three signature schemes). In addition to the finalists, eight alternate candidates were selected.

The first four algorithms selected are a key establishment algorithm named CRYSTALS-Kyber; and three digital signature algorithms named CRYSTALS-Dilithium, FALCON, and SPHINCS+. The first three of these algorithms are based on structured lattices; the last one, SPHINCS+ is a hash-based signature scheme. These four algorithms will therefore be used as the basis for writing U.S. federal standards. The scope of the NIST announcement is international with strong involvement of the cryptography research community, which will make the future US standards also used as de facto international industry standards. Beside the four winners, an extension of the NIST PQC standardization campaign (4th round) is planned for four key establishment algorithms: BIKE [START_REF] Aguilar Melchor | BIKE: Bit Flipping Key Encapsulation[END_REF], HQC [START_REF] Melchor | Hamming Quasi-Cyclic (HQC)[END_REF], Classic McEliece [START_REF] Albrecht | Classic McEliece[END_REF] (all three based on error-correcting codes), and SIKE [START_REF] Jao | SIKE-Supersingular isogeny key encapsulation[END_REF] (isogeny graphs-based). Classic McEliece was the first selected finalist as a key encapsulation mechanism, while BIKE and HQC were alternative candidates. The latter two use special codes to reduce the size of the public key, which is considered the main drawback of code-based cryptosystems.

Classic McEliece is a code-based scheme using binary Goppa codes, the same codes that McEliece originally proposed in [START_REF] Mceliece | A public-key cryptosystem based on algebraic[END_REF]. During Round 2 the scheme merged with NTS-KEM, which was using the same codes. The Classic McEliece scheme uses the dual of McEliece's scheme, as proposed by Niederreiter [START_REF] Niederreiter | Knapsack-type cryptosystems and algebraic coding theory[END_REF], and tightly turns this OW-CPA PKE into an IND-CCA2 KEM. BIKE (Bit Flipping Key Encapsulation) is a key encapsulation mechanism (KEM) based on quasi-cyclic codes with moderate density parity check matrices. The code structure in BIKE is public and allows to reduce the size of the public key. Bit flipping corrects errors by repeatedly flipping the input bits that, given the secret moderate-density parity checks, seem most likely to be errors.

HQC (Hamming Quasi-Cyclic) uses error-correcting codes built from Reed-Muller and Reed-Solomon. The public key includes a random h and s = x+h•y, where x, y are secretes and small Hamming weights. The ciphertext includes u = r 1 + r 2 • h and v = M + s • r 2 + e, where r 1 , r 2 , e are small Hamming weights and M is a message encoded using an error-correcting code. The receiver computes vu • y = M + s • r 2 + eu • y, which is close to M since x, y, r 1 , r 2 , e are small, and decodes the error-correcting code to recover M. SIKE (Supersingular Isogeny Key Encapsulation) is a key encapsulation mechanism based on the hard problem of pseudo-random walks in supersingular isogeny graphs. SIKE is a relatively new problem in the cryptographic arena and currently undergoing several attacks like its instantiation SIDH (Supersingular Isogeny Diffie-Hellman key exchange protocol) [START_REF] Castryck | An efficient key recovery attack on SIDH (preliminary version)[END_REF][START_REF] Fouotsa | A new adaptive attack on SIDH[END_REF][START_REF] Maino | An attack on SIDH with arbitrary starting curve[END_REF][START_REF] Tasso | Resistance of isogenybased cryptographic implementations to a fault attack[END_REF][START_REF] Zhang | Side-channel analysis and countermeasure design on ARM-based quantum-resistant SIKE[END_REF]. These are key recovery attacks, reduces of the level security, side-channel attacks, and fault injection.

Some of these algorithms could therefore later join the same standardization process as the four algorithms already selected. The final objective of NIST is indeed to be able to standardize a varied range of algorithms to cover a majority of use cases. Most of these constructions based on error-correcting codes use matrix-vector products in the decryption, as in Niederreiter's scheme (Table 1).

The goal of our attack is to find the secrete matrix Q. But first, let's assume that we already have the result of the product Q -1 z in N using the same technique as in [START_REF] Colombier | Profiled Sidechannel Attack on Cryptosystems based on the Binary Syndrome Decoding Problem[END_REF]. Our attack on Niederreiter's decryption is now described. It consists in directly finding the secret matrix Q. In the following, we note

An (n -k) × (n -k) invertible matrix Q Compute H pub = QHP pk = (H pub , t) sk = (Q
S = Q -1 ∈ F (n-k)×(n-k) 2 .
We assume that, we can recover the result of the matrix-vector product Sz in N (z * ∈ N n-k ) with side-channel information as in [START_REF] Colombier | Profiled Sidechannel Attack on Cryptosystems based on the Binary Syndrome Decoding Problem[END_REF].

Profiled side-channel measurements. We performed side-channel measurements during the computation of the product Sz in the Niederreiter decryption implementation. The vector z * is computed as the matrix-vector product Sz = z * . We have recorded a single trace that will be sufficient to form the training set for the profiled attack. This trace is composed of n samples and stored as a row vector. We chose the ciphertexts z i in F (n-k) 2 linearly independent as the inputs of the matrix-vector product algorithm. In addition, we stored a second trace, used as a test set when training the classifier. For both traces, we also stored the Hamming weights of the intermediate value resulting from the matrix-vector product Sz i .

After the traces acquisition, we performed an adequate preprocessing for reducing the dimension (eight dimensions since there are nine possible values for the Hamming weight of a byte) of the data by Linear Discriminant Analysis (LDA) to make it easier to handle by the classifier. We chose the random forest algorithm, used previously for side-channel analysis with good results [START_REF] Hettwer | Applications of machine learning techniques in side-channel attacks: a survey[END_REF], to recover the Hamming weight of z * .

We obtained the Hamming weight of the intermediate value of the product Sz i , we derived the entries of z * in N with 98.65% accuracy. Indeed, the Hamming distance between two consecutive intermediate values is exactly the number of 1's found in the bitwise AND between the row of the matrix S and the byte of the ciphertext z. Computing the value of the integer z * entry is equivalent to counting those ones, which in turn is equivalent to summing the Hamming distances between consecutive intermediate values (the absolute value of the difference of their Hamming weights). In our implementation (n = 6, 960, k = 5, 413 and t = 119), the Hamming distance between two consecutive intermediate values is small and satisfies the condition in [START_REF] Colombier | Profiled Sidechannel Attack on Cryptosystems based on the Binary Syndrome Decoding Problem[END_REF]Equation (3)] to recover the entries of z * in N with good accuracy (82% for Hamming distance).

Course of the attack. We propose a chosen-ciphertext attack that essentially consists of 4 steps:

Step 1. We choose the ciphertexts or vectors z i in F

(n-k) 2

linearly independent. We therefore define a matrix Z = (z 1 , • • • , z n-k ) which is invertible.

Step 2. For each vector z i thanks to the side-channel attack, we have the vector z * i = Sz i in N (in reality we get the Hamming weight for each component). This gives us a new matrix

Z * = (z * 1 , • • • , z * n-k ).
Step 3. We solve a matrix system SZ = Z * with S the unknown matrix in

F (n-k)×(n-k) 2
. Since Z is invertible, we multiply on the right-hand side by its inverse and we obtain

S = Z * Z -1 . (1)
Then we just have to read the entries of the right matrix to get the values of S and thus the secret matrix Q. A toy example is discribed in the Appendix A.

The attack as presented above allows to find the secret matrix directly. However, in Step 2, we can raise two questions:

1. Can we know if the matrix Z * is not correct? 2. If so, how can we correct the errors and find the secret matrix S?

We will discuss question 2 in Section 3.2. For the first point, we assume that in Step 2 we obtain the matrix Z * * instead of Z * . So we have

S ′ = Z * * Z -1 (2)
and

Z * * = Z * + E r (3)
where E r is an error matrix. How to distinguish S from S ′ ? We know that S is an invertible matrix in F 2 . Thus it's enough to look directly at its coefficients and compute its determinant.

We have shown that our attack allows us to directly find the secret matrix S in the case where there is no error in Step 2. Otherwise, we know how to detect it. We have two levels of optimization of this attack either minimize the risk of errors when recovering the matrix Z * or reduce its coefficients modulo an integer number to speed up the computations. We can choose judiciously the z i at step 1, for example, taking z i of low Hamming weights allows having regular words with a "1" for each z i . This can considerably reduce the risk of errors during the acquisition of the traces. Moreover, in this case, we would have Z = I n-k and we obtain the secrete matrix S directly without computing Z -1 . We can also suppose that the victim does not accept to decrypt n -k ciphertexts for example, but with the choice of ciphertexts with low regular weights we avoid this problem.

We will now see how to correct the errors in the matrix Z * and find the correct matrix S.

Error Correction

In this section, we will provide an answer to question 2 and show that we can indeed find the matrix S in some cases. We consider the case where we have E r in the matrix Z * at Step 2, equations ( 2) and (3). We consider two assumptions h1 et h2 about E r :

1. The matrix E r has coefficients 0 or 1, (h1 ). 2. The matrix E r has, at most, a 1 on each row, (h2 ). These two assumptions are not restrictive, we will see that we can deduce the general case and we assume that the error can be controlled to some extent, i.e., Z * * does not differ "too much" from Z * .

According to the above assumptions, there exist two finite sets I and J such that:

E r = (i,j)∈I×J E i,j (4)
with E i,j the square matrix of order n -k where all coefficients are zero, except the one of row i and column j which is 1.

We will need the following lemma:

Lemma 1. Let 1 ≤ a, b ≤ n. Let M = (m i,j
) be a square matrix of order n -k.

then :

E a,b M =                0 • • • 0 0 • • • 0 . . . • • • . . . m b,1 • • • m b,n-k . . . • • • . . . 0 • • • 0 0 • • • 0                ← a-th row
In other words:

[E a,b M ] i,j = 0 if i ̸ = a m b,j if i = a
To find S despite the error in Z * . We have

SZ = Z * + E r = Z * + (i,j)∈I×J E i,j Hence S = Z * Z -1 + E r Z -1 = Z * Z -1 + (i,j)∈I×J (E i,j Z -1 )
From the above we deduce the following theorem:

Theorem 1. For any i ∈ [1, n -k], there exists j ∈ [1, n -k] and ε ∈ {0, 1} such that Z * * Z -1 [i] -εZ -1 [j] is binary and Z * * Z -1 [i] -εZ -1 [j] = S[i].
Proof. Let us suppose |I × J| = 1, let I × J = (a, b). Then we have

S = Z * Z -1 + E a,b Z -1 .
According to Lemma 1, only the line a of E a,b Z -1 is nonzero.

We deduce that for all i ̸ = a,

Z * * Z -1 [i] = S[i].
It is ,therefore, sufficient to take ε = 0 and any j.

According to Lemma 1,

E a,b Z -1 [a] = Z -1 [b]
and so it suffices to take ε = 1 and j = b.

Let (a, b) ∈ I × J, (c, d) ∈ (I × J) \ (a, b) and if |I × J| ≥ 2, the hypothesis (h2 ) implies that c ̸ = a.
According to Lemma 1, we have

E c,d Z -1 [a] = 0. thus Z * Z -1 [a] + (i,j)∈I×J (E i,j Z -1 [a]) = Z * Z -1 [a] + E a,b Z -1 [a].
We thus have, on each line, at most one contribution. We deduce then

1. if i ̸ ∈ I, then Z * * Z -1 [i] = S[i],
it is enough to take ε = 0 and any j.

2. if i ∈ I, then there exists j such that (i, j) ∈ I × J, we have

Z * * Z -1 [i] = S[i] + Z -1 [j]
, so it is enough to take the j given previously and ε = 1.

We deduce that in some cases, it is possible to find S by the following algorithm 1.

Algorithm 1 Finding S with errors in Z * 1: We assume that we can determine the erroneous line(s) of S (for example, S is not binary). 2: Let Li be an erroneous row of S, we subtract from Li the rows of the matrix Z -1 and keep those that are binary. 3: We thus obtain a list of possible candidates for S and for each candidate matrix, we compute its determinant to check its invertibility in F2. 4: In particular, if this list is reduced to one element, we obtain S.

We will now lighten the assumptions by deleting (h2 ). According to Theorem 1, we can deduce the general case of our approach.

Corollary 1. Let i ∈ [1, n-k], let r i be the Hamming weight of the vector E r [i], there exists a sequence (j k ) 1≤k≤ri of distinct pairwise elements and a sequence (ε k ) 1≤k≤ri such that

Z * * Z -1 [i]- Z * * Z -1 [i]- ri k=1 ε k Z -1 [j k ] is binary and Z * * Z -1 [i]- ri k=1 ε k Z -1 [j k ] = S[i].
Proof. We use the same notations as before. We suppose that |I × J| ≥ 2 (if not, we are in the previous case). Let (a, b) ∈ I × J.

We try to count the couples (a, t) with t in J. There are as many as the amount of "1" on the row Er[a], in other words, there are r a couples (a, t).

Let

G a = { (a, t) | t ∈ J} = {(a, j 1 ), (a, j 2 ), • • • , (a, j ra )}.
By Lemma 1 , for all c ̸ = a and d ∈ J, the a-th row of E c,d Z -1 is zero. We deduce that Thus we have the following result :

S[a] = Z * * Z -1 [a] = Z * Z -1 [a] + E r Z -1 [a] = Z * Z -1 [a] + (i,j)∈I×J (E i,j Z -1 )[a] = Z * Z -1 [a] + ra k=1 Z -1 [j k ].
1. If i ̸ ∈ I, then the sequence (ε k ) is null and it is enough to take any (j k ) and two by two distinct. 2. If i ∈ I, considering the set G i defined above and according to what precedes, it is enough to take the sequence (ε k ) constant equal to 1 and the (j k ) given by G i .

In the case where the matrix E r has negative coefficients or is not binary, we can adapt Theorem 1. It is sufficient to allow the sequence ε k to take the value -1 and to be able to subtract (or add) the same row several times. We notice that, unlike in the previous case, we cannot give an algorithm to determine a list of possible candidates for the matrix S. A similar approach as the one presented above would be too expensive. We do not know, a priori, the number of rows to remove from each erroneous row. Although correcting the error in the matrix Z * at Step 2 is theoretically possible, it may be difficult in practice.

Comparison to other Attacks

Recall that the goal of our attack is to find the secret matrix Q in the Niederreiter scheme. It has been shown in various previous works that it is possible to obtain secret information about the decryption in the McEliece scheme (Table 2). Falko Strenzke proposed in [START_REF] Strenzke | A timing attack against the secret permutation in the McEliece PKC[END_REF][START_REF] Strenzke | Timing attacks against the syndrome inversion in code-based cryptosystems[END_REF] two attacks on each of the two calls of the Extended Euclidean Algorithm (EEA) in McEliece decoding with the parameters n = 2, 048 and t = 50. This vulnerability in the Patterson algorithm [START_REF] Patterson | The algebraic decoding of Goppa codes[END_REF] in the error correction phase allows an attacker to gather information about the secret n × n permutation matrix P through a timing side channel.

The first attack [START_REF] Strenzke | A timing attack against the secret permutation in the McEliece PKC[END_REF] targets the second call of the EEA in the Patterson algorithm to determine the polynomials forming the error locator polynomial σ(x). The polynomial of deg(σ(x)) ≤ t consists of two polynomials a(x) and b(x) whose degrees have a direct impact on the number of iterations of the EEA. This variation in the number of iterations implies a difference in the execution times and makes possible a timing attack. The attacker creates ciphertext (chosen-ciphertext attack) using random error vectors with Hamming weight wt = 4 and then lets the decryption routine decrypt the ciphertext. It evaluates whether zero or one iteration occurred in the EEA. If an iteration has occurred, deg(b(x)) = 1, nothing is done, and if there are no iterations, deg(b(x)) = 0, the error vector is added as a new row of a matrix over F 2 having n columns. Each time a row is added, a Gauss-Jordan elimination is performed and the rank is determined. Once the maximum rank is reached (here 2,036), the attack is completed with 7,848,229 ciphertexts. However, such an approach only recovers the secret permutation matrix P when the Hamming weight of error e is small (2 or 4).

The second attack [START_REF] Strenzke | Timing attacks against the syndrome inversion in code-based cryptosystems[END_REF] is based on the vulnerabilities that are present in the inversion of the error syndrome through the extended euclidean algorithm (first call in the Patterson algorithm). Strenzke showed the existence of a timing side channel vulnerability in the syndrome inversion that allows the attacker to gain knowledge of the zero-element of the secret support. It is based on the analysis of the key equation to deduce the relations between the degrees of the polynomials involved in it. As in the first attack, this approach only works for Hamming weights of 2, 4, or 6 of the error vector e to recover the secret permutation matrix P . Despite the improvements in [START_REF] Bucerzan | Improved timing attacks against the secret permutation in the McEliece PKC[END_REF], the main problem with these two previous attacks is the number of cases (depending on the Hamming weight of e) that can be exploited to find the secret. In our attack, we have no constraints on the Hamming weights of the ciphertext (or error vector e) to find the secret matrix Q in the case that we correctly construct the matrix Z * with a random forest. Moreover, we attack the least complex step of Niederreiter decoding (first step). The Niederreiter PKE is slightly different from the McEliece PKE (Table 1 and Table 2). However, here too, an error vector is chosen during the encryption and decryption features, and since these features are the prerequisites of our attack, it is also applicable to McEliece's PKE. Unlike the attack of Strenzke in [START_REF] Strenzke | A timing attack against the secret permutation in the McEliece PKC[END_REF] we only need 2,048 ciphertexts instead of 7,848,229 to find the secret n × n

Conclusion

This article presents a key-recovery attack against the Matrix-Vector Product Problem, which is a new formalism that we introduce in based-code cryptography. We have also shown that with a side-channel attack on this operation, we can recover secret information on based-code cryptosystems without solving the hardness of the binary syndrome decoding problem. In addition to the side-channel information, we performed a chosen-ciphertext attack, which, with careful choice of the ciphertexts, can find the secret matrix without errors. When noise (errors) is present during the attack, we have shown that in some cases it is possible to find the secret matrix. Our attack can be applied to code-based cryptosystems with matrix-vector product operations. By performing the same operation as before, we obtain the following S ′ matrix according to equation 2

A Simple Version of Our Attack

S ′ =     0 1 2 0 1 0 1 1 0     .
This matrix S ′ is not binary, so we deduce that the first row contains a fault. We then apply the algorithm 1 of Section 3.2 to determine the list of possible candidates for the matrix S

S ′ [1] -Z -1 [1] = [1 0 1] S ′ [1] -Z -1 [2] = [-1 2 2] S ′ [1] -Z -1 [3] = [-1 1 3].
Only the first case gives a binary vector, we deduce that S ′ [1] = S[1] + Z -1 [START_REF] Aguilar Melchor | BIKE: Bit Flipping Key Encapsulation[END_REF] and we have directly S. We then obtain the following matrix

S ′ =     1 0 1 1 0 0 1 1 0     .
Here the matrix is binary, moreover its determinant det(S ′ ) = 1. In this case, we cannot detect that the matrix is not correct.

  , H, P ) Encrypt(m, pk) = z Encode m → e \\ error vector of wt(e) = t z = H pub e Decrypt(z, sk) = m Compute z * = Q -1 z = Sz \\ target of our attack z * = HP e e * = Decode(z * , H) Retrieve m from P -1 e * 3 Our Attack and Result 3.1 On the Decryption of Niederreiter

Table 2 :

 2 McEliece PKE scheme KeyGen(n, k, t) = (pk, sk) G-generator matrix of C that corrects t errors An n × n permutation matrix P An k × k invertible matrix S Compute G pub = SGP pk = (G pub , t) sk = (S, G, P ) Encrypt(m, pk) = z Encode m → c = mG pub z = c + e \\ e is an error vector of wt(e) = t Decrypt(z, sk) = m Compute z * = zP -1 z * = mSG + eP -1 m * = Decode(z * , G) Retrieve m from m * S -1

Case 1 :From the equation 1 ,Case 2 :

 112 no errors in Z * Let the matrices Z invertible in F 2 and Z * in N respectively constructed in Step 1 and recovered in Step 2 of our attack. One chooses here n -k = 3, S = Z * Z -1 , we find few errors in Z * We keep the same Z, its inverse and Z * matrices in case 1. According to equation 3, Z * * = Z * + E r , we have

Case 3 :

 3 errors in Z * We consider the error matrix

Table 1 :

 1 Niederreiter PKE scheme KeyGen(n, k, t) = (pk, sk) H-parity-check of C that corrects t errors An n × n permutation matrix P

Table 3 :

 3 Attacks to recover the secret matrix P . The Table3give more details on the comparison with Strenzke's attacks to find the secret matrix P .

	Attack	Hamming weight	Number of	Target
			ciphertexts
	Timing attacks of	2, 4 or 6	7,848,229 in [30]	EEA in Patterson
	Strenzke			algorithm
	Our attack against	no constraints	2,048	First step in
	MVPP			McEliece's
				decryption
	permutation matrix P

https://nvlpubs.nist.gov/nistpubs/FIPS/NIST.FIPS.186-4.pdf

https://nvlpubs.nist.gov/nistpubs/SpecialPublications/NIST.SP. 800-56Ar2.pdf

https://nvlpubs.nist.gov/nistpubs/SpecialPublications/NIST.SP.


Mélanie Debiais-Thibaud
email: melanie.debiais-thibaud@umontpellier.fr

Sylvain Marcellini
email: smarcellini@udec.cl

Samuel A Cushman

Frédéric Marin
email: frederic.marin@u-bourgogne.fr

Editorial: The evolution of biomineralization in metazoans

Keywords: biomineralisation, evolution, metazoan, corals, mollusks, vertebrates, sponges (porifera)

The evolution of biomineralization in metazoans

.

 illustrates the variety of contributions covering a large range of metazoan phyla. Of course, encompassing the whole metazoans was a challenge and some groups are missing in this volume: siliceous sponges for example, but also brachiopods, bryozoans, calcifying annelids, or, among the ecdysozoans, crustaceans. However, the fifteen articles presented here constitute a significant sampling of the research on evolutionary processes that have been driving animal biomineralization for an eon.

Regarding non-bilaterians metazoans, Voigt et al. highlight the evolution of carbonic anhydrases in calcareous sponges (Calcarea), distinguishing the different cellular locations of this enzyme in the two subclasses Calcinea and Calcaronea. For cnidarians, Zaquin et al. show how the skeletal matrix proteins of scleractinian corals result from the independent co-option of ancestral genes predating cnidarian diversification, but also from posterior recruitment of duplicated genes that are species-specific. Finally, Conci et al. review different cellular and molecular aspects of the skeletal biomineralization of another cnidarian clade, Octocorallia,

which can 'handle' the two calcium carbonate polymorphs, i.e., calcite-the predominant form-but also, more rarely, aragonite.

The four papers covering the protostome world deal with the molecular aspects of shell biomineralization in molluscs, namely, three bivalves and one gastropod. Sato et al. emphasize the seasonal plasticity of the shell, by describing the microstructure of the winter layer and its associated molecular markers, in the shell of the Japanese pearl oyster, while McDougall et al. associate changes in gene expression in the pearl sac to the pearl structural characteristics in the south sea pearl oyster Pinctada maxima. Fundamental aspects of shell evolution are examined with the contribution of Takeuchi et al. who identify the complete protein shell repertoire of the Japanese crocus clam. At last, Jackson reveals the modularity and plasticity of the calcifying mantle of the terrestrial grove snail, Cepaea nemoralis, by looking at the spatial expression pattern of some key shell-forming genes.

A total of eight articles relate to biomineralisation in the skeleton of deuterostomes, with a major focus on the derived calciumphosphate based skeleton of vertebrates. In particular, Pears et al. While the Research Topic of these research and review articles demonstrates the diversity of biomineralization forms in animals, they also surprisingly underline major genetic variation at all levels of comparisons, from the deepest nodes of metazoans, to more recent divergent groups within vertebrates, molluscs, corals or sponges. We hope that this Research Topic will inspire and encourage the community to further explore this complex and fascinating aspect of animal biology.

  describe in unprecedented detail the delicate mineralisation in the cartilaginous skeleton of Holocephalans, while Clark et al. uncover the various architectures of highly mineralised jaws in batoids (rays). More general aspects of cartilaginous fish skeletal mineralisation are exposed in the review of Atake and Eames and a research article proposed by Berio et al. Specific aspects of skeletal gene evolution are analysed for two vertebrate gene families: Mgp/Bgp Leurs et al., and SPARC/SPARCL Romero et al., Deuterostomian calcium carbonate-based biomineralizations are not to be outdone, with the contribution of Le Roy et al. on the molecular toolkit for building avian eggshell and with the review of Khor and Ettenson on the skeletal function of Alx homeobox genes in echinoderms and all deuterostomes in general.

Frontiers in Genetics frontiersin.org

Author contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

FIGURE 1

The phylogenetic relationships of the 15 contributions of this Research Topic.

Frontiers in Genetics frontiersin.org

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.


Hrishikesh Pingulkar

Jean-Baptiste Salmon

Confined directional drying of a colloidal dispersion: kinetic modeling

We derive a model to describe the dynamics of confined directional drying of a colloidal dispersion. In such experiments, a dispersion of rigid colloids is confined in a capillary tube or a Hele-shaw cell. Solvent evaporation from the open end accumulates the particles at the tip up to the formation of a porous packing that invades the cell at a rate dl dt . Our model based on a classical description of fluid mechanics and capillary phenomena, predicts different regimes for the growth of the consolidated packing, l versus t. At early times, the evaporation rate is constant and the growth is linear, l ∝ t. At longer times, the evaporation rate decreases and the consolidated packing grows as l ∝ √ t. This slowdown is either related to the recession of the drying interface within the packing thus adding a resistance to evaporation (capillary-limited regime), or to the Kelvin effect which decreases the partial pressure of water at the drying interface (flow-limited regime). We illustrate these results with numerical relations describing hard spheres, showing that these regimes are a priori experimentally observable. Beyond this description of the confined directional drying of colloidal dispersions, our results also highlight the importance of relative humidity control in such experiments.

Introduction

The drying of colloidal dispersions, allowing by simple evaporation of a solvent to form a dense and consolidated material from colloids dispersed in a liquid, is a key step of many industrial processes, from coatings to electrode manufacturing. [START_REF] Keddie | Fundamentals of latex film formation. Processes and Properties[END_REF][START_REF] Li | [END_REF] The detailed understanding of the mechanisms at play remains a scientific challenge and involves multiple aspects: transport phenomena, physical chemistry of colloids, mechanical instabilities, etc. [3][4][START_REF] Goehring | Desiccation Cracks and their Patterns: Formation and Modelling in Science and Nature[END_REF][START_REF] Bacchin | [END_REF] In this context, different studies have used a simplified geometrical configuration known as confined directional drying to provide a better understanding. Figure 1 shows these experiments schematically: a horizontal capillary tube or a Hele-Shaw cell, with a cross-sectional dimension h typically less than 100 µm, is filled by an aqueous colloidal dispersion and open on one side to the ambient atmosphere with relative humidity RH. As long as capillary forces pin the air-water interface at the open end of the cell, evaporation causes a flow at a volumetric rate J (m/s):

J = 1 Aρ w dM dt , (1)
to replenish the lost water mass dM dt (kg/s), A (m 2 ) being the cross-sectional area of the cell and ρ w (kg/m 3 ) the mass density of liquid water. This flow continuously accumulates the colloids a CNRS, Solvay, LOF, UMR 5258, Univ. Bordeaux, F-33600 Pessac, France at the tip until the formation of a consolidated close-packed network. For rigid colloids, the porosity of this layer does not prevent evaporation, and thus its growth at a rate dl dt , see Fig. 1. This configuration shares various similarities with the more general directional drying often observed for liquid films deposited on a substrate, 7,8 but with the notable difference that evaporation occurs here only through the end of the tube, and thus only through the colloidal packing. In this sense, confined directional drying is fundamentally closer to configurations such as drying of quiescent films, 9,10 i.e., without lateral flows induced by edge effects, or drying of 2D and 3D drops, for which evaporation also occurs only through the consolidated colloids. [11][12][13][14] The same remark also applies to the classical case of sessile drops, for which several ingredients presented in Fig. 1 also apply, but with the difference of the finite volume of dispersion and the complexity associated with the possible deformation of the evaporation surface during drying. 15 Confined directional drying has been widely used to study the fracture patterns induced by drying, [16][17][18][19][20][21] the phenomenon of delamination, 22,23 or the appearance of shear bands. 24,25 Beyond these issues concerning mechanical instabilities, this technique has also been used to study mass transport in colloid dispersions, [26][27][28][29] the formation of an amorphous solid prior to consolidation, 30 or the mechanisms of colloidal ordering. 31,32 In spite of all these works, one fundamental question remains unanswered: that of the growth kinetics of the colloidal packing, l as a function of t. Experiments have shown that the growth is linear, l ∝ t, for colloids with large radii (a > 100 nm 33 ), as expected from mass conservation 34,35 for a constant evaporation rate J. For smaller colloids, typically a < 20 nm, a transition towards a square root growth, l ∝ √ t, has often been reported at longer times, suggesting a slowdown of the evaporation rate. 17,18,20,22 Dufresne et al. were the first to propose a law able to fit all the l(t) observations, and explained the transition from l ∝ t to l ∝ √ t, by invoking a transition from an evaporationlimited to a flow-limited regime. 17 Their model also showed that the crossover between these two regimes occurred for a critical length l c scaling linearly with the colloid size l c ∼ a, in conflict with the experimental data. 17 Furthermore, the derivation of this law was empirical, and the same authors later proposed 18 that the square root slowdown was actually probably due to the recession of the air-water interface within the porous colloidal packing over a small length ξ , as shown in Fig. 1(b). Wallenstein and Russel then suggested theoretically that Knudsen diffusion of water vapor through the dry porous layer results in a very strong resistance to mass transfer, leading to a slowing of the evaporation rate J and a square root growth of the front, l ∝ √ t, but with a tiny recession, ξ ≪ l. 36 Later, Lidon et al. showed that the Kelvin effect, i.e., the decrease in the partial pressure of water due to highly curved menisci at the drying interface, see Fig. 1, should probably also be taken into account, especially for colloids with small radii a. 37 Their model, including also the possible recession of the drying interface, now predicted that the crossover length scales as l c ∼ a 2 , possibly explaining the discrepancies between experiments and models pointed out by Dufresne et al. 17 Importantly, their model incorrectly assumed that the ambient relative humidity RH was strictly zero. As shown below, this model cannot describe with this assumption the existence of the flow-limited regime suggested in Ref. 17 This error was recently revealed by Tatsumi et al. who developed a model to describe the drying of quiescent colloidal films. 38 This model, of which some relations are resumed below, does not take into account the possible recession of the interface in the packing, and therefore cannot describe part of the regimes detailed in the present work.

(b)

0 h ξ J r c -ξ 0 Ψ ext x Ψ -l -l l Ψ cap h J (a) r ϕ d ϕ 0 0 l Ψ ext x Ψ Ψ cap RH RH dl dt dl dt
The primary motivation of our work is to clarify the description of the kinetics, l versus t, using a simple model of confined directional drying of a colloidal dispersion. Our model takes into account both the Kelvin effect and the possible recession of the evaporation interface, and applies continuously to colloids with radii ranging from a few nanometers to a few hundreds of nanometers. As shown below, this model highlights the importance of RH on the growth kinetics of the dense packing, and reveals the existence of two distinct regimes to explain the square root growth l ∝ √ t. Beyond the very specific geometry of confined directional drying described in Fig. 1, the ingredients of our model also apply a priori to any drying configuration of colloidal dispersions (including also those of finite volume such as the drying of sessile droplets and films), when the Kelvin effect, but also capillary effects come into play by modifying the evaporation conditions. To our knowledge, these ingredients are rarely taken into account, and the evaporation rate is often described as not depending on the presence of closely-packed colloids.

The present paper is organized as follows. In Sec. 2, we first present the main assumptions of our model and the equations governing the growth rate dl dt . In Sec. 3, we then present the main results of this model including the different growth regimes, and discuss the role played by relative humidity. Finally, we illustrate the model on a relevant experimental case using standard relations valid for monodisperse hard-spheres. In Sec. 4, we conclude our work by suggesting possible experiments and research directions to improve the model.

Model

Geometry, assumptions, and water potential

We consider the configuration shown in Fig. 1(a). The cell is initially filled by a dispersion of colloids with radii a ≪ h, at a volume fraction ϕ 0 . We consider that colloids cannot adhere to the cell walls in the dispersed state. We also assume that colloids are rigid and cannot deform upon accumulation. The volume fraction ϕ d of the dense packing is thus considered uniform and constant. To model the growth rate of the packing, we consider for simplicity a one-dimensional description and isothermal conditions. We also assume the local thermodynamic equilibrium and quasi-static conditions for writing the different fluxes driven by evaporation.

In the following, we use the water potential Ψ (Pa) because it allows a convenient description of water independently of its thermodynamic state: liquid in the dispersion and in the pores of the colloidal packing, to vapor in the ambient air. The latter is defined by:

Ψ = µ -µ 0 (T ) V m , (2)
where µ = µ(P, T ) is the water chemical potential expressed in units of J/mol at pressure P (Pa) and absolute temperature T (K), µ 0 (T ) is its standard value at atmospheric pressure P 0 and sea level, and V m (P, T ) is the liquid water molar volume (V m ≃ 1.805 × 10 -5 m 3 /mol at P 0 and room temperature, T = 22 • C).

In ambient air, far from the outlet, the water potential is given by:

Ψ ext = RT V m log (RH) , (3)
R being the universal gas constant. Because colloidal interactions in the dilute dispersion lead to a negligible water potential compared to the other potentials considered later, [START_REF] Russel | Colloidal dispersions[END_REF] one has Ψ ≃ 0 upstream x = -l. Water evaporation results from the drop of the water potential from Ψ ≃ 0 to Ψ ext , see Fig. 1.

Importantly, we will assume below a classical description of fluid mechanics and capillary phenomena down to nanometersized colloids. This assumption is based on the results of Vincent et al. on evaporation-induced flows through nano-porous media, [START_REF] Vincent | [END_REF] who showed that a classical description remains valid even with pore radii r p ≃ 1.7 nm. Furthermore, the problem studied by Vincent et al. shares many similarities with the case shown in Fig. 1, but in a static configuration (the size of the porous medium being fixed), and several equations described below have already been discussed in Ref. [START_REF] Vincent | [END_REF]

Evaporation rate

Two configurations emerge depending on the strength of the capillary forces: fully wetted colloidal packing, Fig. 1(a), or recession of the air-water interface, Fig. 1(b).

No recession

In Fig. 1(a), capillary forces are strong enough to avoid the recession of the air-water menisci in the porous network. In this regime, the global mass balance imposes a flow in the cell and through the whole porous packing. This flow depends on two resistances in series for the drop in water potential: vapor removal in ambient air and liquid water flow through the porous layer. Evaporation rate in ambient air is described assuming quasi-steady conditions by: 41

J = kc sat ρ w (RH ⋆ -RH) , (4)
where RH ⋆ is the relative humidity at the evaporation interface (x = 0) and c sat (kg/m 3 ) the water concentration at saturation in air. k (m/s) is a mass transfer coefficient that accounts for the transport of the vapor by diffusion and convection in air.

The evaporation-induced flow leads to a drop of water potential 4) and ( 8) versus 1 -RH, for l/l k = 10 -2 , 10 -1 , 1, 10, 10 2 , and 10 3 . Inset: zoom for the high values of l/l k . The dotted lines are the approximations Eq. ( 11) for l ≫ l k .

due to the viscous friction given by the Darcy's law:

Ψ ⋆ = - η w l κ J , (5)
where η w (Pa s) is the water viscosity and κ (m 2 ) the permeability of the colloidal packing. We neglect in this description the drop in water potential due to the hydraulic resistance of the cell itself, an assumption fully justified since a ≪ h.

Kelvin's relation, i.e., the continuity of the water potential at the air-water interface, relates RH ⋆ to Ψ ⋆ :

RH ⋆ = exp V m Ψ ⋆ RT , (6)
and the radius of curvature r > 0 of the menisci, see Fig. 1(a), is given by the Kelvin-Laplace equation:

Ψ ⋆ = - 2γ r , (7)
where γ (N/m) is the air-water surface tension. Equations ( 6) and (7) show that RH ⋆ deviates from RH ⋆ = 1 only for water potential drops -Ψ ⋆ that are not negligible compared to RT /V m (≃ 136 MPa at T = 22 • C), and thus radii of curvature typically less than r ≤ 10 nm considering γ ≃ 70 mN/m.

Using Eqs. ( 4), ( 5), and ( 6), one can show that:

RH ⋆ = l k l W l l k exp RH l l k , (8)
with the length scale l k defined by :

l k = κ(RT /V m ) η w kc sat ρ w , (9)
and W the Lambert W -function, 42 the inverse function of f (w) = w exp(w), i.e., z = w exp(w) ⇔ w = W (z).

Figure 2 shows the evaporation rate given by Eq. ( 4) with RH ⋆ calculated using Eq. ( 8) versus (1 -RH) and for various values of l/l k . For l ≪ l k , the Taylor series 42 of W in Eq. ( 8) leads to RH ⋆ ≃ 1, so that evaporation is limited by the rate of vapor removal in the ambient air and thus:

J ≃ kc sat ρ w (1 -RH) . (10)
For l ≫ l k , the asymptotic expansion 42 of W in Eq. (8) shows that RH ⋆ ≃ RH, and Eq. ( 5) leads to:

J ≃ - κ η w l Ψ ext = - kc sat ρ w l k l log(RH) . (11)
In this case, evaporation is limited by the flow through the porous layer and does not depend anymore on the mass transfer coefficient k in ambient air. 40

Meniscus recession

When capillary forces are no longer able to maintain the pinning of the menisci at x = 0, the air-water interface recedes in the porous packing, see Fig. 1(b). This recession occurs when the radii of curvature reach a critical value r c , or equivalently when the drop in water potential Ψ ⋆ reaches the maximal capillary pressure

Ψ cap = -2γ/r c .
In this regime, there is an additional resistance to evaporation due to the transport of the water vapor in the dry porous layer of thickness ξ . To model this case, we first assume that the rate of recession is small compared to the rate of water mass loss dM dt . Equation (1) thus remains valid for estimating the evaporationinduced flow J, and we will discuss later this assumption in more detail.

We also assume homogeneous recession over the cross-section of the porous layer. This corresponds to neglecting the pore size distribution (equivalently, the critical radii r c ) of the colloidal packing which would lead to the existence of a partially saturated layer, due to the preferential emptying of the largest pores. [START_REF] Lehmann | [END_REF] For simplicity, we do not consider this possibility often observed for the drying of porous media, 44,45 and our model actually applies only to cases where the partially saturated layer has negligible width.

In the ambient air, Eq. ( 4) still relates the vapor removal rate to the evaporation-induced flow J. The transport of the water vapor through the dry porous layer is limited by diffusion, and the global mass balance imposes:

J = D p w c sat ρ w ξ (RH cap -RH ⋆ ) , (12)
where D p w is the effective diffusion coefficient of the water vapor in the porous medium, RH ⋆ the relative humidity at x = 0, and RH cap the relative humidity at x = -ξ linked to Ψ cap by the Kelvin relation. Finally, the flow J leads to a drop in water potential across the porous layer given by:

J = - κΨ cap η w (l -ξ ) . ( 13
)
Equations ( 4), (12), and ( 13) can be used to show that:

ξ l = RH cap -RH -(D p w Γ)/(kl) RH cap -RH + Γ , (14)
with the dimensionless parameter:

Γ = - κΨ cap η w D p w c sat ρ w . (15)
Γ compares the flux of liquid water driven by the capillary pressure Ψ cap across a given thickness L of the porous layer [κΨ cap /(η w L)] to the flux of liquid water driven by the diffusion of water vapor over the same thickness [D p w c sat /(ρ w L)]. As shown in Sec. 3.3.2 with numerical applications, Γ ≫ 1 regardless of the colloid size a. Because the numerator in Eq. ( 14) is smaller than RH cap -RH ≤ 1, Γ ≫ 1 implies ξ ≪ l, and Eq. (13) shows that the evaporation rate can be approximated by:

J ≃ - κ η w l Ψ cap . ( 16
)

Growth of the porous layer

The flow J driven by evaporation continuously accumulates colloids to the porous packing, which consequently invades the cell at a given rate. The growth dynamics, l versus t, depends on multiple phenomena such as the transport of the colloids by advection up to the packing, diffusion which opposes the formation of concentration gradients upstream x = -l, but also on the consolidation dynamics of colloids in the packing. Despite this complexity, the growth rate is often described in the literature, 17,20,[36][37][38] by the simplified colloid mass balance:

(ϕ d -ϕ 0 ) dl dt ≃ ϕ 0 J . (17)
This relation neglects, however, any temporal evolution of the colloid concentration profile in the dispersion upstream x = -l.

In Sec. 3.3.1, we show using a simple transport model for the colloidal dispersion, that these variations are small, and that Eq. ( 17) correctly describes the growth rate of the dense packing. Moreover, Eq. ( 17) combined with the weak recession of the air-water menisci, ξ ≪ l (coming from the assumption Γ ≫ 1), leads to dξ dt ≪ J. This confirms that the rate of menisci recession makes a negligible contribution to the mass loss dM dt , and that the evaporation-induced flow J remains correctly given by Eq. ( 1) in this regime.

Dimensionless model

For clarity, we define the dimensionless quantities:

l = l l k , t = t t k , and j = J J 0 , (18)
with:

J 0 = kc sat ρ w , t k = l k J 0 ϕ d -ϕ 0 ϕ 0 , ( 19
)
and l k defined in Eq. ( 9). With these dimensionless variables, our model for the growth dynamics is simply written:

d l dt = j . ( 20
)
In the regime of no recession, i.e., when Ψ ⋆ > Ψ cap , Eqs. ( 4) Fig. 3 Critical length lrec (in log scale) given by Eq. ( 25) in the plane α-RH. The white line given by α =log(RH) delimits the zone of menisci recession, see Eq. (24).

and ( 8) lead to the normalized evaporation rate:

j = 1 l W l exp RH l -RH . (21)
When Ψ ⋆ < Ψ cap , the air/water interface recedes, and j given by Eq. ( 16) is now simply written:

j = α l , (22)
with

α = - Ψ cap RT /V m , (23)
a dimensionless parameter comparing the capillary pressure Ψ cap to RT /V m , the characteristic pressure scale of the Kelvin effect.

Results and Discussion

Condition for recession of the interface

In the kinetic model described in the previous section, the recession of the menisci never occurs when the critical capillary pressure is smaller than the water potential in the ambient air, i.e., when Ψ cap < Ψ ext , see Fig. 1. With the notation of Eq. ( 23), this condition reads:

α > -log(RH) . (24)
When Eq. ( 24) is not fulfilled, the menisci recede within the porous layer when Ψ ⋆ reaches Ψ cap . Equations ( 5) and ( 21) then show that it occurs for the critical length:

lrec = α exp (-α) -RH . ( 25
)
Figure 3 shows lrec in the plane α-RH. The regime for which the menisci never recede is expected only at high relative humidity and for large α values. For standard conditions, this regime corresponds to high values of the capillary pressure Ψ cap , since RT /V m ≃ 136 MPa at T = 22 • C. According to the Kelvin-Laplace 20) RH = 0.5. The dashed lines show the approximation Eq. ( 27) at low l, and the flow-limited regime Eq. ( 28) at larger l. The other curves are given by Eq. ( 30) for α = 2 × 10 -1 , 2 × 10 -2 , 2 × 10 -3 , and 2 × 10 -4 (from orange to brown). equation Eq. ( 7), this corresponds to very small critical radii, typically r c ≃ 1.5 nm for RH = 0.5 assuming the surface tension is γ ≃ 70 mN/m. Such cases are therefore only expected for colloids with radii a ≤ 10-20 nm, as in some experimental studies cited above. 17,18,37 As mentioned previously, Vincent et al. studied the drying-induced flows through a porous silicon layer with pore radius r p ≃ 1.7 nm. [START_REF] Vincent | [END_REF] These experiments share various similarities with the problem studied here, but in a static configuration. The experiments of Vincent et al., corresponding to α ≃ 0.6 with our notations, clearly revealed drying-induced flows given by Eq. ( 11) for high RH, and recession of the menisci at lower RH, leading to capillarity-limited flows, Eq. ( 16).

In the case of larger colloids, typically a ≫ 20 nm, one has α ≪ 1, and Eq. ( 25) turns into real units in:

l rec ≃ - κ η w Ψ cap J 0 (1 -RH) . (26)
As expected, this relation does not involve anymore the characteristic pressure scale RT /V m , and could have been derived without invoking the Kelvin effect.

Growth dynamics

Figure 4 shows the numerical solutions of the ordinary differential equation Eq. ( 20) for RH = 0.5 and different α values, with the initial condition l = 0 at t = 0. When the condition Eq. ( 24) is fulfilled, j is given by Eq. ( 21), while for the other cases, Eq. ( 20) is solved with j given by Eq. ( 21) for l ≤ lrec , and with j given by Eq. ( 22) for l > lrec .

We first focus on the regime of no recession (yellow curve in Fig. 4). For l ≪ 1, the evaporation rate is constant and wellapproximated by Eq. ( 10) ( j ≃ 1 -RH with the dimensionless variables), so that the length of the porous layer grows linearly following:

l ≃ (1 -RH) t . ( 27
)
At longer times scales, the evaporation rate decays due to the decrease in water potential at the evaporation interface. For l ≫ 1, Fig. 5 Case RH = 0.5. code for the relative humidity at the evapinterface (RH ⋆ when no recession, RH cap in the other case). The line is Eq. ( 25) and delimits the regime of menisci recession from that of no recession.

Eq. shows that the rate follows approximately j ≃ -(1/ l) log(RH), leading to the root growth:

l2 ≃ -2 log(RH) t . ( 28
)
The transition from the evaporation-limited regime given by Eq. ( 27) to the flow-limited regime given by ( 28) occurs at the length:

lkel ≃ -2 log(RH) 1 -RH , ( 29
)
with the subscript kel to emphasize that the origin of the square root growth at long times is due to the Kelvin effect. When the condition Eq. ( 24) is not fulfilled, the growth first follows the same l versus t curve up to l = lrec given by Eq. ( 25). At later times, the length of the porous layer grows as:

l2 = l2 rec + 2α(t -trec ) , (30)
and thus:

l2 ≃ 2α t , (31)
for l ≫ lrec . Equation ( 30) has already been derived by Wallenstein and Russel, 36 but without taking into account the Kelvin effect, and thus the possible decrease of the evaporation-driving force for small curvature radii. Later, Lidon et al. found again Eq. ( 30) now taking into account the Kelvin effect, 37 but assuming RH = 0 and thus missing the description of the flow-limited regime described above, as revealed recently by Tatsumi et al. 38 Figure 5 summarizes on a single diagram in the plane l-α, the different asymptotic regimes of growth of the packing for a fixed RH, as well as the values of the relative humidity at the air-water interface. This diagram helps to reveal that the model always predicts a linear growth l ∝ t at early times followed by a square root slowing down l2 ∝ t but with two different origins.

For α <log(RH), the slowing down occurs due to the minute recession of the interface that adds a resistance to evaporation. In this capillarity-limited regime, the kinetics follows with real units:

l 2 ≃ -2 ϕ 0 ϕ d -ϕ 0 κΨ cap η w t , (32)
and does not depend neither on the ambient relative humidity (RH) nor on the external mass transfer (J 0 ). For α >log(RH), the square root slowing down is due to the decrease of the evaporation driving force (Ψ ⋆ ≃ Ψ ext ) leading to a flow-limited regime which follows with real units:

l 2 ≃ -2 ϕ 0 ϕ d -ϕ 0 κΨ ext η w t . ( 33
)
In this regime, the growth rate does not depend again on external mass transfer conditions (J 0 ), but still depends on the imposed relative humidity through Ψ ext .

A model case: monodisperse hard-sphere colloids

To illustrate the results of our model by numerical applications and to validate the assumptions made previously, we now use relations for monodisperse hard-sphere colloids of radius a. For such colloidal interactions, directional drying accumulates the colloids at the tip of the cell up to the random close-packing at volume fraction ϕ d ≃ 0.64. We therefore implicitly assume, for simplicity, that colloids do not crystallize into an ordered network, either because of a slight polydispersity or because the crystallization kinetics is slower than the colloid accumulation rate. For such randomly-packed assembly of hard-spheres, the mean pore radius r p is defined as twice the ratio of pore volume to surface area: 46,47

r p = 2 3 1 -ϕ d ϕ d a, (34)
r p ≃ 0.375 a for ϕ d ≃ 0.64. Different expressions exist for estimating the maximal capillary pressure Ψ cap before recession of menisci, and we will assume as in Ref. 48 for the case of fully wetted particles:

Ψ cap = - 2γ r p = - 3ϕ d 1 -ϕ d γ a , (35)
leading to Ψ cap ≃ -5.33γ/a for ϕ d ≃ 0.64. Finally, we chose to estimate the permeability of the randomlypacked assembly of colloids by the Carman-Koseny relation: 47

κ = a 2 (1 -ϕ d ) 3 45ϕ 2 d . (36)

Validity of Eq. (17)

Peppin et al. described the dynamics of the solidification of hardsphere colloidal suspensions in the experiment described in Fig. 1 taking into account the transport of the colloids by advection and diffusion. 34,35 In these models, Eq. ( 17) can be demonstrated strictly from the colloid mass conservation but for a constant evaporation rate J. This is due to the existence of a steady solution of the advection-diffusion model governing the colloid con-centration, in the reference frame of the growth of the packing. For a non-constant J, there is no such steady solution a priori, and Eq. ( 17) is only an approximation.

Nevertheless, Eq. ( 17) remains a good approximation as long as the temporal variation of the colloid concentration field for x < -l can be neglected. For quasi-static conditions and negligible buoyancy, 49 the colloid concentration upstream the growing front can be approximated by the exponential decay ϕ(x,t) ≃ ϕ 0 + (ϕ dϕ 0 ) exp[(x + l)/l g ], with l g = D/J and D the collective diffusion coefficient of the dispersion, well approximated by the Stokes-Einstein relation for hard-spheres. 34 The colloid mass balance, now taking into account the temporal variation of this colloid concentration profile leads to:

(ϕ d -ϕ 0 ) dl dt + dl g dt = ϕ 0 J . (37)
The approximate colloid mass balance Eq. ( 17) comes from Eq. ( 37) when dl g dt ≪ dl dt , equivalently when:

D ≪ - κΨ ⋆ η w , (38)
using l g = D/J, Eqs. ( 11) and ( 16), with either Ψ ⋆ = Ψ ext or Ψ ⋆ = Ψ cap depending on the regime considered. Relations (35) and (36) then show that Eq. ( 38) is always fulfilled for hardspheres over a wide range of radii a, and that Eq. ( 17) can be safely used to describe the growth rate. Moreover, these numerical applications show that the term κΨ ⋆ η w in Eq. ( 38) is very large. It is thus very likely that Eq. ( 17) remains valid even for colloids with long-range repulsive interactions, such as those often studied in the works cited previously, 17,18,20,22,37 and for which colinteractions are known to enhance the diffusion coefficient D. 27,50

Validity of Γ ≫ 1

The capillarity-limited growth described in particular by (16) and (32) follows from Γ ≫ 1 with Γ given by Eq. The calculation of Γ requires nevertheless to estimate the effective coefficient of the vapor in the porous packing D p w .

The regime of vapor diffusion within the pores of the colloidal packing depends on the value of the Knudsen number defined as Kn = λ /r p , with λ the mean free path in air (λ ≃ 64-68 nm at room conditions 51 in the range RH = [0-1]). Using Eq. ( 34), the transition Kn = 1 corresponds to radii a ≃ 170-180 nm. For Kn ≪ 1, there are different theoretical estimates of the effective gas diffusion coefficient in a porous medium depending on its 3D morphology. 52 We chose the following relation:

D p,n w ≃ 2(1 -ϕ d ) 2 + ϕ d D air w , (39)
derived by Neale and Nader, 53 that gives good estimates for homogeneous and isotropic assemblies of spheres, D air w being the diffusion coefficient in air. In the Knudsen regime, Kn ≫ 1, Huizenga and Smith 46 showed using measurements on randomly closepackings of monodisperse silica colloids, that the effective diffu- Fig. 6 Γ defined by Eq. ( 15) versus colloid radius a, estimated with the relations for hard-spheres given in Sec. 3.3. Inset: effective diffusion coefficient D p w of the vapor in the porous medium versus a, see Eq. ( 41). The dotted line at low a corresponds to Knuden diffusion Eq. ( 40), and the one at large a to effective diffusion given by Eq. (39).

sion coefficient follows:

D p,k w ≃ 1 -ϕ d τ 4r p 3 2RT πM w , (40)
with the tortuosity τ ≃ 1.47 ± 0.1 and M w the molar mass of water.

To cover a wide range of colloid radius a, and thus of Knudsen number, we finally estimate the effective diffusion coefficient D p w by:

1 D p w = 1 D p,n w + 1 D p,k w , (41)
as done for instance in Ref. 54 , see the inset of Fig. 6. Figure 6 shows Γ estimated using Eqs. (35), (36), and (41) versus the colloid radius a for T = 22 • C. These data show that Γ ≫ 1 regardless of the size of the colloids. This numerical application validates the assumption done in Sec. 2.2.2 of minute recession of the menisci within the dense packing, and the existence of a capillarity-limited regime. Interestingly, Fig. 6 also shows that this assumption is a priori not due to Knudsen diffusion, at least for large colloids (a > 100 nm), a hypothesis often mentioned to explain the slowdown observed in confined drying experiments. 18,36,37

Numerical application of the model

We now illustrate the results of our model on a relevant experimental configuration: a microfluidic capillary of diameter h = 25 µm and centimetric length. For diffusion-limited transport from a flat air-water interface flush just with the outlet of the capillary, the mass transfer coefficient in Eq. ( 4) is: 55 k ≃ 8D air w /(πh), leading to J 0 ≃ 50 µm/s at T = 22 • C. Similar values of evaporation rates could be obtained with cells of larger cross-sectional dimensions, but with convection enhancing evaporation in the ambient air.

We use the relations given in Sec. 3.3 for hard-spheres to compute the crossover length l c between the linear regime l ∝ t at early the root growth l ∝ √ In regime, l l kel l kel given by Eq. ( 29), while l c = l rec with l rec given by Eq. ( 25) in the capillarity-limited regime. Figure 7 reports these data versus the colloid size a, for various RH (more precisely the minimal length between l kel and l rec ). Figure 7 also shows for each RH, the critical radius a below which menisci are never expected to recede, see the criterion given by Eq. ( 24). This figure shows that most of the regimes predicted by the model are a priori expected to occur for colloids with radii ranging from ≃ 5 to ≃ 100 nm in a cell of centimetric length for the conditions of evaporation considered here. Interestingly, these data also reveal the important role played by RH on the transition between flow-limited and the capillarity-limited regimes, but also on the l c values particularly in the latter regime. Some experiments 17,18,20,22,37 have measured l c , in particular for silica colloids of small size a ≤ 20 nm (Ludox dispersions). Unfortunately, the relative humidity was not a control parameter in these experiments, and was not often even reported. Dufresne et al. have also reported non-linear variations of l c with a in the range a = 6-26 nm. 17 The scaling laws shown in Fig. 7, l c ∼ a 2 in the flow-limited regime and l ∼ a in the capillarity-limited regime, could maybe explain these results. As shown by our model, only the control of the relative humidity in these experiments would help to test this assumption.

Conclusions

In this work, we have developed a model to describe the dynamics of the consolidated front in confined directional drying experiments of colloidal dispersions. This model predicts at long times, a square root slowdown, l ∝ √ t, with two different origins: a flowlimited regime linked to the Kelvin effect and a capillarity-limited regime due to the recession of the evaporation interface. Our model also clarifies various assumptions commonly used in this context, and shows that Knudsen diffusion is not strictly at the origin of the capillarity-limited regime, as the latter can a priori also be observed for large colloid sizes, see Fig. 6.

Beyond these results, our model highlights the importance of controlling the evaporation conditions in these experiments, and in particular the relative humidity RH. A control of the relative humidity would indeed allow to discriminate between the regime limited by capillarity, whose growth rate does not strictly depend on RH, and the regime limited by the flow showing a dependence with RH through Ψ ext , see Eq. (33). Interestingly, the evaporation rate through the colloidal packing in the capillarity-limited regime is independent of RH, see Eq. ( 16). This was also demonstrated by Vincent et al. with experiments on drying-induced flows through a nano-porous silicon layer. [START_REF] Vincent | [END_REF] This humidity-insensitive water evaporation is reminiscent of recent results for molecular mixtures, 56,57 but originates here in subtle couplings between evaporation, capillarity and vapor diffusion, as for example also in Ref. 58 Importantly, no experiment have reported, to our knowledge, even the minimal recession of the drying interface within the colloidal packing in the regime of square root growth l ∝ √ t, observed particularly with colloids of small radii, a ≤ 20 nm. 17,18,20,22,37 This result could suggest that the flow-limited regime is at the origin of these observations, but would deserve more precise local experimental measurements to confirm it. Furthermore, the capillarity-limited regime shares many similarities with the drying of a porous medium, for which the Kelvin effect can also play an important role. 45,59 Local measurements of the water content in the colloidal packing (saturation) as done for instance in Ref. 45 using magnetic resonance imaging would also be a major asset for confined directional drying experiments.

Also, the high values of water potential Ψ ⋆ due to the flow induced by evaporation are associated with tensile stresses that are often released in the experiments by the appearance of fractures, [16][17][18][19][20] or the delamination of the packing from the cell walls. 22,23 The model discussed in this article does not take these different phenomena into account at all, and it seems crucial to us to evaluate their role on the transport of water by evaporation and on the different regimes predicted by the model.

Finally, the confined directional drying configuration described in this work shares strong similarities with the important problem of colloid accumulation and clogging during drying of a porous medium [60][61][62][63] , for which Fig. 1 could locally describe the porescale dynamics. Here again, the main ingredients of our model and some of its results could a priori also apply to these cases, and in part explain the global drying dynamics of colloidal dispersions in porous media.

Fig. 1

 1 Fig. 1 Confined directional drying of a colloidal dispersion. (a) Capillary forces hold the air-water menisci (radius of curvature r) at the outlet of the cell. (b) Recession of the evaporation interface within the colloidal packing, leading to a dry layer of thickness ξ (radius of curvature r c ). The water potential profiles are shown below for both scenarios.

Fig. 2

 2 Fig.2Evaporation rate given by Eqs. (4) and (8) versus 1 -RH, for l/l k = 10 -2 , 10 -1 , 1, 10, 10 2 , and 10 3 . Inset: zoom for the high values of l/l k . The dotted lines are the approximations Eq. (11) for l ≫ l k .

n o r e c e s s i o n 10

 10

Fig. 4

 4 Fig.4Yellow: numerical solution of Eq. (20) RH = 0.5. The dashed lines show the approximation Eq. (27) at low l, and the flow-limited regime Eq. (28) at larger l. The other curves are given by Eq. (30) for α = 2 × 10 -1 , 2 × 10 -2 , 2 × 10 -3 , and 2 × 10 -4 (from orange to brown).

Fig. 7

 7 Fig.7 Critical length l c for the crossover between the linear growth l ∝ t and the square root growth l ∝ √ t for the experimental conditions given in Sec. 3.3.3, and various relative humidity, RH = 0.2 (blue, ■), 0.5 (red, ⋄), and 0.8 (dark, •). The symbols show the critical radii below which the menisci do not recede, Eq.

1-10 J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,

[START_REF] Keddie | Fundamentals of latex film formation. Processes and Properties[END_REF][START_REF] Li | [END_REF][3][4][START_REF] Goehring | Desiccation Cracks and their Patterns: Formation and Modelling in Science and Nature[END_REF][START_REF] Bacchin | [END_REF][7][8][9][10]

| 1-10 J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,

J o u r n a l N a me , [ y e a r ] , [ v o l . ] ,

Acknowledgements

We acknowledge Solvay, CNRS and the ANR program Grant No. ANR-18-CE06-0021 for financial support.

Notes and references