J Garcia-Alfaro
email: joaquin.garcia-alfaro@acm.org

M Barbeau
email: barbeau@scs.carleton.ca

E Kranakis
email: kranakis@scs.carleton.ca

Proactive Threshold Cryptosystem for EPC Tags ⋆

Keywords: Radio Frequency Identification (RFID), Electronic Product Code (EPC), Wireless Security, Threshold Cryptography, Privacy

We define three privacy-preserving solutions to the problem of distributing secrets between manufacturers and vendors of items labeled with Electronic Product Code (EPC) Gen2 tags. The solutions rely on the use of an anonymous threshold secret sharing scheme that allows the exchange of blinded information between readers and tags. Moreover, our secret sharing scheme allows self-renewal of shares with secret preservation between asynchronous shareholders. The first two solutions address the eavesdropping and rogue scanning threats. The third solution mitigates as well tracking threats.

INTRODUCTION

The EPCglobal class-1 generation-2 (Gen2 for short) specification [START_REF] Epcglobal | EPC Radio-frequency Identity Protocols Class-1 Generation-2 UHF[END_REF], approved as ISO18000-6C, has been reported vulnerable to privacy attacks in previous studies [START_REF] Juels | RFID Security and Privacy: A Research Survey[END_REF][START_REF] Garcia-Alfaro | Analysis of Threats to the Security of EPC Networks[END_REF]. Consumer privacy is indeed an important concern about Gen2 applications. For instance, the use of Gen2 tags for item-level passive tagging [START_REF]EPC Item Level Tagging Joint Requirements Group[END_REF] of end-user goods, allows customers to enjoy the benefits of RFID technology, but anyone with a compatible Gen2 reader can access a consumer's purchase data. The readability of tag identification in the Gen2 protocol clearly violates consumer privacy. The issue must be properly handled before releasing this technology for item-level tagging. A radical solution is the use of the kill feature that disables Gen2 tags at purchase time [START_REF] Epcglobal | EPC Radio-frequency Identity Protocols Class-1 Generation-2 UHF[END_REF]. This solution is far from being effective because it requires spending more time at checkout stands and voids the benefits of the RFID technology offered to customers, such as processing of returns and automated recycling. Our goal is to provide non destructive lightweight alternatives that preserve consumer privacy. In this paper, we survey related works and present an original scheme for the construction of a lightweight threshold cryptosystem [START_REF] Karnin | On Secret Sharing Systems[END_REF] that can be deployed on low-cost Gen2 systems. The scheme protects EPC Gen2 tag data against eavesdropping, rogue scanning, and tracking by malicious readers.

Tag Identification (TID) disclosure on the Gen2 Protocol

The memory of an EPC Gen2 tag is separated into four independent blocks: reserved memory, EPC data, Tag Identification (TID), and User memory. Gen2 tags communicate this information by accumulating power from reader interrogations [START_REF] Epcglobal | EPC Radio-frequency Identity Protocols Class-1 Generation-2 UHF[END_REF]. Figure 1 shows the steps of the EPC Gen2 protocol for product inventory. In Step 1, a reader queries the tag and selects one of the following options: select, inventory, or access [START_REF] Epcglobal | EPC Radio-frequency Identity Protocols Class-1 Generation-2 UHF[END_REF]. Figure 1 represents the execution of an inventory query. It assumes that a select operation has been previously completed in order to singulate a specific tag from the population of tags. When the tag receives the inventory query, it returns a 16-bit random string denoted as RN16. This random string is temporarily stored in the tag memory. As a response to the inventory query, the tag enters in the ready state, and backscatters in Step 2 the random string RN16. In Step 3, the reader replies to the tag a copy of the random string, as an acknowledgment. If the echoed string matches the copy of the RN16 squence stored in the tag memory, the tag enters in the acknowledged state and returns its corresponding tag identification (TID).

Security features on Gen2 tags are minimal. They protect message integrity via 16-bit Cyclic Redundancy Codes (CRC) and generate 16-bit pseudo random strings. Let us observe that any compatible Gen2 reader can access the TID. This is due to the lack of authentication between a Gen2 reader and a Gen2 tag. To overcome this problem, a number of solutions have been proposed in the literature. Two solutions proposed rely on the use of cryptographic primitives to encrypt TIDs and the use of pseudonyms for the TIDs. Both solutions require that the reader and tag share a common secret (either a key to decrypt a protected TID or a property to map a pseudonym to the true TID). Therefore, an effective mechanism for the distribution of secrets among the entities, readers and tags, of a supply chain must be introduced.

We present in this paper a threshold cryptosystem that provides both consumer privacy and distribution of secrets. Our solution addresses the following three threats: [START_REF] Garcia-Alfaro | A Proactive Threshold Secret Sharing Scheme Handling Gen2 Privacy Threats[END_REF] Eavesdropping: an adversary listening passively to the RF (Radio Frequency) communications between a tag and reader aiming to intercept the TID; (2) Rogue Scanning: an adversary interacting actively with a tag to access the TID; [START_REF] Juels | RFID Security and Privacy: A Research Survey[END_REF] Tracking: an adversary correlating RF communications to either passively or actively identify an instance of a given tag. We present three different variants of our solution for the exchange of secrets between manufacturers and vendors of Gen2 labeled items. The two first variants handle the eavesdropping and rogue scanning threats. The proactiveness of the third variant addresses, in addition, the tracking threat.

The main properties of our approach are: (1) low-cost Gen2 tag renewal with secret preservation and without the need to synchronize to a reader performing an inventory process or any other tags holding shares of the same secret; (2) size of shares compact enough to fit into the inventory responses of low-cost EPC Gen2 tags (i.e., less than 528 bits, as suggested by EPCglobal in [START_REF] Epcglobal | EPC Radio-frequency Identity Protocols Class-1 Generation-2 UHF[END_REF]); (3) secret sharing construction that guarantees strong security; (4) reconstruction of the secret does not require the identity of the shareholders, e.g., the TIDs. The remainder of the paper is organized as follows. Section 2 surveys privacy-preserving solutions for low-cost RFID systems. Section 3 presents the formalization of our proposal. Section 4 provides some simulation and experimental details. Section 5 concludes the paper.

RELATED WORK

The high constraints of EPC Gen2 tags makes challenging the use of standard cryptography-based solutions for the design of privacy preserving mechanisms. In this section, we survey solution and trends recently published in the literature.

Use of Standard Cryptography-based Solutions

MAC (Message Authentication Code) based security protocols are among the first solutions discussed in the literature for securing low-cost RFID applications. In [START_REF] Takaragi | An Ultra Small Individual Recognition Security Chip[END_REF], for example, Takaragi et al. present a solution based on CMOS technology that requires less than four thousand gates to generate MACs using 128 bit identifiers stored permanently in tags at manufacturing time. Each identifier relies on an initial authentication code concatenated with chip manufacturer data. The result of this concatenation is posteriorly hashed with a given secret to derive a final MAC. This MAC is communicated from manufacturers to clients and shared by readers and tags. The main benefit of this approach is that it increases the technical difficulties of performing eavesdropping and rogue scanning. However, the use of static identifiers embedded in tags at manufacturing time does not solve the tracking threat. Moreover, brute force attacks can eventually reveal the secrets shared between readers and tags. The discovery of secrets could lead to counterfeit tags.

An enhanced solution relies on the use of lock-based access control. In [START_REF] Weis | Security and Privacy Aspects of Low-cost Radio Frequency Identification Systems[END_REF], Weis et al. propose a mechanism to prevent unauthorized readers from reading tag contents. A secret is communicated by authorized readers to tags on a secure channel. Every tag, using an internal function, performs a hash of this secret and stores the result in its internal memory. Then, the tag enters into a locked state in which it responds to any query with the stored hash value. Weis et al. also describe a mechanism for unlocking tags, if such an action is needed by authorized readers (i.e., to temporarily enable reading of private data). Regarding the tracking threat, Ohkubo et al. propose in [START_REF] Ohkubo | Efficient Hash-chain Based RFID Privacy Protection Scheme[END_REF] the use of hash chains achieving on-tag evolution of identifiers. Avoine and Oechslin discuss in [START_REF] Avoine | A Scalable and Provably Secure Hash Based RFID Protocol[END_REF] limitations of the aforementioned approach. They propose an enhanced hash-based RFID protocol to address eavesdropping, rogue scanning, and tracking by using timestamps. Similarly, Henrici and Müller discuss in [START_REF] Henrici | Hash-based Enhancement of Location Privacy for Radio-frequency Identification Devices Using Varying Identifiers[END_REF] some weaknesses in the lock-based schemes presented in [START_REF] Weis | Security and Privacy Aspects of Low-cost Radio Frequency Identification Systems[END_REF][START_REF] Ohkubo | Efficient Hash-chain Based RFID Privacy Protection Scheme[END_REF] and present an improved mechanism. Several other improvements and lock-based protocols, most of them inspired by lightweight cryptography for devices such as smart cards, can be found in [START_REF] Molnar | Privacy and Security in Library RFID: Issues, Practices, and Architectures[END_REF][START_REF] Dimitriou | A Lightweight RFID Protocol to Protect Against Traceability Cloning Attacks[END_REF][START_REF] Peris-Lopez | An Efficient Authentication Protocol for RFID Systems Resistant to Active Attacks[END_REF].

Hardware Limitations

Note that the approaches reviewed in Subsection 2.1 require the implementation of efficient one-way hash primitives within low-cost RFID tags. It is the main challenge of these proposals. Resource requirements of standard oneway hash functions, such as MD4, MD5, and SHA-128/SHA-256, exceed the constraints of low-cost Gen2 tags [START_REF] Epcglobal | EPC Radio-frequency Identity Protocols Class-1 Generation-2 UHF[END_REF]. The implementation of these functions require from seven thousand to over ten thousand logic gates; and from six hundred to over one thousand two hundred clock cycles [START_REF] Peris-Lopez | An Efficient Authentication Protocol for RFID Systems Resistant to Active Attacks[END_REF]. The complexity of these standard one-way hash functions is therefore an obstacle for their deployment on Gen2 tags.

The use of standard encryption engines for the construction of hash operations has also been discussed in the literature. For example, the use of Elliptic Curve Cryptosystem (ECC) [START_REF] Menezes | Handbook of Applied Cryptography[END_REF] for the implementation of one-way hash primitives on RFID tags has been studied in [START_REF] Wolkerstorfer | Is Elliptic-Curve Cryptography Suitable to Secure RFID Tags? In Ecrypt Workshop on RFID and Lightweight Crypto[END_REF]. Its use of small key sizes is seen as very promising for providing an adequate level of computational security at a relatively low cost [START_REF]Networked RFID Systems and Lightweight Cryptography -Raising Barriers to Product Counterfeiting[END_REF]. An ECC implementation for low-cost RFID tags can be found in [START_REF] Batina | An Elliptic Curve Processor Suitable for RFID-tags[END_REF]. In [START_REF] Feldhofer | Strong Authentication for RFID Systems Using the AES Algorithm[END_REF], on the other hand, Feldhofer et al. present a 128-bit implementation of the Advanced Encryption Standard (AES) [START_REF] Daemen | The Design of Rijndael: AES-the Advanced Encryption Standard[END_REF] on an IC of about three thousand five hundred gates with a power consumption of less than nine microampers at a frequency of 100 kHz. Although this implementation is considerably simpler than any previous implementation of the AES algorithm, its requirements are still seen as too high for low-cost RFID tags.

Towards Secret Sharing Strategies

The use of secret sharing schemes [START_REF] Karnin | On Secret Sharing Systems[END_REF] is proposed by Langheinrich and Martin in [START_REF] Langheinrich | RFID Privacy Using Spatially Distributed Shared Secrets[END_REF][START_REF] Langheinrich | Practical Minimalist Cryptography for RFID Privacy[END_REF] as a key solution for addressing authentication in Gen2 scenarios (e.g., supply chain applications of the retail industry). The work presented in [START_REF] Langheinrich | RFID Privacy Using Spatially Distributed Shared Secrets[END_REF] simplifies the lookup process performed on back-end databases for identifying tags, while guaranteeing authentication. TIDs, seen in this work as the secrets that must be shared between readers and tags, are encoded as a set of shares, and stored in the internal memory of tags. The authors propose the use of a Perfect Secret Sharing (PSS) scheme, in which the size of the shares is equivalent to the size of the secret, based on the (t,n)-threshold secret sharing scheme introduced in [START_REF] Shamir | How to Share a Secret[END_REF]. The combination of shares at the reader side leads to the reconstruction of original TIDs. To prevent brute-force scanning from unauthorized readers -trying to obtain the complete set of shares -the authors propose a time-limited access that controls the amount of data sent from tags to readers. At the same time, a cache based process ensures that authorized readers can quickly identify tags. In [START_REF] Langheinrich | RFID Privacy Using Spatially Distributed Shared Secrets[END_REF], the authors extend the previous proposal to spread the set of shares across multiple tags. Still based on Shamir's perfect secret sharing scheme, this new approach aims at encoding the identifier of an item tagged with multiple RFID devices by distributing it as multiple shares stored within tags. Authentication is achieved by requiring readers to obtain and combine the set of shares.

A different use of secret sharing schemes is presented by Juels, Pappu, and Parno in [START_REF] Juels | Unidirectional Key Distribution Across Time and Space with Applications to RFID Security[END_REF]. They propose the use of a dispersion of secrets strategy, rather than the aggregation strategy used by Langheinrich and Marti. In this approach, a secret used to encrypt Gen2 TIDs is split into multiple shares and distributed among multiple tagged items. Construction and recombination of shares is based on a Ramp Secret Sharing (RSS) scheme, in which the size of each share is considerably smaller than the size of the secret, at the price of leaking out secret information for unqualified sets of shares. To identify the tags, a reader must collect a number of shares above a threshold. At the manufacturer facility, large quantities of items of the same product, initially tagged together with shares of the same secret, guarantee that authorized readers can always reconstruct the secret and, therefore, decrypt the TID of the tagged items. At the consumer side, the items get isolated. Without the space proximity of other items holding the remainder shares of the secret, an unauthorized reader cannot obtain the sufficient number of shares to reconstruct the key that allows identifying the TID. Privacy is, therefore, achieved though the dispersion of secrets and encrypted identifiers. Moreover, the proposed scheme helps to improve the authentication process of tags. Assuming that t shares are necessary for readers to obtain the EPC data assigned to a pallet, a situation where the number of shares obtained by readers is below t leads to conclude that unauthorized tags are present in the pallet.

The main limitation of this approach is that a critical privacy threat to consumers, i.e., the tracking threat defined in Section 1.1, is not addressed. It is a requirement stated by most authors, such as Juels and Weis in [START_REF] Juels | Defining Strong Privacy for RFID[END_REF]. Privacypreserving solutions for RFID applications must guarantee both anonymity and untraceability. In the sequel, we show that it is possible to improve privacy-preserving using secret sharing strategies. We have developed a proofof-concept threshold cryptosystem that provides, in addition to eavesdropping and rogue scanning, tracking protection.

ON THE CONSTRUCTION OF A PROACTIVE THRESHOLD SE-CRET SHARING SCHEME FOR EPC GEN2

The construction of our proactive (t,n)-threshold cryptosystem relies on the computation of the Moore-Penrose pseudoinverse of a homogeneous system of n linear equations with t unknowns (where t < n) over a finite field Z p ,

a 11 x 1 + a 12 x 2 + a 13 x 3 + • • • + a 1t x t = 0 (mod p), a 21 x 1 + a 22 x 2 + a 23 x 3 + • • • + a 2t x t = 0 (mod p), . . . a n1 x 1 + a n2 x 2 + a n3 x 3 + • • • + a nt x t = 0 (mod p),
in which t and n are positive integers, and p is a prime number. The vector columns of the coefficient matrix A associated to the system of linear equations are linearly independent, i.e., matrix A has rank t and so the vector columns of A span an inner-product subspace in Z n×t p of dimension t. The Moore-Penrose pseudoinverse (also called the generalized inverse) of a non-square matrix A ∈ Z n×t p , hereinafter denoted as A † , is the closest representation that A can get to its inverse (since non-square matrices, i.e., n = t, do not have an inverse). Let us notice that if rank(A) = t = n, i.e., A is a full rank square matrix, the Moore-Penrose pseudoinverse of A is certainly equivalent to the inverse matrix A -1 , i.e.,

A † = A -1 | A ∈ Z n×t p ∧ rank(A) = t = n (1)
Otherwise, the Moore-Penrose pseudoinverse of a rectangular matrix A exists if and only if the subspaces Ker A (null space of matrix A) and Im A (range space of matrix A) have trivial intersection with their orthogonals. In the case that A ∈ Z n×t p has rank(A)= t, it can be proved that A † exists and it can be computed as follows:

A † = (A ⊥ A) -1 A ⊥ ∈ Z t×n p | A ∈ Z n×t p ∧ rank(A) = t = n, (2)
in which A ⊥ denotes the transpose of matrix A. It can also be proved, cf. [START_REF] Meyer | Matrix Analysis and Applied Linear Algebra[END_REF],

that if A ∈ Z n×t p | rank(A) = t, A †
is the unique solution that satisfies all of the following four equations defined by Penrose in [START_REF] Penrose | A Generalized Inverse for Matrices[END_REF]:

(A A † ) ⊥ = A A † , A † A A † = A † , (A † A) ⊥ = A † A, and A A † A = A (3)
For our specific construction, we are interested in showing that the resulting matrix A † keeps the orthogonal projection property required in [START_REF] Penrose | A Generalized Inverse for Matrices[END_REF]. Indeed, we are interested in showing that the resulting matrix P A computed as

P A = A A † ∈ Z n×n p | A ∈ Z n×t p ∧ rank(A) = t = n (4)
is indeed an orthogonal projector that satisfies the idempotent property (meaning that

P k A = P A for all k ≥ 2. Certainly, if P A = A A † , then P 2 A = (A A † ) (A A † ), i.e., P 2 A = (A A † A) A † ).
From Equation (3), we obtain that P 2 A = A A † , i.e., P 2 A = P A , and so P k A = P A for all k ≥ 2. Therefore, if A ∈ Z n×t p and rank(A)= t, then A A † ∈ Z n×n p is an orthogonal projector. Figure 2 shows how the orthogonal projector P A can be used to project a vector v onto the column space of matrix A.

The Moore-Penrose pseudoinverse is a very useful technique used in many engineering fields such as error correction, identification, control design, and structural dynamics. For an over-determined system of linear equations without solution, the Moore-Penrose pseudoinverse finds the least squares solution (i.e., projection of the solution onto the range space of the coefficient matrix of the system). It is also helpful to find the infinite set of solutions in the range space of under-determined set of equations (i.e., fewer constraints than unknowns). The computation of the Moore-Penrose pseudoinverse of a homogenous system of t linear equations with n unknowns (e.g., the computation of the pseudoinverse of matrix A ⊥ ∈ Z t×n p ) is hence a valid alternative for the construction of our proactive threshold secret sharing.

p = PA " v q = (I -PA) " v v A A " p = PA q = (I -PA) " v v " v A " FIGURE 2
Orthogonal Projection of a Vector v onto the Subspace Spanned by the Column Vectors of Matrix A.

Basic (t,n)-Threshold Secret Sharing Scheme Based on the Invariance Property of Orthogonal Projectors

Orthogonal projectors have already been used for the construction of threshold secret sharing schemes. In [START_REF] Bai | A Reliable (k,n)-Image Secret Sharing Scheme[END_REF][START_REF] Bai | A Strong Ramp Secret Sharing Scheme Using Matrix Projection[END_REF], for example, the invariance property of orthogonal projectors is used for the redundant storage of computer images. Indeed, an asynchronous proactive (t,n)-threshold secret sharing scheme can be constructed based on the same observation -meaning that the invariance property of orthogonal projectors can be used to allow shareholders to renew their share without synchronization with other parties and without altering the secret. The key idea of the proposed approach is that the orthogonal projector P A computed from Equation ( 4) and a random matrix A ∈ Z n×t p with rank t is always equivalent to the projector P B obtained from the same equation and any t independent random range images spanned from A.

Before going any further, let us start with a simple example that depicts the basic idea of our approach. It exemplifies the construction of a (2, 3)threshold, non-proactive yet, cryptosystem; and the reconstruction process by three independent reconstruction processes. Given two matrices such that A is a random matrix composed of two linearly independent column vectors a 1 , a 2 ∈ Z 3×1 31 , i.e., rank(A) is equal to 2; and X is a random matrix composed of three linearly independent column vectors x 1 , x 2 , x 3 ∈ Z 2×1 31 , i.e., rank(X) is equal to 3. Note that we simplify the notation, assuming A = [a 1 , a 2 , . . . , a t ], where each a i is the i-th column vector of matrix A; and X = [x 1 , x 2 , . . . , x n ] where each x i is the i-th column vector of matrix X. Let   be the resulting matrix obtained by multiplying matrices A and X. We assume hereafter that the column vectors a ′ 1 , a ′ 2 , and a ′ 3 in matrix A ′ are indeed the shares of our cryptosystem; and that P A ∈ Z 3×3 31 is the secret of the cryptosystem, in which P A is the orthogonal projector obtained by applying Equation (4) to matrix A.

A ∈ Z 3×2 31 , X ∈ Z 2×3 31 , A =   7
A ′ ∈ Z 3×3 31 =   19
Let us now assume that a distribution process δ disseminates the shares a ′ 1 , a ′ 2 , a ′ 3 ∈ A ′ to three independent shareholders α, β , and γ. We define the following three column vectors:

V α =   19 20 16   , V β =   15 28 16   , V γ =   27 2 24  
as the corresponding shares held respectively by α, β , and γ. We also assume that a reconstruction process ρ 1 requests to shareholders α and β their respective shares (notice that our example describes a (2, 3)-threshold cryptosystem and so only two shares suffice to reconstruct the secret). A second reconstruction process ρ 2 requests to shareholders α and γ their respective shares. Finally, a third reconstruction process ρ 3 requests to shareholders γ and β their shares. Processes ρ 1 , ρ 2 , and ρ 3 build, independently, three re- construction matrices B 1 , B 2 , and B 3 (by simply joining the share vectors they collected from each shareholder): Therefore, the three processes ρ 1 , ρ 2 , ρ 3 may successfully reconstruct the secret (i.e., P A ) by performing the same operation described by Equation (4). The following theorem establishes the corretness of the approach for the general case.

B 1 =   19
Theorem 1 Let A ∈ Z n×t p be a random matrix of rank t. Let A ′ ∈ Z n×n p be the result of multiplying matrix A with a set of n linearly independent column vectors x 1 , x 2 , . . . ,

x n ∈ Z t×1 p , i.e., A ′ = Ax i (mod p) ∀x i ∈ [x 1 , x 2 , . . . , x n ].
Let B be any submatrix from A ′ with exactly t column vectors. Then, the orthogonal projectors P A and P B derived from Equation ( 4) are identical.

Proof Note that P A = A A † and P B = B B † are the orthogonal projectors obtained by applying Equation (4) to both A and B. Since B is any submatrix derived from A ′ with exactly t column vectors, we can also denote B as the resulting matrix obtained by multiplying A ∈ Z n×t p times a given matrix X ∈ Z t×t p . Therefore, P B = B B † is equal to P B = (A X)(A X) † and so to P B = A X X † A † . We know from Equation (1) that X † = X -1 when X is a square matrix. Therefore, P B = A X X -1 A † . Since matrix X gets cancelled, we obtain that P B = A A † and so identical to P A .

Efficiency

The efficiency of a secret sharing scheme can be evaluated in terms of the information entropy of its shares and secret of the cryptosystem [START_REF] Iwamoto | Strongly Secure Ramp Secret Sharing Schemes For General Access Structures In Information Processing Letters[END_REF]. A secret sharing scheme is said to be perfect if it holds that the entropy of the shares is greater that or equal to the entropy of the secret. As a consequence, the size of each share of a perfect secret sharing scheme must be equal or greater than the size of the secret. This is an inconvenient to the hardware limitations of the RFID model discussed in Section 2.2. RSS may considerably improve this efficiency, by allowing a trade-off between security and size of the shares [START_REF] Blakley | Security of Ramp Schemes[END_REF]. This is the case of the approach presented in the previous section (cf. Section 3.1). Notice that the size of each share a ′ i ∈ A ′ of our construction is considerably smaller than the size of the secret P A . More precisely, every share a ′ i is a column vector in Z n×1 p , while the size of the secret is a matrix in Z n×n p , i.e., the size of every share is n times smaller than the secret.

To analyze the robustness of a RSS scheme, in terms of its security, it is necessary to quantify the amount of information about the secret that an intermediate set of shares, smaller than the threshold t, may leak out. This leakage of secret information represents the size of the ramp, in which a small ramp provides stronger security to the scheme than a larger ramp. Yakamoto proposed in [START_REF] Yamamoto | On Secret Sharing Systems Using (k,l,n)-threshold Scheme[END_REF] to quantify the exposure of secret information from each share by defining a second threshold t ′ , where 0 < t ′ ≤ t. By definition, a qualified coalition of t shares may reconstruct the secret. An unqualified coalition of tt ′ shares cannot reconstruct the secret, but leaks out information about it. Less than t ′ shares may not reconstruct the secret and does not reveal any information about the secret. The amount of information leaked out from the secret by an unqualified coalition of tt ′ shares can be quantified in terms of information entropy. Yakamoto proved in [START_REF] Yamamoto | On Secret Sharing Systems Using (k,l,n)-threshold Scheme[END_REF] that the security of a ramp secret sharing scheme is strong enough when the following equivalence applies:

H(S|C) = t -t ′ t H(S), (5)
in which H(S) is the information entropy of the secret, and C is an unqualified coaliation of tt ′ shares. We prove in the sequel that the security of the threshold cryptosystem presented in Section 3.1 is, according to [START_REF] Yamamoto | On Secret Sharing Systems Using (k,l,n)-threshold Scheme[END_REF], strong enough.

Theorem 2 Let A ∈ Z n×t p be a random matrix of rank t. Let P A ∈ Z n×t p be the orthogonal projector obtained by applying Equation ( 4) to matrix A. Let A ′ ∈ Z n×n p be the result of multiplying matrix A with a set of n linearly independent column vectors x 1 , x 2 , . . . , x n ∈ Z t×1 p . Then, the basic (t,t ′ , n)threshold secret sharing scheme constructed from the invariance property of the orthogonal projector P A , in which matrix P A is the secret of the cryptosystem, and the column vectors a ′ 1 , a ′ 2 , . . . , a ′ n ∈ A ′ are the shares of the cryptosystem, is equivalent to Equation [START_REF]EPC Item Level Tagging Joint Requirements Group[END_REF].

Proof Since the information provided by matrix A derives P A by simply applying Equation ( 4), we know that H(P A |A) = 0. Using some information entropy algebra manipulation, we can use this result to decompose H(P A ) as

H(P A ) = H(P A |A) + H(A) -H(A|P A ) = H(A) -H(A|P A ) (6)
Notice that matrix A is any full rank matrix chosen uniformly at random from the sample space in Z n×t p . It is proved in [START_REF] Migler | Weight And Rank Of Matrices Over Finite Fields[END_REF] that there are exactly ∏ t-1 i=0 (p np i ) random matrices of rank t in Z n×t p . Therefore, we can compute H(A) as follows:

H(A) = log 2 t-1 ∏ i=0 (p n -p i ) (7)
Knowing A and P A easily leads to H(A|P A ). From Equations ( 3) and ( 4), we have that P A times A is equivalent to A, meaning that A is an eigenvector matrix of P A . Hence, the decomposition of

P A into t eigenvectors [v 1 , v 2 , . . . , v t ] = V ∈ Z n×t p
provides information about A. More precisely, matrix A can be obtained from V by using a transformation matrix W ∈ Z t×t p . Since the sample space from which matrix W can be uniformly chosen is exactly of size ∏ t-1 i=0 (p tp i ), we have that H(A|P A ) can be obtained as follows:

H(A|P A ) = log 2 t-1 ∏ i=0 (p t -p i ) (8)
Using Equations ( 7) and ( 8) we can now compute H(P A ) = H(A)-H(A|P A ):

H(P A ) = log 2 t-1 ∏ i=0 (p n -p i ) -log 2 t-1 ∏ i=0 (p t -p i ) (9)
Let us now quantify, in terms of entropy, the information about P A provided by an unqualified coalition

A ′ of t ′ shares, s.t., A ′ = [a ′ 1 , a ′ 2 , . . . , a ′ t ′ ],
and where 0 < t ′ < t. Since matrix A ′ can be seen as a random matrix of rank t ′ chosen uniformly from the sample space ∏ t ′ -1 i=0 (p np i ), we have that H(A ′ ) can be denoted as follows:

H(A ′ ) = log 2 t ′ -1 ∏ i=0 (p n -p i ) (10)
Matrix A ′ is also an eigenvector matrix of P A . The decomposition of

P A into t eigenvectors [v 1 , v 2 , . . . , v t ] = V ∈ Z n×t
p provides information about A ′ . Indeed, matrix A ′ can be obtained from V by using a transformation matrix W ′ ∈ Z t×t ′ p . Since the sample space from which matrix W ′ can be uniformly chosen is exactly of size ∏ t ′ -1 i=0 (p tp i ), we have that H(A ′ |P A ) can be obtained as follows:

H(A ′ |P A ) = log 2 t ′ -1 ∏ i=0 (p t -p i ) (11)
We can quantify the amount of information about P A provided by A ′ , i.e., H(P A |A ′ ), using the results from Equations ( 9), [START_REF] Avoine | A Scalable and Provably Secure Hash Based RFID Protocol[END_REF], and (11):

H(P A |A ′ ) = H(P A ) -H(A ′ ) + H(A ′ |P A ) = log 2 t-1 ∏ i=0 (p n -p i ) -log 2 t-1 ∏ i=0 (p t -p i ) - log 2 t ′ -1 ∏ i=0 (p n -p i ) + log 2 t ′ -1 ∏ i=0 (p t -p i ) ( 12
)
When p is a large number, we can simplify the logarithmic expressions in Equations ( 9) and [START_REF] Molnar | Privacy and Security in Library RFID: Issues, Practices, and Architectures[END_REF] to derive H(P A ) and H(P A |A ′ ) as the following approximations:

H(P A ) ≈ t(n -t) log 2 p H(P A |A ′ ) ≈ (t -t ′ )(n -t) log 2 p
We observe that the information entropy of P A , knowing A ′ , is approximatively t-t ′ t times the information entropy of P A :

H(P A |A ′ ) ≈ t -t ′ t H(P A ), (13)
which, according to Equation ( 5) provided in [START_REF] Yamamoto | On Secret Sharing Systems Using (k,l,n)-threshold Scheme[END_REF], guarantees that the security of the ramp threshold secret sharing scheme is strong enough.

Let us conclude this section by determining a value of t, in terms of n, that guarantees that t -1 shares cannot reconstruct the secret. Given that the secret is the orthogonal projection P A derived from the computation of Equation ( 4) and matrix A, and observing again that the projection of A onto the subspace spanned by its range space remains in the same place, i.e., P A • A = A, it is therefore trivial to observe that the projection of any share onto the same subspace does not change either. This effect can be used by a malicious adversary in order to discover P A by solving n consecutive equations of (t -1) shares. Since, by definition, a (t, n)-threshold secret sharing scheme must prevent any coalition of less than t shares from reconstructing the secret, the parameter t of our construction shall be bounded in terms of n as follows:

(t -1)n < n(n + 1) 2 -1, t < 3 + n 2 ( 14
)
From Theorems 1 and 2, we conclude that if t < 3+n 2 , the scheme presented in Section 3.1 is a strong ramp threshold secret sharing scheme in which exactly t shares may reconstruct the secret, but t -1 or fewer shares cannot.

Pseudo-Proactive Threshold Secret Sharing Scheme Based on the Invariance Property of Orthogonal Projectors and Multiplicative Noise for the Renewal of Shares

We significantly improve in this section the results presented in Section 3.1 by showing that the introduction of multiplicative noise in the coefficients of matrix A ′ does not affect the reconstruction phase. By multiplicative noise we assume independent scalar multiplication of column vector shares a ′ i ∈ A ′ and scalar random numbers r 1 , . . . , r k for stretching these vectors. Indeed, we show that the introduction of multiplicative noise into the column vectors of any reconstruction matrix B i obtained from t column vectors in A ′ does not affect the results.

The following example shows the key idea of this new version. Assuming again a (2, 3)-threshold secret sharing scheme based on the orthogonal projectors of matrices A ∈ Z 3×2 31 , X ∈ Z 2×3 31 , and If we now generate three matrices B 1 , B 2 , and B 3 as combinations of vector columns from

A ′ = AX ∈ Z 3×3 31 : A =    7
A ′ = [a ′ 1 , a ′ 2 , a ′ 3 ]
and multiplicative noise, such as is the result of multiplying matrix A with n linearly independent column vectors x 1 , x 2 , . . . , x n ∈ Z t×1 p , i.e., A ′ = Ax i (mod p) ∀x i ∈ [x 1 , x 2 , . . . , x n ]; then, any submatrix B derived from exactly t column vectors in A ′ , but streched by multiplicative noise, can still be factorized as B = A X ′ , where X ′ ∈ Z t×t p is a square random matrix resulting from the set of t linearly independent column vectors in X, but stretched by a specific scaling random number r modulo p. We know from Equation ( 1) that (X ′ ) † = (X ′ ) -1 when X ′ is square. Therefore, X ′ gets cancelled during the reconstruction phase, i.e., P B = A X ′ (X ′ ) -1 A † , and we obtain that P B = P A = A A † .

B 1 ∈ Z 3×2 31 = [5 • a ′ 1 , 17 • a ′ 2 ] (mod 31), B 2 ∈ Z 3×2 31 = [7 • a ′ 1 , 13 • a ′ 3 ] (mod 31), and B 3 ∈ Z 3×2 31 = [9 • a ′ 3 , 22 • a ′ 2 ] (mod 31): B 1 =   2

Proactive Threshold Secret Sharing Scheme Based on the Invariance Property of Orthogonal Projectors and Both Multiplicative and Additive Noise for the Renewal of Shares

We have seen in the previous section that every share in the set of shares derived from matrix A ′ can be independently transformed by adding multiplicative noise, and so generating numerically different shares, but still guaranteeing the invariance property of orthogonal projectors to always reconstruct the initial secret (i.e., the orthogonal projector P A derived from matrix A). However, even if the new shares are numerically different, any malicious adversary can successfully observe that the shares are always linearly dependent, since the transformation process is simply stretching the initial share by some scaling random factor r.

We solve this problem by combining both multiplicative and additive noise in the transformation process. The only requirement is to provide to the process in charge of reconstructing the secret a reference used in the transformation process. We assume that this reference is the last column vector in matrix A ′ . We also assume that the generation process in charge of the construction of A ′ guarantees that the last column vector is an un-ordered collection of distinct elements. Then, shareholders are given access to this reference to renew their shares with a linear combination of this reference column. Note that this reference column must be also known a priori by the reconstruction process, but not by any malicious adversary that has access to the renewed shares. Let us illustrate with an example the key idea of this version. Assuming a (2, 3)threshold secret sharing scheme based on matrices A ∈ Z 3×2 31 , X ∈ Z 2×3 31 , and Every shareholder is given column a ′ 3 and either column a ′ 1 or column a ′ 2 . Let us assume two shareholders α and β in the system, each holding one of the following two share pairs V α and V β :

A ′ ∈ Z 3×3 31 = Ax i (mod p) ∀x i ∈ X: A =    7
V α =   19 27 20 2 16 24   , V β =   15 27 28 2 16 24  
Let us assume that a reconstruction process ρ 1 requests to each share- holder their share combination. Both α and β return to ρ 1 a linear transfor- mation from the column vectors in their share pairs. Shareholder α generates a random value r α = 15, transforms v α1 into v α1 • 15 (mod31), and returns

b α ∈ Z 3×1 31 = v α1 + v α2 . Similarly, β generates a random value r β = 14, transforms v β 1 into v β 1 • 14 (mod31) and returns b β ∈ Z 3×1 31 = v β 1 + v β 2 .
Two other reconstruction processes ρ 2 and ρ 3 request to each share holder their shares. Shareholders α and β return to ρ 2 and ρ 3 two different linear com- binations from the column vectors in their share pairs. Shareholder

α returns b ′ α ∈ Z 3×1 31 = 28 • v α1 + v α2 to process ρ 2 , and b ′′ α ∈ Z 3×1 31 = 5 • v α1 + v α2 to process ρ 3 . Shareholder β returns b ′ β ∈ Z 3×1 31 = 19 • v β 1 + v β 2 to process ρ 2 , and b ′′ β ∈ Z 3×1 31 = 21 •V β 1 + V β 2 to
B i = r α • a ′ 1 + a ′ 3 , r β • a ′ 2 + a ′ = r α • Ax 1 + Ax 3 , r β • Ax 2 + Ax 3 = A (r α • x 1 + x 3 ) , A (r β • x 2 + x 3 ) = A x ′ 1 , A x ′ 2 ) = A X ′ i ( 15
)
in which r α and r β are the random factors introduced by each shareholder on every interrogation as multiplicative noise; and X ′ i ∈ Z 2×2 31 is a random full rank square matrix derived from A ′ , and so from A X, plus the multiplicative and additive noise introduced by the shareholders on every interrogation. Since matrix X ′ i is a square matrix, the equivalence defined in Equation (1) applies, i.e., X † i = X -1 i . Therefore, the computation of any orthogonal projector P Bi based on Equation (4) cancels matrix X ′ i and so P Bi is always identical to matrix P A . This establishes the general case of the new approach based on the proof of Theorem 1.

Let us also observe that if processes ρ 1 , ρ 2 , and ρ 3 are executed by a qualified entity Ψ 1 with knowledge of reference a ′ 3 , the returned set of column vectors b α , b ′ α , b ′′ α , and so forth, are clearly linked:

b α =   2 23 20   , b ′ α =   9 1 17
  = r 1 b α +   27 2 24   , . . .
Conversely, if we assume that processes ρ 1 , ρ 2 , and ρ 3 were executed by a malicious adversary Ψ 2 who is trying to link the shares returned by either α or β , for tracking purposes, but not having access to the column vector reference a ′

SIMULATION AND EXPERIMENTAL RESULTS

In Section 3, we have seen the formalization of our proposal for the reconstruction of a pre-distributed secret once a sufficient number of shares are collected. We present in this section the results obtained with an experimental setup that simulates EPC Gen2 adapted shares generation and reconstruction of secrets. Our prototype system allows to experiment the exchange of shares with a regular EPC Gen2 reader and simulated Gen2 tags. The objective of this setup is to demonstrate the practical viability of our proposal.

Experimental Setup

Figure 3 pictures our experimental setup. It is based on the execution of standard EPC Gen2 inventory queries, but enabled with TIDs that are enhanced by our proposed threshold cryptosystem, between a regular EPC Gen2 reader and several Gen2 tag instances simulated by the IAIK UHF demo tag [START_REF]UHF RFID Demo Tag[END_REF].

The IAIK UHF demo tag is a programmable device intended for developing new extensions to the EPC Gen2 standard. The demo tag consists of an antenna, an RF front-end, a programmable microcontroller, and a firmware library. The antenna captures the energy emitted by the reader and powers up the RF front-end of the tag. The RF front-end demodulates the information encoded in the signal. The encoded data is processed by the programmable microcontroller to compute a response. To compute the response, the programmable microcontroller executes a software implementation of the EPC Gen2 protocol, implemented in the firmware library. The response is then modulated by the RF front-end and backscattered to the reader. More details can be obtained in [START_REF]UHF RFID Demo Tag[END_REF][START_REF] Aigner | BRIDGE -Building Radio frequency IDentification for the Global Environment[END_REF]. Our share renewal scheme has been implemented in ANSI C using the Crossworks IDE for AVR from Rowley Associates [START_REF] Crossworks | Crossworks v1.4 and v2.0 for AVR[END_REF]. The theoretical construction detailed in Section 3 has been adapted to be executed over the Atmel AVR ATmega128 [37] microcontroller of the IAIK UHF demo tag. The AT-mega128 is an 8-bit microcontroller based on the AVR architecture. It has (a) (b)

Experimental Setup. In (a), we can see the CAEN A829EU Reader, the AVR JTAG MKII Programmer, the IAIK Graz UHF Demo Tag, and some regular EPC Gen2 tags. In (b), we can see the Java graphical front end that summarizes the process of collecting the secret shares and reconstruction of secrets.

32 registers of 8-bits that can act as the destinations of standard arithmetic operations. In addition, the ATmega128 microcontroller contains 128KB of flash memory and 4KB of data memory that can be addressed by three independent registers of 16-bits. Since the response of inventory queries is a mandatory operation specified in the EPC Gen2 protocol, an existing response function implemented for the ATmega128 microcontroller is already included in the original firmware stored on the IAIK UHF demo tag. By using the Crossworks IDE, we code and merge the new functionality with the general firmware library to adapt the existing inventory response process to the renewal scheme of shares. The AVR JTAG MKII programmer [37] is used to transfer and to debug the updated firmware merged with the adapted inventory routine. On the reader side, the short-range reader CAEN A829EU [START_REF] Caen Rfid | [END_REF] emits the inventory queries. The reader is controlled by a back end computer over a USB serial port and a Java application. The Java application is in charge of generating the inventory queries and processing the reconstruction of secrets.

Collection of Shares and Reconstruction Rates

Four different populations of EPC Gen2 tags are simulated and tested. All four simulations are built according to the item-level inventory scenarios reported in [START_REF]EPC Item Level Tagging Joint Requirements Group[END_REF][START_REF] Juels | Unidirectional Key Distribution Across Time and Space with Applications to RFID Security[END_REF]. Our objective is to show how our construction can be used in order to maximize the item traceability rate at the upper levels of a supply chain, i.e., at the manufacturer, distributor and retailer sides, while minimizing the traceability rate at the lower levels of a supply chain, i.e., at the consumer side. The study presented in [START_REF]EPC Item Level Tagging Joint Requirements Group[END_REF] shows that items that are initially assembled and tagged together within large collections at the manufacturer side, i.e., top level of the supply chain, get progressively dispersed into very small subsets when they reach the bottom level of the supply chain, i.e., the consumer side. Two appropriate item examples analyzed in [START_REF]EPC Item Level Tagging Joint Requirements Group[END_REF] are personal hygiene tools and pharmaceuticals products. According to [START_REF]EPC Item Level Tagging Joint Requirements Group[END_REF], personal hygiene tools like, for instance, razor blades, are initially assembled and tagged together at the manufacturer side of the supply chain in large populations of more than 6,000 tagged items. They are later dispersed in the supply chain until being picked up by consumers in groups of less that five items. Similarly, for pharmaceutical items assembled initially in large quantities of more than 7,000 tagged items at the manufacturer side, we should only expect that no more than six items from the initial population can end in possession of a single consumer at the same time. In accordance with these observations, we simulate four different populations of EPC Gen2 tagged items. The tags of each population are initialized with four independent sets of secret sharing schemes constructed according to our proactive threshold secret sharing scheme in GF(2 5 -1). More precisely, we initialize the tags of the first population with a (13, 24) scheme that produces tag inventory responses of 120 bits; the tags of the second population with a (10, 18) scheme that produces tag inventory responses of 90 bits; the tags of the third population with a (7, 12) scheme that produces tag inventory responses of 60 bits; and the tags of the fourth population with a (5, 8) scheme that produces tag inventory responses of 40 bits. Figure 4(a)-(d) pictures the average and the 95% confidence intervals of the reconstruction rates obtained with the collection of less than 35 shares from each simulated population. We recall that these simulations take into account the evaluation reported in [START_REF]EPC Item Level Tagging Joint Requirements Group[END_REF]. We, therefore, consider the upper bounds of five to thirteen items as the sizes of groups of items picked up by consumers (i.e., lower level of a supply chain). Above these bounds, it is straighforward that authorized readers at the store, warehouse or manufacturer facilities will always reach the necessary threshold to reconstruct the secret and access the appropriate TIDs. For each experimental test of each population, the inventory query emitted by the EPC Gen2 reader is responded by exactly m random tags, where 35 < m < 0. We recall that the use of a (t,n) scheme means that of the n available shares, we need to collect, at least, t different shares to successfully reconstruct the distributed secret. Each population of tags is initialized by randomly allocating shares from each of threshold scheme. Each interrogation is executed 100 times with random series of simulated tags. The results we show are therefore the average and 95% confidence intervals computed after each series of interrogations.

The results confirm that while the reconstruction rate minimizes the traceability of tagged items as soon as these items get dispersed in small quantities on the consumer side (amount of tagged items below quantities of less than twenty tags), it guarantees the identification of these items at the upper levels of the supply chain (amount of tagged items above quantities of more than thirty tags). From these results, we may also conclude that the compact size of shares of all four schemes are appropriate enough to fit on the inventory responses suggested on the EPC standard. Note that the resulting inventory responses that are containing the shares (i.e., 120 bits for the first population; 90 bits for the second population; 60 bits for the third population; and 40 bits for the fourth population) do successfully fit within the maximum response size of 528 bits suggested by EPCglobal in [START_REF] Epcglobal | EPC Radio-frequency Identity Protocols Class-1 Generation-2 UHF[END_REF].

CONCLUSIONS

We presented a proactive secret sharing procedure to provide consumer privacy and distribution of secrets. Our solution addresses the eavesdropping, rogue scanning, and tracking threats. The main properties of our approach are: (1) low-cost share renewal with secret preservation and without a need of synchronization; (2) compact size of shares; (3) secret sharing construction that guarantees strong security; (4) reconstruction of a secret does not require the identity of the shareholders. We have also presented the implementation of a practical experiment with our proposed cryptosystem in a real EPC Gen2 scenario. By means of a compatible Gen2 reader, and a programmable Gen2 tag implementing our proactive share renewal process, we have shown that a standard EPC Gen2 reader can reconstruct an appropriate pre-disctributed secret dispersed over a set of Gen2 tags. The set of tags communicate the renewed shares to the reader by using a standard inventory response operation, enhanced by our proposed proactive share renewal.



  We can now observe that the orthogonal projector obtained by applying Equation (4) to either B 1 , B 2 , or B 3 is equivalent to the orthogonal projector obtained by applying Equation (4) to matrix A: , P B 1 = P B 2 = P B 3 =



  process ρ 3 . Finally, the process ρ 1 assembles with b α , b β the reconstruction matrix B 1 ∈ Z 3×2 31 ; the process ρ 2 builds with b ′ α , b ′ β the reconstruction matrix B 2 ∈ Z 3×2 31 ; and the process ρ 3 produces with b ′′ α , b ′′ β the reconstruction matrix B 3 ∈ Z 3×2 31 :We observe that the orthogonal projectors obtained by applying Equation (4) to matrices B 1 , B 2 , and B 3 are identical to the orthogonal projector obtained by applying Equation (4) to matrix A:, P B 1 = P B 2 = P B 3 = Notice that each matrix B i = [b i1 , b i2], s.t. i ∈ {1 . . . 3}, can be decomposed as follows:

3 ,

 3 the returned set of column vectors b α , b ′ α , and b ′′ α , as well as column vectors b β , b ′ β , and b ′′ β , are viewed as unlinked.

4

 4 FIGURE 4 Simulation results. (a) First population results with a[START_REF] Dimitriou | A Lightweight RFID Protocol to Protect Against Traceability Cloning Attacks[END_REF][START_REF] Juels | Unidirectional Key Distribution Across Time and Space with Applications to RFID Security[END_REF] proactive threshold secret sharing scheme of shares; (b) second population results with a (10,18) scheme; (c) third population results with a (7,12) scheme; (d) fourth population results with a (5,8) scheme.

Acknowledgments

-This work has been supported by the Natural Sciences and Engineering Research Council of Canada (NSERC), the Mathematics of Information Technology and Complex Systems (MITACS), the Spanish Ministry of Science (grants TSI2007-65406-C03-03 E-AEGIS and CONSOLIDER-INGENIO 2010 CSD-2007-00004 ARES), and the Institut TELECOM (Futur et Ruptures program).

We would also like to thank J. Melia-Segui and J. Herrera-Joancomarti for fruitful discussions on the experimental setup used in Section 4.


Alexandre Lebouvier

Clarisse Delalondre

Francois Fresnet

Valerie Boch

Vandad Rohani

Francois Cauneau

Laurent Fulcheri
email: laurent.fulcheri@mines-paristech.fr

Three-Dimensional Unsteady MHD Modeling of a Low Current -High Voltage Non-Transferred DC Plasma Torch Operating With Air

Keywords: MHD unsteady modeling, Low current DC arc discharge, non-transferred plasma torch

We present in this paper the MHD modeling of a DC plasma torch operating with air under very peculiar high voltage-low current conditions. The model developed is 3D, timedependent and assumes LTE. The study has been carried out considering an axial injection of air with flow rates varying in the range 0.16 -0.5 g/s and currents varying in the range 300 -600 mA.

The numerical modeling has been developed using Code Saturne, CFD software developed by EDF R&D which is based on co-located finite volume.

After a detailed description of the model, results are presented, analyzed and discussed. The influence of current and air flow rate over the arc characteristics are studied in terms of temperature, velocity, electrical potential, Joule heating and arc root motion. Regarding numerical issues, the MHD modeling of low current -high voltage arc discharge is particularly tricky since, below 1 A, the self-induced magnetic field becomes negligible and the convection effects induce a highly irregular and unstable motion of the arc column. However, despite these difficulties, the numerical model has been successfully implemented. Numerical results have shown good correlation and good trends with experimental ones despite a discrepancy probably due to the LTE assumption. The model gave fruitful and significant information on parameters that could hardly be obtained experimentally. This preliminary work is likely to open the way towards a better understanding of low current arc discharges, which technologies are currently encountering an important development in many application fields.

I. INTRODUCTION

F OR many years, numerical modeling became a unique and powerful tool for the development and the optimization of plasma processes. With computing power improvement, models get more and more sophisticated and close to real conditions thus allowing today the simulation of three-dimensional and unsteady systems with detailed geometry description. Numerous fluid models have been reported in the literature mostly on direct current (DC) arc discharge modeling using high current above 100 A [START_REF] Gleizes | Thermal plasma modelling[END_REF], [START_REF] Trelles | Arc Plasma Torch Modeling[END_REF] including: free burning arcs [START_REF] Chemartin | Three dimensional simulation of a DC free burning arc. Application to lightning physics[END_REF], [START_REF] Gonzalez | Numerical modelling of an electric arc and its interaction with the anode: part II. The three-dimensional model-influence of external forces on the arc column[END_REF], transferred arcs [START_REF] Chau | Modeling and experimental validation of a 1.2 MW DC transferred well-type plasma torch[END_REF]- [START_REF] Douce | Numerical modelling of an anodic metal bath heated with an argon transferred arc[END_REF] and non-transferred arcs [START_REF] Moreau | Modeling the restrike mode operation of a DC plasma spray torch[END_REF]- [START_REF] Selvan | Numerical and experimental studies on DC plasma spray torch[END_REF]. Many DC arc plasma torches models have been developed for different applications [START_REF] Boulos | Thermal Plasmas-Fundamentals and Applications[END_REF] as: cutting [START_REF] Zhou | Comparative study of turbulence models on highly constricted plasma cutting arc[END_REF], welding [START_REF] Murphy | The effects of metal vapour in arc welding[END_REF], deposition and spraying [START_REF] Moreau | Modeling the restrike mode operation of a DC plasma spray torch[END_REF]- [START_REF] Selvan | Numerical and experimental studies on DC plasma spray torch[END_REF], [START_REF] Li | Three dimensional modeling of the plasma spray process[END_REF], [START_REF] Baudry | Numerical modeling of a DC non-transferred plasma torch: Movement of the arc anode attachment and resulting anode erosion[END_REF], steelmaking [START_REF] Douce | Numerical modelling of an anodic metal bath heated with an argon transferred arc[END_REF], waste disposal [START_REF] Freton | Prediction of the cathodic arc root behaviour in a hollow cathode thermal plasma torch[END_REF] and ultra fine particle production [START_REF] Tang | Threedimensional modelling of a dc arc plasma in a twin-torch system[END_REF]. These thermal plasmas generally work under low-medium voltage between 30 and 500 V [START_REF] Trelles | Arc Plasma Torch Modeling[END_REF]. In nontransferred plasma torch, the arc behavior can generally be classified into three operation modes : steady, takeover and restrike modes respectively [START_REF] Trelles | Arc Plasma Torch Modeling[END_REF]. The steady mode corresponds to a fixed position of the arc attachment whereas the takeover mode corresponds to a quasi-periodic fluctuation of the voltage drop and hence of the arc movement. The restrike mode is characterized by a highly unstable movement of the arc where the reattachment phenomenon plays a significant role. The first models were two-dimensional (2D) or 2D axisymmetric thus, the vortex inflow and the arc root attachment could not be perfectly reproduced [START_REF] Douce | Numerical modelling of an anodic metal bath heated with an argon transferred arc[END_REF], [START_REF] Kaddani | Threedimensional modeling of unsteady high-pressure arcs in argon[END_REF], [START_REF] Freton | Comparison between a two-and a three-dimensional arc plasma configuration[END_REF]. With the improvement of computing power, many research teams moved towards threedimensional models which are closer to real configurations [START_REF] Gonzalez | Numerical modelling of an electric arc and its interaction with the anode: part II. The three-dimensional model-influence of external forces on the arc column[END_REF]. Most of these models assume the local thermodynamic equilibrium (LTE) except e.g. Trelles et al. [START_REF] Trelles | Non-equilibrium modelling of arc plasma torches[END_REF] and Park et al. [START_REF] Park | Two-Dimensional Numerical Modeling of Direct-Current Electric Arcs in Nonequilibrium[END_REF] who developed non-LTE (NLTE) models in argon. The calculations are often realized in a simple gas, mainly argon [START_REF] Gonzalez | Numerical modelling of an electric arc and its interaction with the anode: part II. The three-dimensional model-influence of external forces on the arc column[END_REF], [START_REF] Li | Three dimensional modeling of the plasma spray process[END_REF], [START_REF] Tang | Threedimensional modelling of a dc arc plasma in a twin-torch system[END_REF], [START_REF] Freton | Comparison between a two-and a three-dimensional arc plasma configuration[END_REF] but also nitrogen [START_REF] Chau | Modeling and experimental validation of a 1.2 MW DC transferred well-type plasma torch[END_REF], pure oxygen [START_REF] Zhou | Comparative study of turbulence models on highly constricted plasma cutting arc[END_REF], or a gas mixture: air [START_REF] Freton | Prediction of the cathodic arc root behaviour in a hollow cathode thermal plasma torch[END_REF], Ar/H 2 [START_REF] Moreau | Modeling the restrike mode operation of a DC plasma spray torch[END_REF], Ar/N 2 [START_REF] Selvan | Three-Dimensional Numerical Modeling of an Ar-N2 Plasma Arc Inside a Non-Transferred Torch[END_REF]- [START_REF] Selvan | Numerical and experimental studies on DC plasma spray torch[END_REF]. The main focuses of arc plasma torch modeling are the arc root displacement and attachment [START_REF] Baudry | Numerical modeling of a DC non-transferred plasma torch: Movement of the arc anode attachment and resulting anode erosion[END_REF], [START_REF] Freton | Prediction of the cathodic arc root behaviour in a hollow cathode thermal plasma torch[END_REF], and the metal electrode evaporation [START_REF] Murphy | The effects of metal vapour in arc welding[END_REF], [START_REF] Baudry | Numerical modeling of a DC non-transferred plasma torch: Movement of the arc anode attachment and resulting anode erosion[END_REF]. Although the laminar flow regime is commonly used for modeling plasma torch [START_REF] Blais | Three-dimensional numerical modelling of a magnetically deflected dc transferred arc in argon[END_REF], [START_REF] Moreau | Modeling the restrike mode operation of a DC plasma spray torch[END_REF], [START_REF] Selvan | Comparisons Between Two Different Three-Dimensional Arc Plasma Torch Simulations[END_REF], [START_REF] Freton | Prediction of the cathodic arc root behaviour in a hollow cathode thermal plasma torch[END_REF], [START_REF] Bernardi | Comparison of different techniques for the FLUENT c -based treatment of the electromagnetic field in inductively coupled plasma torches[END_REF], literature reports few studies using turbulent k-ǫ model [START_REF] Li | Three dimensional modeling of the plasma spray process[END_REF], [START_REF] Tang | Threedimensional modelling of a dc arc plasma in a twin-torch system[END_REF], Rij-ǫ and low-Re models [START_REF] Delalondre | Turbulence Modelling in Electric Arc[END_REF], comparing laminar to standard k-ǫ model [START_REF] Bauchire | Modeling of a DC plasma torch in laminar and turbulent flow[END_REF], or comparing several turbulent models (standard k-ǫ, RNG k-ǫ, realizable k-ǫ) [START_REF] Zhou | Comparative study of turbulence models on highly constricted plasma cutting arc[END_REF].

For several years, researches carried out at CEP (Center for Energy and Processes, MINES ParisTech) in the field of hydrocarbons reforming led to the development of a very peculiar DC plasma torch technology operating under low current -high voltage conditions: typically in the range 400 mA and 2 kV respectively. This type of torch is a non-thermal atmospheric plasma torch and, in some operating points, can be assimilated to a LCGA (Low Current Gliding Arc) as we work under 1 A [START_REF] Fridman | Non-thermal atmospheric pressure discharges[END_REF]. The high degree of non-equilibrium of this torch is promising for chemical applications to support chemical reactions selectively [START_REF] Fridman | Non-thermal atmospheric pressure discharges[END_REF], but suffer from a lack of knowledge of its physics and chemistry. Our researches led to the development of basic knowledge on both experimental [START_REF] Petitpas | Ethanol and E85 Reforming Assisted by a Non-thermal Arc Discharge[END_REF], [START_REF] Rollier | Experimental Study on Gasoline Reforming Assisted by Nonthermal Arc Discharge[END_REF] and numerical [START_REF] Gonzalez-Aguilar | Three Stages Modeling of n-Octane Reforming Assisted by a Nonthermal Arc Discharge[END_REF], [START_REF] Rollier | Thermodynamics and kinetics analysis of gasoline reforming assisted by arc discharge[END_REF] sides. Unfortunately, whereas many studies have been dedicated to high current plasma discharges, the literature on MHD modeling of low current -high voltage plasma torches is very poor and the domain remains almost unexplored. One of the main reason explaining the lack of studies dedicated to this purpose is probably due to the fact that, at low current, the self-induced magnetic field becomes negligible and the convection effects induce a highly irregular motion of the arc column. These arc motion instabilities imply very unstable physical phenomena. Thus, one of the most challenging issue in the field of plasma modeling for the next coming years is to be able to simulate very low-current arc discharges owing to important numerical instabilities deriving from physical instabilities.

Different numerical models developed for plasma reforming addressing thermodynamics and kinetics issues have been reviewed in Ref. [START_REF] Petitpas | A comparative study of non-thermal plasma assisted reforming technologies[END_REF]. On the other hand, Bromberg et al. [START_REF] Bromberg | CFD modeling of Plasmatron Methane Reformer[END_REF] and Fridman et al. [START_REF] Chirokov | Reversevortex plasma stabilization: experiments and numerical simulation[END_REF] used computational fluid dynamics (CFD) codes to simulate the hydrodynamic flows in the torch with considering neither kinetic reactions nor MHD equations. In the first paper, the plasma is not considered and the study has been focused on the methane/air mixing at the location of the plasma. In the second one, the plasma is considered as a source term of power and the reverse vortex flow influence on the plasma is investigated.

In the first part of this paper, the model assumptions, boundary conditions, equations and simulation parameters are successively presented. In the second part of this paper, the low current arc discharge behavior is presented, analyzed and discussed in terms of temperature, velocity, electrical potential, Joule heating and arc root velocity with a special attention being paid to evaluate the influence of arc current and air flow rate.

II. MATHEMATICAL MODEL

A. Assumptions

The three-dimensional model studied is based on the following main assumptions:

• The plasma is considered as a single continuous fluid (Air).

• The plasma is optically thin and at Local Thermodynamic Equilibrium (LTE).

• The gas is treated as semi-compressible because the Mach number is lower than 0.2.

• The gas flow is laminar and time-dependent.

• Gravitational effects are taken into account (-x direction). Let us point out that whereas high current arcs or high power-density plasmas [START_REF] Tanaka | Thermally and chemically non-equilibrium modelling of Ar-N 2 -H 2 inductively coupled plasmas at reduced pressure[END_REF] are commonly assumed to be in LTE. For low current or low power density arcs, which has a high degree of non-equilibrium, this hypothesis is generally not correct. It has been demonstrated that this assumption can imply great differences in the arc behavior and characteristics [START_REF] Trelles | Non-equilibrium modelling of arc plasma torches[END_REF]. However, even if this is a strong assumption, the LTE assumption has been considered due to highly challenging numerical issues linked with the low current. The objective of this study was to get preliminary information on transport phenomena including mass, momentum, and energy as well as on the gas flow and temperature fields. Obviously, it would be interesting subsequently, in a later phase, to develop more sophisticated NLTE models such as two-temperature models.

B. Governing equations

The model can be defined by the set of following equations based on fluid dynamics Navier-Stokes and Maxwell electromagnetic equations respectively. All the equations are written in a Cartesian system (Ox, Oy, Oz). In thermal plasmas, which are assumed to be electronically neutral, the only resulting force applied on a volume element is the Lorentz Force (N/m 3 ):

F L = J × B ( 1
)
where J and B correspond to the electric current density vector (A/m 2 ) and the magnetic field (T) respectively. The generalized Ohm's law approximation commonly used in electric arc modeling [START_REF] Trelles | Arc Plasma Torch Modeling[END_REF] is :

J = σE (2)
where σ and E represent the electric conductivity (S/m) and the electric field (V/m) respectively. Neglecting the Neumann's electromotive field, the electrical field can be linked to the electrical potential φ(V) by the following relation:

E = -∇φ (3)
In absence of external magnetic field, the electric potential equation is approximated to:

-∇ • (σ∇φ) = 0 (4)
The magnetic field and the magnetic vector potential (T.m) are related by the following relation:

B = ∇ × A (5)
The Maxwell-Ampere equation is:

∆A = -µ 0 J (6)
where µ 0 is the vacuum permeability. The set of Navier-Stokes fluid conservation equations is presented below.

Mass conservation

∂ρ ∂t + ∇ • (ρV) = 0 (7)
where ρ and V are the mass density and the velocity vector.

Energy conservation equation

∂ρh ∂t + ∇ • ρhV -∇ • λ C p ∇h = J • E -S rad (8)
where h, λ, C p and S rad are the gas enthalpy, the thermal conductivity, the specific heat and the radiation losses respectively. Momentum equation

∂ρV ∂t + ∇ • (ρV ⊗ V) = -∇p + ∇ • τ + J × B + ρg (9)
where p, τ and g are the pressure, the shear stress tensor and the gravity vector respectively.

The study was carried out for an arc discharge operating with air. The properties data have been taken from [START_REF] D'angola | Thermodynamic and transport properties in equilibrium air plasmas in a wide pressure and temperature range[END_REF] for enthalpy, density, specific heat, electrical conductivity, dynamic viscosity and thermal conductivity. For the source radiation, the net emission coefficients have been taken from [START_REF] Naghizadeh-Kashani | Net emission coefficient of air thermal plasmas[END_REF].

The resistive MHD equations have been solved using Code Saturne v. 1.3, developed by EDF R&D. The Code Saturne electric arc module has been used. This CFD software is based on co-located finite volume [START_REF] Archambeau | Code Saturne: a finite volume code for the computation of turbulent incompressible flows[END_REF]. Code Saturne's capabilities enable the handling of incompressible and expandable flows with heat transfer using SOLU convective scheme (Second Order Linear Upwind) and SIM-PLEC pressure correction scheme (Semi-Implicit Method for Pressure-Linked Equations, Consistent) [START_REF] Van Doormaal | Enhancements of the SIMPLE method for predicting incompressible fluid flows[END_REF].

C. Computational grid and boundary conditions

The torch geometry corresponding to the real experimental set-up, is presented on Fig. 1. It is composed of two separated zones: the plasma zone and the post-discharge zone with 70 mm / 4 mm and 100 mm / 11 mm length / inner radius respectively. Let us point out that, under real reactive operating conditions, the plasma zone is the part where the arc plasma really takes place whereas the post-discharge zone is an active or passive zone, located downstream the plasma zone where, reactions ignited in the plasma zone continue to take place.

Usually, in DC plasma torches, the tip electrode is the cathode. In our case, the system operates under inversed polarity. The tip electrode, on which a positive electrical potential is applied, being the anode whereas the cylinder electrode is grounded and acts as the cathode. Indeed, previous works have demonstrated this configuration was the most suitable for low current systems [START_REF] Rollier | Thermodynamics and kinetics analysis of gasoline reforming assisted by arc discharge[END_REF].

The grid mesh, shown on Fig. 2, is exclusively composed of hexa-elements and contained 339 000 cells. Cubic cells of ∆x = 0.044 mm compose the center of the computational grid. The mean axial and radial spacing grid are 0.7 mm and 0.21 mm respectively. These values are very similar to the ones used by Freton et al. [START_REF] Freton | Prediction of the cathodic arc root behaviour in a hollow cathode thermal plasma torch[END_REF]. This very fine grid was chosen to compute correctly the gradients of all transport variables into the arc column. The grid is axially finer near the anodic and the cathodic arc root attachment position.

The boundary conditions are detailed in Table I. The atmospheric pressure is imposed at the outlet. The temperature of the injected air is 600 K. The walls are considered adiabatic. At the first time step, an atmospheric pressure and a temperature of 600 K is applied to the whole domain.

The potentiel to the cathode is imposed to 0 V, and the potential to the anode is calculated from the previous time step. The voltage drop in the sheath is not included in the model. The potential results correspond to the voltage drop in the arc column. In our case, the electrical current is very low, so it is difficult to assume a value of electrical voltage drop in the cathodic and anodic sheaths. For high intensity electric arcs, authors generally used an imposed value of 3 to 10 V for the potential drop in the sheaths [START_REF] Gonzalez | Numerical modelling of an electric arc and its interaction with the anode: part II. The three-dimensional model-influence of external forces on the arc column[END_REF], [START_REF] Moreau | Modeling the restrike mode operation of a DC plasma spray torch[END_REF], [START_REF] Sanders | Measurement of anode falls and anode heat-transfer in atmospheric-pressure high-intensity arcs[END_REF].

A null-flux condition is applied for the heat conduction at the electrodes. Even if this assumption is strong, it can be justified by the fact that the electrodes are not water-cooled so the losses at the electrodes are not as significant as they are in thermal torches.Nevertheless, in a next step model, we will take into account the electrode sheath by introducing a global model based on the energy exchange in the cathodic and anodic sheath [START_REF] Moreau | Modeling the restrike mode operation of a DC plasma spray torch[END_REF], [START_REF] Trelles | Modelling of the arc reattachment process in plasma torches[END_REF].

D. Simulation parameters

The time step ∆t is set to 10 µs as Chemartin et al. [START_REF] Chemartin | Three dimensional simulation of a DC free burning arc. Application to lightning physics[END_REF], [START_REF] Chemartin | 3D simulation of electric arc column for lightning aeroplane certification[END_REF]. This time step allows obtaining fair plasma characteristics with a reasonable computation time. Freton et al. [START_REF] Freton | Prediction of the cathodic arc root behaviour in a hollow cathode thermal plasma torch[END_REF] and Park et al. [START_REF] Park | Three-Dimensional Modeling of Arc Root Rotation by External Magnetic Field in Nontransferred Thermal Plasma Torches[END_REF] have used a 5 µs time step after a parametric study.

At the first time step, a hot temperature channel is artificially created to ignite the arc such as proposed, for example, by Moreau et al. [START_REF] Moreau | Modeling the restrike mode operation of a DC plasma spray torch[END_REF]. The channel has a 0.8 mm radius. The temperature imposed in the initial channel is 6000 K for starting calculation (cf. Fig. 3). The initial anode voltage drop is set to 1000 V. For the first 50 time steps i.e. 500 µs, the inflow is 1 % of the total inflow in order to allow the correct arc ignition. Subsequently, for the 150 following time steps, a linear velocity ramp is applied at the inlet to increase the inflow from 1 to 100 %.

The current target is set between 300 mA and 600 mA. The current is mesured by the integral of the current density in a z-normal plane S z (cf. Eq. 10), located 1.5 mm after the anode. At each time step, the anodic electrical potential φ n is fitted to maintain the current target from Eq. 11.

I = Sz J • dS (10)
φ n = I target I mesured φ n-1 (11)
2000 time steps have been considered, corresponding to a numerical time of 20 ms. The calculations have been realized on 8 processors (Intel Xeon 2.66 GHz) for 30 hours. That corresponds to a calculation time of around 55 s for each time step. The number of under iteration is around 2300 for pressure, 1300 for potential vector, and under 100 for enthalpy and velocity.

III. RESULTS

A. Results for I = 400 mA and Q air = 0.32 g/s

In this section, we present results for an electrical current of 400 mA and an air flow rate of 0.32 g/s.

Fig. 4 represents the evolution of the temperature. We can observe that the arc root moves along the same plane as the ignition plane (Oy, Oz) due to the axial injection. Once the arc root has reached the cathode tip, it remains attached to the cathode tip at approximately 1.30 mm above the exit. The arc root path follows the cathode shape at the exit. Indeed, in this zone, the velocity magnitude is close to zero (cf. Fig 5) because of the dead zone created by the sudden section enlargement. In this zone, the cathodic arc root is less disturbed by the gas inflow and the dead zone constitutes a favourable arc root attachment point.

The magnetic field is produced by the local curvature of the current path and the self-induced magnetic field. The arc dynamics inside the torch is the result of the imbalance between the electromagnetic forces and the flow drag. As expected, the magnetic field magnitude is very low. Its maximum value reaches 0.7 mT at the cathode tip where tip effect is the highest. As a result, the Lorentz forces are also very low (up to 120 N/m 3 ). Around the arc core, the mean value of the magnetic field magnitude is around 0.1 mT. These very low values explain mainly the instabilities observed in low current plasma torches compared to high current plasma torches, much more stable. The magnetic field is so low that the arc column is only stabilized by the hydrodynamic effects, especially by the wall or vortex stabilization commonly used in low current Fig. 4. Representation of the temperature field displaying the evolution of the arc shape in a cross-section for three specific times (t = 2, 5 and 20 ms). I = 400 mA. Q air = 0.32 g/s. -high voltage torch.

The maximum velocity is 67 m/s and is located at the cathode tip. The velocity magnitude distribution is not symmetric due to the arc curve. Once the inlet velocity has reached its maximum, the arc root velocity can be estimated to 15.2 m/s along the cathode inner surface. This calculated velocity is higher than the average velocity in the nozzle (w = 11 m/s) because of the thermal effects. On the other hand, the Lorentz force is so low that the induced convection effect can be neglected.

The anodic and cathodic arc root temperature values are around 5500 and 5100 K respectively. This temperature corresponds to the temperature of the plasma at the interface between the electrode sheath region and the arc column. The arc column is very stable inside the cathode essentially due to the axial injection and the laminar flow assumption. One can observe on Fig. 6 that the arc core temperature in the cathode reaches 5300 K with a significant temperature gradient between the center of cathode and its borders (∼ 1500 K/mm). Fig. 7 presents a real photographs of the plasma from experiments in close conditions. This picture confirms the phenomenon predicted by the simulation, previously presented on Fig. 4. The plasma plume observed by the experiment has the same shape as that numerically obtained. The cathodic arc root becomes attached to a point located above the exit.

We can also observe on Fig. 7 that the arc root radiations are higher than the arc core radiations. This difference is not observed by the simulation (cf. Fig. 8), the arc root is as bright as the arc core. This phenomenon can probably be explained by the fact that the model does not take into account the electrode metal evaporation. Indeed, the cathode being composed of stainless steel, it is composed by more than 60 % of iron and it is well known that iron net emission is very high [START_REF] Menart | Net emission coefficients for argon-iron thermal plasmas[END_REF]. We can also see on this picture that the plasma plume length, which is about 55 mm numerically, is close to the experimental one, which is around 45 mm.

As assumed by Selvan and coworkers [START_REF] Selvan | Numerical and experimental studies on DC plasma spray torch[END_REF], the arc core radius is the radial distance along the anode from the centerline to the point where the current density is zero. In practice, the zero current density is supposed being achieved when the current density is 2 orders of magnitude smaller than the maximum current density which is 115 A/cm 2 . The current density is shown on Fig. 9. By this technique, the arc core radius can be estimated to 0.5 mm. This information is hardly accessible by the experiment.

On the other hand, the radius of the most emitting zone, has been experimentally estimated, to the naked eye, to 1 mm at the torch exit. It can be compared to the radius value of the most emitting zone found by the simulation. From Fig. 8, this value is around 0.83 mm, so demonstrating the good correlation between the numerical and the experimental results.

The voltage drop and the inlet velocity ramp are presented on Fig. 10. After 500 µs, the arc is blown by the air flow rate. The higher the arc length, the higher the voltage. Indeed, it is well known that the arc voltage varies quasi-linearly with the arc length. We can observe few low instabilities (f = 20 kHz) in the range 1.3 -6 ms. These instabilities are probably numerical instabilities due to the sudden high stretch of the low current arc.

Once the arc root has reached the cathode tip, the mean voltage is equal to 2.7 kV and corresponds to a deposited plasma power of 1.1 kW. For that case, a typical experimental oscillogram is displayed on Fig. 11 and gives a mean voltage 30 % (1.9 kV) lower than the simulated one. As demonstrated by Trelles et al. [START_REF] Trelles | Non-equilibrium modelling of arc plasma torches[END_REF], the total voltage drop for a NLTE model of a high current plasma torch can be 50 % lower than for a LTE model. Indeed, for NLTE models, the electron temperature and the electrical conductivity are higher than in LTE models. Hence the resistance to the current flow decreases and leads to a decrease of the voltage drop decrease. Thus, the LTE assumption is probably the main reason of this discrepancy, even if some other factors could intervene e.g. the metal vapor near the electrodes.

After the arc has reached the tip, we can also observe numerically, on Fig. 10 the apparition of voltage instabilities. These instabilities arise from a coming and going motion of the arc bend, which change the arc length and thus the voltage drop. Instabilities are also observed experimentally and can be seen on Fig. 11. The normalized FFT of both experimental and simulated voltage oscillograms (for t > 10 ms) are presented on Fig. 12, in a semilogarithmic scale. The experimental spectrum is broad, whereas the simulated is thinner but shows several frequency peaks. The oscillation frequency peaks are in the same range of frequencies. The highest peak is for a frequency of 1.85 kHz and 3.2 kHz for the numerical and experimental spectrum, respectively. Let us point out that the 43 kHz peak directly comes from the power supply frequency which is based on a resonance converter. Hence the instabilities over 10 ms could be due mainly to hydrodynamic forces. Even if the arc root stays at a fixed position (similar to the steady mode) these periodic fluctuations, coming from the arc plume length, remind the takeover mode [START_REF] Duan | Arc instabilities in a plasma spray torch[END_REF].

The Joule heating (J • E) balances the radiation losses and the convective cooling effect driven by the Lorentz forces !" # Fig. 12. Comparison between the normalized FFT of the experimental and the numerical voltage oscillograms, plotted in a semilogarithmic scale. I = 400 mA. Q air = 0.32 g/s. [START_REF] Tang | Threedimensional modelling of a dc arc plasma in a twin-torch system[END_REF], almost negligible in our case. In low current plasma torches, the convective cooling effect mostly comes from the gas inflow. The Joule heating is strong (27 kW/cm 3 ) in the anodic zone, where the current density is maximum, but lower than 1 kW/cm 3 in the rest of the arc due to the low current used.

B. Influence of the current

The current influence has been investigated in the range of 300 mA -600 mA for an air flow rate of 0.32 g/s.

Fig. 13 shows photographs of the plasma plume vs. electrical current. The higher the current, the higher the power and the stability. Indeed, for I = 300, 400 and 600 mA, the experimental deposited power is 640, 770 and 950 W respectively.

Fig. 14 shows the numerical oscillograms for I = 300, 400 and 600 mA and on Figs. 11 and 15 the associated experimental oscillograms. On Fig. 14, we can observe that, before the arc root reached the cathode tip (t < 10 ms), the higher the current, the lower the instabilities magnitude.

The mean voltage drop, showed on Fig. 16, and the standard deviation for the three cases, both experimental and simulated, have been reported on Table II. The model follows the experimental voltage trends: the higher the current, the lower the voltage but the higher the deposited power. This is a typical characteristic of non-thermal plasmas. For I = 300, 400 and 600 mA, the numerical deposited power is 1000, 1100 and 1200 W respectively.

Moreover, one can observe that the higher the current, the higher the instabilities magnitude. The oscillation frequencies are in the same frequency range for the three currents.

We can also notice that the lower the current, the higher the difference between simulated and experimental mean voltage. This phenomenon is probably due to the non-equilibrium effects which increase when decreasing the current.

The radial profile of the temperature in the middle of the cathode is shown on Fig. 17. The higher the current, the higher the temperature and the wider the arc core temperature profile. As a result, the higher the temperature, the higher the arc root

C. Influence of the air mass flow rate

The influence of the air flow rate has been investigated in the range of 0.16 -0.5 g/s for a current of 400 mA.

Fig. 18 shows that the increase of the air flow rate results in the increase of voltage instabilities. For an air flow rate of 0.5 g/s, the arc shape is very disturbed by the flow instabilities appearing in the post-discharge zone. An example is shown on Fig. 19. The simulated oscillations of the jet implies a periodic and significant arc voltage increase in order to maintain the arc. The voltage shapes remind the three operation modes [START_REF] Duan | Arc instabilities in a plasma spray torch[END_REF]: steady, takeover and restrike for 0.16, 0.32 and 0.5 g/s respectively. Even if those are simulated results, two candidate mechanisms for such instabilities should need further studies to confirm or not their role in this case: the instability of a jet in a cavity, and the Kelvin-Helmholtz instability.

The first process involves the instability of the pressure field surrounding the jet near the inlet. This is a classical version of the Coanda instability, and has been reported even in the case of thermal jets [START_REF] Muller | The role of inlet location on the dynamics of thin two-dimensional buoyant surface jets[END_REF]. In our case, it results in an oscillation of the jet, but observational evidence reports that in similar cases, wall-hugging jets can be also observed.

Table III shows that, as observed experimentally, the higher the air flow rate, the higher the voltage drop and the higher the instabilities. We can also observe that the higher the mass flow rate, the higher the discrepancy between simulated and numerical mean voltage. Indeed, the higher the mass flow rate, the higher the cooling of the arc by convection heat transfer, and the higher the non-equilibrium effects.

The radial profiles of the temperature and velocity in the The lower the mass flow rate, the wider the temperature profile. We can notice that the maximum temperature reached in the center of the arc core for the three mass flow rates are very close (between 5140 and 5230 K). Hence, the arc core temperature depends mainly of the current and the mass flow tends to cool down the arc core outline. The maximum velocity is equal to 30 m/s for Q air = 0.16 g/s and 62 m/s for Q air = 0.5 g/s. At low flow rates, the velocity profile is relevant of the pure laminar parabolic case. It can also been seen that at higher air flow rates, deviations to this reference profile are retrieved, even if no turbulence model is activated in the code. At the corresponding regimes (Re = 2000, typically) transitional mechanisms toward turbulence are strictly impossible to inhibit, so this distortion of the velocity profile can be linked with several candidate processes, from purely numerical effects to physical mechanisms, such as the Kelvin-Helmholtz instability that generally develops in the shear-stress zone surrounding jets. In our case, the vortex creation that should progressively occurs at scales between the grid size and the jet diameter, resulting in an increase of the mean velocity in this mixing layer is not observable, but there remains the evidence of a macro-scale shear-stress effect. This has already been described by Trelles et al. [START_REF] Trelles | Non-equilibrium modelling of arc plasma torches[END_REF].

IV. CONCLUSION

The MHD unsteady modeling of a non-thermal DC nontransferred low current -high voltage plasma torch operating with air has been successfully implemented. The low current arc discharges are characterized by a highly irregular and unstable motion of the arc column which is particularly difficult to simulate due to numerical instabilities. Significant information, experimentally hardly accessible, have been obtained.

Among the most important results to point out, one can quote:

• The numerical arc moves linearly along the same plane and has the same shape as the experimental one. Due to a velocity dead zone at the cathode tip, the arc root is attached few mm above the cathode tip as observed experimentally.

• The magnetic field and the Lorentz force are very low because of the low current applied to the anode. The arc column is very stable in the nozzle because of the wall-stabilization of the arc. At the nozzle tip, the low current leads to instabilities in the arc length and thus in the voltage which could be linked to pure hydrodynamic effects.

• The model overestimates the voltage drop probably because of the LTE assumption. Indeed, the non-equilibrium effects are significant and result in a discrepancy between experimental and numerical voltage values.

• The temperature gradients are significant and lead to a highly non-homogeneous arc discharge.

• The arc root radiation is underestimated because the electrode metal evaporation is not taken into account.

• The arc core radius can be estimated between 0.5 and 0.83 mm by the model.

• For t > 10 ms, the numerical oscillation frequency range, similar to the experimental one, could probably be due to the hydrodynamic forces.

• The higher the current, the lower the voltage, the higher the deposited power, the higher the arc core temperature and radius, the lower the voltage instabilities.

• The higher the current, the lower the voltage discrepancy between the experimental and numerical results. The higher the air flow rate, the higher the voltage difference between experimental and numerical results. These two phenomenon are probably due to the non-equilibrium effects which decrease with increasing the current and decreasing the mass flow rate.

• The higher the air flow rate, the higher the voltage instabilities. Maximum temperatures are similar but the temperature and velocity profiles become narrower.

• For Q air ≥ 0.5 g/s, the jet appears progressively as instable and broadens. Such behavior for the model indicate that experimental confirmation for possible candidate mechanisms may be searched for : Coanda instability and Kelvin-Helmoltz instability.

In a further work, the development of a NLTE model will be considered to take into the high non-equilibrium effects of the plasma torch. The perspective of this preliminary model are also the implementation of a enthalpy flux and a sheath zone at the electrodes and a better description of the inlet zone. The full anode and inlet zone with tangential entry will be meshed in order to take into account the real vortex flow. The laminar flow regime is a strong assumption in our case and turbulent models will be tested to observe their influence on the arc shape and the voltage drop.

The restrike mode will also be implemented to simulate the glidarc behavior, happening for certain values of current and air flow rate. A hot gas column reattachment process, based on the overcome of a fixed value of the electric field, is considered.

This preliminary work is likely to open the way towards a better understanding of low current arc discharges which technologies are currently encountering an important development in many application fields.

Valerie Boch was born in 1963. She received the Master dedree in aerodynamics from University of Poitiers, France. Since 1987, she has been with Renault SAS, France, where she is currently head of advanced fuel team in the R&D department. As a Senior Scientist in the Center for Energy and Process in the same school, he has developed since 1992 several advanced methods in the domain of radar remote sensing for Geophysical flows studies in Meteorology and Oceanography. He is specialised in visualisation techniques and modelling of geophysical and industrial flows, with a special emphasis on decision support systems regarding risk and environmental problems, and has published more than 40 papers in scientific journals and conferences. He has leaded task forces in several international projects for UNEP, IMO and EU, and has participated to the design of the first Numerical Information Systems integrating spatial oceanographic observations for environmental plans. He is working today on adapting the numerical modelling tools of Fluid Mechanics for project management integrating Renewable Energies. He joined Mines-ParisTech (formally Ecole des Mines de Paris) in 1989 as a Senior Researcher, became head of the Plasma Group in 2003 and Maitre de Recherche in 2007. His field of research covers the development of thermal and non-thermal plasma processes for energy and environmental issues. These researches are currently focussing on fossil and renewable hydrocarbon conversion including : reforming, gasification, Gas To Liquid conversion, gas treatment (NOx, COV, particles. . . ) and nanostructures synthesis (carbon black, fullerenes, nanotubes).

Laurent Fulcheri was born in

Laurent Fulcheri has more than thirty papers published in peer review journals, and around ten patents. He participated in many international conferences as invited speaker.

Fig. 1 .

 1 Fig. 1. Schematic of the plasma torch.

Fig. 2 .

 2 Fig. 2. Computational grid of the plasma torch.

Fig. 3 .

 3 Fig. 3. View of the initial hot temperature channel at the first time step.

Fig. 5 .

 5 Fig. 5. Representation of the velocity magnitude distribution. I = 400 mA. Q air = 0.32 g/s.

Fig. 6 .Fig. 7 .

 67 Fig.6. Numerical temperature profile in the radial direction taken in the middle of the cathode. I = 400 mA. Q air = 0.32 g/s.

Fig. 8 .

 8 Fig. 8. Representation of the radiation emission distribution. I = 400 mA. Q air = 0.32 g/s.

Fig. 9 .

 9 Fig. 9. View of the current density and isopotentials. I = 400 mA. Q air = 0.32 g/s.

Fig. 10 .

 10 Fig. 10. Simulated arc voltage versus time. I = 400 mA. Q air = 0.32 g/s.

Fig. 11 .

 11 Fig. 11. Typical experimental voltage oscillogram. Q air = 0.32 g/s. I = 400 mA.

Fig. 13 .

 13 Fig.[START_REF] Zhou | Comparative study of turbulence models on highly constricted plasma cutting arc[END_REF]. Photographs of the experimental plasma for three currents. Q air = 0.32 g/s.

Fig. 14 .

 14 Fig. 14. Numerical arc voltage vs. time for three different currents. Q air = 0.32 g/s.

Fig. 15 .Fig. 16 .

 1516 Fig. 15. Experimental oscillograms for I = 300 and 600 mA. Q air = 0.32 g/s.

Fig. 17 .

 17 Fig. 17. Radial temperature profile at z = 35 mm (middle of the cathode) for three different currents. Q air = 0.32 g/s.

Fig. 18 .

 18 Fig. 18. Numerical arc voltage versus time for three different air flow rates. I = 400 mA.

Fig. 19 .

 19 Fig. 19. Temperature distribution and velocity vectors field illustrating the instabilities. Q air = 0.5 g/s. I = 400 mA. t = 20 ms.

Fig. 20 .

 20 Fig. 20. Numerical radial temperature profile at z = 35 mm for different air flow rates. I = 400 mA.

Fig. 21 .

 21 Fig. 21. Numerical radial velocity profile at z = 35 mm for different air flow rates. I = 400 mA.

Vandad

  Rohani was born in Paris, France, in 1980. He received the Ph.D. degree in plasma physics from University of Paris-Sud XI, France, before joining, in 2008, the Center for Energy and Processes of MINES ParisTech, France, as a researcher in plasma physics and part-time lecturer. His research activities are focused on the development of thermal and nonthermal plasma sources for energy applications and chemistry. Francois Cauneau obtained a DEA in Astrophysics in 1985 at Observatoire de Paris, and a Ph.D. in 1988 in Remote Sensing at University Paris 7. He is Maitre de Recherches at Ecole des Mines de Paris -Mines ParisTech, where he teaches in fluid dynamics, wind turbines aerodynamics and engines technology.

  France in 1960. Graduated from Ecole Nationale Superieure d'Arts et Metiers (ENSAM-Paris) in 1983, he got a DEA from INP-Grenoble in 1984, a Ph.D. and Habilitation degree in 1989 and 2003, respectively.

TABLE I BOUNDARY

 I CONDITIONS OF THE 3D MHD MODEL.

IEEE TRANSACTIONS ON PLASMA SCIENCE, VOL. 39, NO. 9, SEPTEMBER 2011 © 2011 IEEE Digital Object Identifier 10.1109/TPS.2011.2160208

ACKNOWLEDGMENT

The authors gratefully acknowledge the Code Saturne development team for their technical support.