Olivier Maurice
email: olivier.maurice@airbusafran-launchers.com

Philippe Durand
email: philippe.durand@cnam.fr

Modelling cobotic systems

Keywords:

come

T he method of Kron allows constructing the equations of any coupled system. It starts from the definitions of all the components of the system. We consider here the system as being a robot. Using a serie of transformations, it leads to the conception of the whole robot. The system of equations is a system of coupled equations multiphysic. This system of equations can be modified to become a manifold which is a mathematical image of the robot. The hyper-surface associated with the manifold can be transformed thanks to a matrix that includes probabilities coming from a game. This game occurs by involving the robot and a human. Some observables of the system serve as indicators for the actions of the robot. These actions are seen like decisions taken by a gamer [START_REF] Fudenberg | Game theory[END_REF]. By postponing these decisions in a payoff matrix of a game theory, we study the interaction between the robot and a human. From these interactions, the human takes some probabilistic decisions which influence the expression of the manifold. We first describe the tensorial analysis of networks applied to the robot definition. Then we show how the system of equations deduced from the first task can be seen in differential geometry. This leads to the manifold definition. We give then the technique modifying the manifold by probabilities associated with the human decisions. To end we show how the robot can influence human decision inside of the game theory, coupled with TAN for modelling of the human machine interaction. To conclude we speak of various project to be done in the future.

Primitive parts

Robots are made of various parts (we supposed here that these parts are solid parts). Before to be assembled, each part can be studied independently and characterized. We speak here of parts models. The modelling the parts comes from measurements applied on these parts. The experiment consists in submitting the part to an effort, then to measure the flux created by this effort. But this flux can be measured in various conditions of environment. So, the law obtained that gives the relation between the effort w and the flux f is a summation of many functions giving w = z(f ) for each condition of environment [START_REF] Maurice | Proposition d'un formalisme comme support pour les études théoriques en systémique[END_REF]. If T is a variable of the environment, we make the assumption that it exists a function D defined by: n(T ) D(n) = 1. For each value of n that corresponds to one interval of values associated with T , a law w = z(f ) is defined. When the environmental conditions change, the law is changed also, continuously from one interval to the other. In general, the laws have the form: w = D 1 .z 1 (f ) + D 2 .z 2 (f ) + . . . + D n .z n (f ). The logistic function can be used for F (with a, b ∈ R):

T ∈ [a, b] ⇒ [a,b] D T = 1 1 + ce -α(T -a) - 1 1 + ce -α(T -b) (1)
For the parameter T that belongs to the interval [a, b], the law associated with the part is given by:

[a,b] D T z 1 (f ) = w. For a next interval of values [b, c] the law becomes [b,c] D T z 2 (f ) = w.
The global behaviour of the part is so given by:

w = [a,b] D T z 1 (f ) + [b,c] D T z 2 (f )
With this decomposition in domains [START_REF] Maurice | PROPOSITION D'UNE THEORIE POUR LA CEM[END_REF], the set of laws describes the whole behaviour of the part for any condition of the environment. For example taking into account of the temperature T, of the pressure p and of the altitude A the law can be expressed by (where i n are intervals in R):

w = i 1 D T i 1 D p i 1 D A z T 1p1A1 (f ) + i 1 D T i 1 D p i 2 D A z T 1p1A2 (f ) + . . .
which can be generalized following:

w = k Π n in D p k z in (f ) (2)
where p k represents the various parameters.

Representation of the primitive part by a branch in a graph

When the effort is applied on the primitive part, it results the development of a potential associated with the equivalent work [START_REF] Maurice | Analyse tensorielle des réseaux hybridée comme outil théorique pour la CEM et l'électromagnétisme[END_REF]. A classical representation of this process is Kirchhoff's branch. In electromagnetism, the effort (w) becomes the electromotive force, the flux is the current (f ) and the potential (U ) is the potential difference in electricity. The relation between these three quantities is given by Kirchhoff's law:

w, f, V ∈ R, w = z(f ) + U (3)
The sum of products involving the flux in (2) can be seen as an operator ζ applied to f (ζ • f ):

{i 1 ,i 2 ,...} Π p k in D p k z in (f ) ⇒ ζ ≡ {i 1 ,i 2 ,...} Π p k in D p k z in ( ) (4)
so that the relation (3) can be written:

w = ζ • f + U
To make a robot we need several primitive parts [START_REF] Maurice | Proposition d'un formalisme comme support pour les études théoriques en systémique[END_REF]. To each part, an equation of Kirchhoff is created after some physical characterizations. It means that each real object is now represented by an equation with many parameters and one variable. We consider the equations as in bijective correspondances to the real parts of the robot. We obtain finally a set of N equations to cover the whole parts of the robot. This collection is the first step that leads to the robot modelling.

The k parts are associated with k equations of the form:

w k = k ζ kk • f k + U k ( 5
)
Each equation is defined in relation with a branch in a graph. We name C o the collection of equations and T 1 the set of the branches in a graph G. A correspondence is created between C o and T 1 , associating each element of C o to each elements of T 1 . In this approach, T 2 represents the cycles in the graph G and T 0 the nodes.

Constructing a system

The first step to construct a system, is to construct its subsystems [START_REF] Journée | Simulation et Complexité, notions et méthodes[END_REF]. A subsystem is a set of branches connected. We accept the idea that any real object is an assembling of components. Each component has been characterized and associated with a physical law as described before. The components connected on the real object are represented by branches connected in the graph created to represent the object. This construction leads to a set of branches linked through an incidence relation. The incidence matrix gives the relation between the nodes and the branches in the graph, translating by these relations the first step in the construction of the subsystem. Then an operation allows to modify advantageously the equation in order to give it more riemannian characteristics. That's the last operation to do to obtain a tensorial equation representing the robot.

If A is the incidence matrix, B the vector made with the ordered set of the branches {B} ∈ T 1 and N the vector made with the ordered set of the nodes {N } ∈ T 0 , the relation between branches and nodes is given by:

B a = d A d a N d (6)
(6) gives all the information to understand how each subsystem is made. It has a mathematical image. Each branch being defined by its own physical law, grouping the branches in a common subsystem using an incidence matrix leads to realize a direct summation of the physical laws [START_REF]Différence fondamentale entre la description nodale et la description modale de la topologie électromagnétique dans un objectif de Diakoptic[END_REF]. If ζ (a)(a) is the physical law of one component (a) involved in the subsystem Su, {B} being the set of branches pointed out by the incidence matrix A, the global operator Θ associated with the subsystem Su is defined by (in the suite of the article we consider ζ, Θ, χ, Ω, . . . as operators):

Θ = ⊕ dim(B) n=1 ζ (n)(n) (7)
an identical operation being done for the fluxes, we obtain for the group of subsystems one equation given by:

ω α = α Θ αα .f α + U α
In this description misses the radiative interactions or, in general, all interactions not transmitted by solid materials [START_REF] Reineix | On mathematical definition of chords between networks[END_REF]. It means that an operator χ that takes in charge these kind of interactions must be added to Θ. χ is a non diagonal operator. It results from the induced fluxes that belong to one branch B (a) coming from efforts located on another branch B (b) , a = b. The entire subsystem Su is then described by the equation:

ω α = β (Θ αβ + χ αβ ) .f β + U α (8)
This can be resume using Ω = Θ + χ by:

ω α = β Ω αβ .f β + U α (9)

Connecting the subsystems

A system S is then made with N subsystems Su in a similar operation in two steps. In a first step we realize the connection of the branches in cycles through an adequate connectivity. In a second step we add interactions between cycles to the previous operator elaborated in the cycles space. Each subsystem is by this way transformed in an operator defined in the cycles space. The interactions between these connex graphs realized in the cycles space leads to the definition of the whole system.

First step

For any graph of definition (G, A), it exists a spanning tree that reaches all the nodes of the graph [START_REF] Reineix | Kron's method and cell complexes for magnetomotive and electromotive forces[END_REF]. If {N } is the ordered set of the nodes, in the order explored by the spanning tree T, the cycles are constructed by going back, each time it is possible, to the node considered in the tree T. The tree is described by the set of nodes {x 1 , x 2 , . . .} (the tree is first a set of branches that implies through A a set of nodes). Starting from the node x 1 , we travel through the tree until there is a possibility to reach x 1 using a branch that does not belong to the tree. This branch added to the branches passed through constitutes a cycle. Once there are no other possible cycle starting from node x 1 , the process is repeated starting from node x 2 . This mechanism is done until all the branches may be linked with cycles.

The number of cycles dim {M } , {M } ∈ T 2 is determined by Euler-Poincaré's law:

dim {M } = dim {B} -dim {N } + dim {R} (10)
{R} being the set of connex graph, or in other words, the number of subsystems Su in the system S.

As the cycles are determined using one chosen spanning tree, more than one possibility of set of cycles exists. We can wonder if some choice are better than others? Some empirical rules were established:

• the resulting functions associated with a cycle should not reach infinite values;

• when branches are boundaries between two subsystems, and when these branches are associated with filters, the retained cycles should share these branches;

• the cycles should show physically the energy flow with fluxes which start from efforts and go back to the same efforts.

Anyway, whatever a set of cycles retained, if the final operator obtained with this choice in the space of cycles is not well conditioned, another set of cycles may be searched.

Once the cycles are defined, the connectivity Λ giving the relations between branches and cycles is also defined, all the more if the same directions are chosen for branches and cycles. If f k are the fluxes in the branches space and J β the fluxes in the cycles space, we have the relation:

f k = β Λ kβ J β (11)
by replacing [START_REF]Différence fondamentale entre la description nodale et la description modale de la topologie électromagnétique dans un objectif de Diakoptic[END_REF] in ( 9) we obtain:

ω α = β Ω αβ . η Λ βη J η + U α (12)
The relation between f and J can be translated to the relation between w and e, e being the effort at the cycle level:

e κ = α Λ κα ω α
multiplying both members of ( 12) by Λ κα gives:

α Λ κα ω α = α Λ κα β Ω αβ . η Λ βη J η + α Λ κα U α (13)
but α Λ κα U α = 0 and if we make the assumption that:

β Ω αβ . η Λ βη J η = βη Ω αβ Λ βη .J η
we finally obtain:

e κ = αβη Λ κα Ω αβ Λ βη .J η (14)
by noting Λ κα Ω αβ Λ βη = h κη , ( 14) becomes:

e κ = η h κη .J η (15)

Second step

Some interactions can exist between cycles to construct the complete subsystems. To do this, we add an extra-diagonal matrix µ to h that includes the interactions [START_REF]De l'expérimentation sur des systèmes complexes en compatibilité électromagnétique, à leurs représentations et leurs analyses dans un espace géométrique abstrait[END_REF]:

e κ = η (h κη + µ κη ) .J η
Noting H κη = h κη + µ κη we obtain the equations of the subsystems given by:

e κ = η H κη .J η (16)
Now we have a set of subsystems (e κ , H κη ). To construct the whole system, we need coupling using cords or by sharing some branches. We look at these two techniques. First of all, the system is made using as previously the direct summation of the subsystems. This is made for the functions e and for the operators H:

e κ = ⊕ κ e κ H κη = ⊕ (κ, η) H κη
If cords have to be added, a matrix of these links Φ is added to H .

H κη = H κη + Φ κη
Other possibilities exist. One consists in sharing a common frontier which is a branch and its two nodes between the two subsystems. In the next explanation, the branch b ∈ T 1 wears the operator b.

We consider a first subsystem Su 1 with a frontier

[(n 1 , n 2 ) , b] n 1 , n 2 ∈ {n i , n i-1 , . . . , n 2 , n 1 } ⊂ T 0 , b ∈ {b i , b i-1 , .
. . , b} ⊂ T 1 and a second one Su 2 with the same frontier. When we make the direct summation of the two subsystems, it means that we construct an operator H including both operators H 1 and H 2 of the two subsystems. This leads to a matrix that has the form:

H =     . . . . . . . . . z + b 0 0 0 0 0 0 0 0 b + t . . . . . . . . .    
z and t being any operators included in H 1 and H 2 . b is the shared branch. b is the shared branch. So, when connecting the two subsystems that share the branch b and its nodes (it means that b and its nodes belongs to the common sets of branches and nodes of the unified system), the branch b becomes the interaction operator and the matrix H takes the form:

H =     . . . . . . . . . z + b 0 0 -b 0 0 -b 0 0 b + t . . . . . . . . .    
In that case, this is equivalent to say that the matrix of interactions Φ is composed of the shared branches operators (here for dimension 2):

Φ = 0 -b -b 0
Finally, the system is completely described by the equation:

e κ = η H κη .J η (17)

Standardised system definition

The elements of the operator H contain the functions defined on domains (even cords can be defined on domains). The intervals associated with the domain definitions give all the informations on the limit values of the parameters. The system is at this stage completely defined by the equation and the intervals. But it remains difficult to analyse the system because the matrix H can be non symmetric [START_REF] Leman | Use of the circuit approach to solve large EMC problems[END_REF].

For example, if we consider a simple problem of automatic: a feedback loop of direct gain G. The graph of this kind of system is made of three cycles. Noting x k the sources and y m the flux, a first cycle of law H 11 = 1 gives y 1 = x 1 , without any other feedback. A source is reported by a cord on the second cycle: x 2 = Gy 1 . Choosing a law H 22 = 1 for the second cycle also, this implies that: y 2 = Gy 1 . Now for a feedback loop of gain -R, we report on the first cycle a source x 1 = Ry 2 . The output of the system can be extracted with a third cycle (H 33 = 1) giving the result in y 3 , and a cord that links y 2 with y 3 through a function -h. Finally:

   x 1 = y 1 -Ry 2 0 = Gy 1 + y 2 0 = -hy 2 + y 3 (18)
Compared to [START_REF] Maurice | Electromagnetic game modeling through Tensor Analysis of Networks and Game Theory[END_REF] with e ↔ x, y ↔ J, it leads to:

H =   1 -R 0 G 1 0 0 -h 1   which is not symmetric.
But we can consider e = H(J) where the e are the local coordinates of a tangent plan T pS to the hypersurface S described by the curvilinear coordinates J. To do that, we make a local parametrization of S with:

1. a domain I, the parametrization domain;

2. a map ψ : I → R n of class C 1 and which verifies ψ (I) ⊂ S.

Once these definitions are given, we can calculate the base vectors b k associated with the local tangent plan T pS:

b k = ∂ψ ∂J k (J k ) (19)
In the previous example we should obtain:

b 1 =   1 G 0   b 2 =   -R 1 -h   b 3 =   0 0 1  
From these basic vectors we can define a metric G:

G km = b k , b m (20)
Now in the idea to obtain some similar approach as Riemann's approach, we would like to write something like (by using now Einstein's convention in the rest of the document):

T α = G αβ .J β
T being a tensor for inertia and sources/efforts. In [START_REF] Maurice | Electromagnetic game modeling through Tensor Analysis of Networks and Game Theory[END_REF] when there are time derivations of fluxes, their derivation versus the fluxes leads to zero. For the electric for example, we have the electromotive forces L∂ t i (L is the inductance, i the electrical current). Making ∂ i L∂ t i gives 0. So, to retrieve ( 17) using (20) we need to add some terms each time the derivation versus the fluxes or any other operation doesn't allow to obtain [START_REF] Maurice | Electromagnetic game modeling through Tensor Analysis of Networks and Game Theory[END_REF] with only (20).

Remarking that ∂e κ ∂Jη gives the jacobien matrix J η κ , we can first write (17) under the form:

e κ = η J η κ .J η + N η κ ∂J η ∂t + . . .
under the assumption that this new form gives back equation ( 17) thanks to added terms like N . We can multiply each member by the transposed matrix of J : Γ, which leads to:

κ Γ κ σ e κ = κ Γ κ σ η J η κ .J η + N η κ ∂J η ∂t + . . .
By making a direct equivalence between the fluxes J η and their contravariant correspondences J η and by noting:

T σ = κ Γ κ σ e κ - κ η Γ κ σ N η κ ∂J η ∂t + . . . knowing that: G ση ≡ κ Γ κ σ J η κ ( 17
) can now be written:

T σ = G ση .J η (21)
This equation, completely identical to 16, presents the advantage that this time, the metric G is symmetric and it gives to Kron's equation more of characteristics of a riemannian space, including eigenvalues. But we must remember that the expression of G is many more complicated than in usual riemannian cases. The relation (4) leads to basic vectors of the same form and to fluxes defined on intervals. By the fact, equation (21) can be seen as the one of a manifold M [START_REF] Maurice | Seconde Géometrisation pour la CEM[END_REF].

At this step we can accept that the complete physical behavior of any robot can be modeled using (21) associated with a set of intervalls given for each variables:

M ≡    T σ = G ση .J η i 1 (p 1 ) = [v 1 , v 2 ] . . . ( 22
)
This manifold is also represented by its associated graph (G, A).

Basics for game theory

The other part of the formalism uses game theory. The game here is between the robot and some living being. Major principles of game theory are incorporated in a tabular which leads to the matrix Θ which includes the payoff matrix D and that bayesian -markovian matrix P . In the robot we identify some particular outputs that correspond to the robot's actions. On another side, some particular sollicitations of the robot sensors imply specific reactions taken by the robot. The reactions of the robot are translated through particular values of fluxes acting on the outputs. The set of fluxes for outputs is noted {J r }. The robot is one of the gamer. But to give to its actions some symbols more understandable, we use a neuronal network Ξ which takes in charge the correspondences between the set of values for outputs {J r } and the action symbols C k . These symbols are seen as some choices for the robot:

C k = o Ξ k o J o , J o ∈ {J r } (23)
The property of the neuronal network [START_REF] Howard | Neural network design[END_REF] is to transform the continuous values of J o into discrete choices for C k . Note that J belongs to a natural space as C which comes from J through the matrix Ξ. It can be considered as a part of the robot itself, like an added layer for the outputs. In front of the robot reactions, we have the reactions of the living being, for an assumption of psychological profile. These reactions are identified in a set of possible discrete actions {π h }. This set makes a dual space to the choices of the robot as complex behavior of the living being cannot be directly transcripted. It can be projected in a vectorial space having directions of psychological behaviors, but needs to be synthesized in real values to identify each family of reactions. The purpose of the mixed probability and payoff matrix Θ ασ is to give, for each couple of choices (C σ , π α ) the gain that the living being should win. For each choice of the living being, we can compute the average hope of earnings depending on the choices of the robot and of their probabilities. But in fact the matrix Θ ασ computes the partial hopes of earnings for the living being [START_REF] Maurice | Electromagnetic game modeling through Tensor Analysis of Networks and Game Theory[END_REF]. If ε is the average hope of earnings:

ε α (π α ) = Θ ασ (P (σ, α) (C σ | π α,t -, E)) C σ (24)
where P (σ, α) (C σ | π σ,t -, E) is the bayesian probability for the robot playing action C σ knowing that the living being has previously played π α,t -and under environmental conditions E. Under some rational assumptions, the living being should make the choice of the highest hope of earnings. And this choice should imply some reaction of the robot, feeling this decision taken by the living being (through its sensors). So, some given choices from the living beings must be translated in signals (commands or signals perceived by the robot's sensors) for the robot. To do that, we use a sophisticated neuronal network ℵ α κ that will translates the hope of earnings ε α in signals for the robot:

e i = ℵ α i P (α, σ) (π α | C σ,t -, E ε α , e i ∈ {e κ } , e κ → {e κ } (25)
{e κ } being the set of components for commands included in the covector e κ , and P (α, σ) (π α | C σ,t -, E as previously the bayesian probability that the living being plays π α depending on the action made by the robot on past times and depending on the environmental conditions.

Coupling the robot and the living being

The psychological profile of the living being is principally translated through the object ℵ. ℵ is a sophisticated matrix which gives to each hope of earnings, a signal in direction of the robot. This signal is linked with some reaction of the living being which is perceived by the robot's sensors. But it can be also a direct command that a person applies to the robot. ℵ depends on the bayesian probability P σ (π α | C σ,t -, E . It means that in the neuronal network ℵ the probabilities are the weighting coefficients on one of its layers. These weightings are partly established by the past behaviour of the robot.

In each case of the probabilities P or P , the past means to go back more or less far in the previous actions depending on the system complexity and depending on its capacity to memorize the past. Both matrices Ξ et ℵ can contain complicated functions to transform the fluxes into actions and the decisions into signals for the robot. The whole system including the material model, the two neuronal network and the matrix of game covers both psychological and material parts of the coupled system.

In both cases of Ξ and ℵ, logistic functions allow to transform continuous values into discrete values. By using thresholds, a weak probability attached to a strong hope of earnings can lead to zero if the threshold for a signal is not reached. This gives a model to understand why the best cooperation is not always obtained between living beings and machines. This comes from the fact that the reactions of the living beings may answer to individual gains rather than to the interest of the whole system.

Model of the living being

In case of a person having a rational behaviour, the fact to take the hope of earning maximum as a criterion to determine his choice is an efficient method (this process is similar to put to one or zero the probability coefficients). But in a more general case, the relation between the vector of hope of earnings and the signal perceived by the robot is many more complex. Various assumptions can be studied including suicide behaviours etc. and knowing that the environment can influence the decisions. So, for each family of psychological profiles, interactions between the robot and the living being will be studied. We consider living beings in general, because all living beings have complex strategies in response to various situations. When taking into account the psychological profiles that represent some probabilities in the whole population, several cases are studied in parallel with in each case different interactions between the robot and the living being. The bayesian probability associated with the robot decision allows to take into account the uncertainties of its sensors, the influence of the environment, the influence of other possible gamers, etc., and, of course, of the previous decision of the living being.

In this model we consider that some quite complicated algorithm are implemented in the robot that justifies the use of bayesian probabilities to model its reactions. Anyway, everyone will be able to adapt the two layers Ξ et ℵ depending on the model he wants to develop. Note that the payoff matrix Θ is already difficult to be constructed. The gain can depend on social considerations, on humanist considerations or simply on technical considerations.

In cobotic we often consider the environment of firms or private residences or agricultural environments. In these cases the higher hope of earnings is a good criterion all the more if it is in relation with the wanted objectives considered as standards. Average and standard deviation in comparisons between computations and experiments may give the same statistics under these assumptions. For animals, we can think that the conditions are ob-tained if the animals are living near to their natural environment.

Illustration

We consider a very simple robot. It's a mobile system that can go ahead or stop. The mobile can get some information from the environment with a sensor that is sensitive to the movement. We begin to model the robot. In this simple case, it is no use to develop any neural network. The input and output of the mobile are directly associate through cords to the game theory mixed matrix Θ. The robot working can be represented by the next graph. We can see a first bloc associated with the robot perception. This sensor can command an electronic in charge of the machine movement. A last bloc is the machine. Figure 1 shows the whole system we consider. For the moment we only focused on the high part of the graph, without looking at the matrix in the low part. First cycle numbered 1 is associated with the perception function.For a signal received on its input, the electronic gives an answer saying if the input signal has changed between two time samples. This sensor says if something has moved in the scene in front of the mobile. We consider that the pedestrian in front of the mobile is static, except if he wants to cross the road. If he crosses the road, the output of the first cycle sends a signal to the second cycle. The transmission is taken in charge by a cord which reports an electromotive force (emf) on the second cycle in relation with the current in the first cycle. This emf comes to suppress the command of the machine existing on cycle 2. The output of the second cycle comes to command the third one which is the machine circuit. We can gives each cycle an arbitrary impedance operator equal to one. Without the coupling function (cords), the impedance matrix is an identity matrix of dimension three. To this matrix we add the interaction operators t ij that link cycles 2 and 1, 3 and 2. This leads to the matrix ζ:

ζ =   1 0 0 t 21 1 0 0 t 32 1   (26)
The equation e α = ζ ασ f σ gives the solution through the fluxes f σ . Depending on the flux values, high or low, the correspondence to know what the robot will do consists in making a link between the flux value and the robot action: to go ahead or stop. These two choices can be seen as the choices of a gamer: the robot. The relation between the machine and the pedestrian can be seen as a game. The second gamer is the pedestrian while the first one is the robot. The robot actions a comes from output fluxes f . As we don't use a neural network in this simple illustration, a simple connection matrix Ξ selects the fluxes f used for the actions a: a σ = Ξ σ σ f σ . We can define the probability for the pedestrian to go ahead knowing that the robot stops. Four probabilities must be define by the same way. We write in general:

P (α, σ) = P (r α | a σ ) (27)
r j being the action of the pedestrian. The output of the game are the partial hopes of earnings for each choices of the pedestrian. This means that we can associate to each combination of choices a value of earning. We note Θ ασ the mixed payoff and bayesian matrix. The partial hope of earnings α are given by:

α = Θ ασ a σ (28)
The α being translated in sources e α for the cycles through the matrix ℵ α α : e α = ℵ α α α , because we don't use neural networks in this simple illustration. Θ ασ = D ασ P (α, σ), D being the gains that the living being can hope for each of his actions r in reponse to the robot actions a and P (α, σ) the probability of each couple of choices.

Finally, the couple of equations that describes the system is given by :

   e α = ζ ασ f σ e α = ℵ α α Θ ασ Ξ σ σ f σ (29)
These two equations can be reduced in an unique equation. In fact, the game part is a special and complex cord but leads to the interaction between the output f 3 of the robot and its command e 2 . The global object ℵ α α Θ ασ Ξ σ σ takes place as the coupling operator ζ 23 , as shown in the graph. In our case, e α is reduced to e 2 . So:

ζ =   1 0 0 t 21 1 ℵ α 2 Θ ασ Ξ σ 3 0 t 32 1   (30)
The payoff matrix can be quite easily established, once some criteria defined. It's more difficult to determine the probabilities. Measurements can be done, or hypothesis, based on previous experiences and psychological arguments. What is clear is that the probability law is a function that depends on a σ and r α . The dependance on r α can be replaced by a dependance on β αα f α where f α in this case is a known flux. This allows to compute after the theoretical analysis of the system by derivating the function e α (f σ ).

We write

ℵ α 2 Θ ασ Ξ σ 3 (f ν ) = η (u [f ν ]
). We can compute the Jacobian J of the system e α = ζ ασ f σ . It gives for each of its columns the base vectors of the parametrized hypersurface e α :

J ν α = ∂ f ν e α =       1 0 0 t 21 1 u (f ν ) η (u) 0 t 32 1       (31)
From the jacobian we can deduce the metric G:

G = J T J =       1 + (t 21 ) 2 t 21 t 21 u η t 21 1 + (t 32 ) 2 u η + t 32 t 21 u η u η + t 32 1 + (u η ) 2       (32)
The source becomes:

T = J T e = E p t 21 E p + E c + u η E m t 32 E c + E m (33)
E x are the source for each mesh. The new equation that represents the problem is now:

T α = G ασ f σ (34)
This is the equation that we study.

Before to go further we need to define the bayesian probabilities and to detail θ mechanism.

If w is the action of waiting for, for the pedestrian, and T the action of going ahead and if S is the action of stop for the machine and A the action of rolling, we can define the probabilities following:

                   P (w|S) = e -ad ωA dtf 3 (t) P (w|A) = 1 -e -ad ωA dtf 3 (t) P (T |A) = e -βvm , v m = f 3 χv 0 P (T |S) = 1 -e -βvm (35)
d ωA is the distance between the car and the pedestrian, v m the car speed. From these equations we can generate the law giving the source reported on mesh 2. The source is given by:

           e ω = γ 1 0 D f 3 θ 00 P (w|S) + 1 D f 3 θ 01 P (w|A) e T = γ 2 0 D f 3 θ 10 P (T |S) + 1 D f 3 θ 11 P (T |A) (36)
γ i are factors for adjusting the values to the request command for the machine. The source reported on mesh 2 is the higher of the two values e ω et e T . We see in that case that the relation linked with the cord ζ 32 is a complex operator involving probabilities and game theory. We have:

u η = ∂ f 3 max e ω (f 3 ) f 3 , e T (f 3 ) f 3 (37)
We compute:

∂e ω ∂f 3 = ėω = 1 D f 3 θ 01 - 0 D f 3 θ 00 γ 1 ad ωA f 3 e -ad ωA dtf 3 (t) (38) and ėT
= 0 D f 3 θ 10 + 1 D f 3 θ 11 βγ 2 χv 0 e -βvm (39)
knowing both e ω and e T we can compute the component u η :

u η = max ėω f 3 -e ω [f 3 ] 2 , ėT f 3 -e T [f 3 ] 2 (40)

Analysis

We can verify how the machine command changes depending on the output of the machine by calculating:

∂T 2 ∂f 3 = ∂ ∂f 3 [G 2γ .f γ ] (41)
We have, as no flux depends on other fluxes:

∂T 2 ∂f 3 = ∂ ∂f 3 u η .f 3 + t 32 (42)
Now, f 3 being a machine command, so a power one, we can make the assumption that f 3 >> 1. In that case,

It appears clearly that e ω > e T , if θ ij > 0, ∀i, j. It means that the pedestrian has has all interest in waiting that the vehicle passes far. The higher the distance, the higher is gain in waiting. This shows that our law for the bayesian probability is logical but not human! But our naive example shows also how the formalism can lead to rigorous and mathematical conclusions, once the good laws for the behaviors are found. It shows also how the technique works to include game theory inside the system of analytical equations that are the lagrangian expression of the problem. Various cases should be improved to test the formalism even if the psychological modeling is not relevant. It will give to the method robustness and reliability, allowing after a deeper work focused on the psychological side. This psychological side is taken in charge by the interactions (or cords) using game theory and neural networks.

Conclusion

The xTAN method and formalism presented seems for us to be a good candidate for modeling complex systems. It can include all physics making links between characterizations under various environments using domains. It can include multiscale problems through cords that take in charge operators of integration that make the change of scales. Using xTAN for the cobotic seems under these facts like a natural conclusion. In the society that begun to be constructed, the cobotic will be a major dimension in cyber-physical systems. In order to analyze and study the behaviors of this complex systems and their relation with humans and animals, we need a mathematical tool that allow to model them. xTAN can be this tool or, at least, merit to be test for this objective.

Future projects will consist in testing xTAN in various situations for the cobotic in cyber-physical systems and societies. More, taking into account the electromagnetic environment, we will be able to study the risk of disturbing these systems in front of different menaces, may them be volontaries or not. We can think in of electromagnetism conflict inside cities where many electromagnetic sources exist. In time of e-cities, this kind of problem will be of major interest for all societies.

Figure 1 :

 1 Figure 1: Complete system

∂T 2 ∂f 3 ≈ 1 D f 3 θ 01 - 0 D f 3 θ- 0 D f 3 θ 10 + 1 D f 3 θ 11 γ 2 β 2 χ 2 v 2 0

 3103012 max (ë ω , ëT ) 00 γ 1 ad ωA (1 -(f 3 ) 2 ) e -ad ωA dtf 3 (t) ëT = e -βvm


Anne Jaffrézic

Annie Soulier

Lucie Carrera

Barbara Le Bot

Impact of agricultural practices on veterinary pharmaceutical occurrence in superficial waters

Impact of agricultural practices on veterinary pharmaceutical occurrence in superficial waters

Jaffrézic A., Soulier A. ., Carrera L., Le Bot B.,., Jardé E Background and objectives Diffuse contamination of surface water by veterinary pharmaceuticals (VP) has been poorly studied in intensive breeding context. Recycling animal waste on soil is supported to fertilize crops and reduce utilization of mineral fertilizer but may have some environmental impact by recycling pollutants on soils and transferring pollutants from soils to rivers. Most of pharmaceuticals are consumed indifferently by human and animal, especially antibiotics. This study aims to quantify veterinary pharmaceuticals (VP) in superficial water in two agricultural catchments. One is dedicated to the production of drinkable water. Agriculture but also waste water treatment plant may contribute to the contamination of water and it is difficult to distinguish animal versus human source of contamination.

Methods and results * Grab and storm events water sampling were carried out from March 2013 to December 2014 on the 4,9 km² Kervidy-Naizin headwater catchment located in Brittany (western France) belonging to the AgrHyS environmental research observatories and part of Network in Drainage Basins. Water quality is impacted by intensive breeding (swine, cattle) with mean nitrate concentration of about 80 mg/l. Grab sampling were performed on different sub catchments from 2 to 80 km² located in Ille et Vilaine, impacted both by agricultural diffuse pollution and a sewage treatment plant. Priorization of molecules was established by a veterinary survey.

On Kervidy-Naizin,Thrimethoprim, Oxytetracyclin and Enrofloxacin are frequently quantified (from 57 to 42%). Mean and max concentrations respectively (12 to 80 ng/l) and (20 to 230 ng/l) are identical in flood and non-flood events. The frequency of quantification of Eprinomectin, a cattle veterinary antimicrobial, is 31% and mean concentration is high (400 ng/l). On the Ille et Vilaine catchment, human pharmaceuticals are frequently quantified Diclofénac (77%), et Carbamazepin (55%). But human/veterinary Flunixin (69 %), Lincomycin (56%) and Fluméquin (55%) or specifically veterinary ,Sulfaméthazin (67%), are also quantified. Cumulated concentrations ranged from 1 to 3178 ng/L on Kervidy Naizin and 8 to 2500 ng/L on Ille et Vilaine site. Maximum VP concentrations are associated to the transfer of one or two molecules,mostly antimicrobial. High temporal variability of the nature of the molecules and the concentration were observed during the seven storm events.

Discussion and conclusion

Specifically veterinary pharmaceutical were quantified in both sites and particularly in headwater catchments above the plant. The impact of agricultural practices on water contamination by VP was confirmed. Presence of human/veterinary molecules suggested that diffuse agricultural pollution may really contribute to the presence of pharmaceutical in water but it is difficult to identify the main sources. Keywords : veterinary pharmaceutical, catchment, agriculture, farming waste, storm event.


C Pean

B Rotenberg

P Simon

M Salanne

Understanding the different (dis)charging steps of supercapacitors: influence of potential and solvation

Keywords:

Supercapacitors are an innovative and promising technology in the field of energy storage. In the present study, we use modeling via computer simulation as a technique to study the dynamic processes of charging and discharging. This methodology is complementary to experiments. Various model systems composed of different structures of carbon electrodes, in contact with either pure or solvated ionic liquid, polarized at positive or negative potentials were investigated. From the data obtained, we identify several characteristic times for the charging mechanism of supercapacitors. Furthermore, we determine the influence of the structure of the electrode material, and the effects of potential and solvation on dynamical processes.

Introduction

Supercapacitors, also known as Electrical Double-Layer Capacitors (or EDLCs), are powerful and promising electrochemical storage systems [START_REF] Simon | Capacitive Energy Storage in Nanostructured Carbon Electrolyte Systems[END_REF]. In particular, it has been shown that nanoporous carbon-based supercapacitors exhibit very high capacitances [START_REF] Chmiola | Anomalous increase in carbon capacitance at pore sizes less than 1 nanometer[END_REF]. Various experimental techniques allow to investigate supercapacitors: among them, the most powerful are, in addition to electrochemical methods, EQCM [START_REF] Levi | Electrochemical Quartz Crystal Microbalance (EQCM) Studies of Ions and Solvents Insertion into Highly Porous Activated Carbons[END_REF], classical NMR [START_REF] Forse | Nuclear magnetic resonance study of ion adsorption on microporous carbide-derived carbon[END_REF][START_REF] Merlet | Lattice simulation method to model diffusion and NMR spectra in porous materials[END_REF], solid state NMR [START_REF] Deschamps | Exploring electrolyte organization in supercapacitor electrodes with solid-state NMR[END_REF], infrared spectroscopy [START_REF] Richey | Ion dynamics in porous carbon electrodes in supercapacitors using in situ infrared spectroelectrochemistry[END_REF], SAXS and SANS [START_REF] Boukhalfa | Small-angle neutron scattering for in situ probing of ion adsorption inside micropores[END_REF][START_REF] Humer | Scanning microwave microscopy and scanning capacitance microscopy on colloidal nanocrystals[END_REF]. For all experimental techniques though, it is difficult to get detailed and quantitative information at the molecular level. In contrast, molecular simulations allow to analyse finely the processes at the molecular scale, since it provides the position of all the atoms over time. The two methodologies are then complementary. Among the numerous computational studies that have already been published on supercapacitors, electrodes can be simulated as planar structures [START_REF] Vatamanu | Molecular Simulations of the Electric Double Layer Structure, Differential Capacitance, and Charging Kinetics for N-Methyl-N-propylpyrrolidinium Bis(fluorosulfonyl)imide at Graphite Electrodes[END_REF][START_REF] Vatamanu | Molecular Insights into the Potential and Temperature Dependences of the Differential Capacitance of a Room-Temperature Ionic Liquid at Graphite Electrodes[END_REF], slit pores [START_REF] Xing | On the Atomistic Nature of Capacitance Enhancement Generated by Ionic Liquid Electrolyte Confined in Subnanometer Pores[END_REF][START_REF] Lee | Charging dynamics of supercapacitors with narrow cylindrical nanopores[END_REF][START_REF] Kondrat | Superionic state in double-layer capacitors with nanoporous electrodes[END_REF], or regular carbon chains [START_REF] Vatamanu | Increasing Energy Storage in Electrochemical Capacitors with Ionic Liquid Electrolytes and Nanostructured Carbon Electrodes[END_REF]. All these models are idealized and the models with complex disordered electrodes are closer to experimental reality. Here, we simulate nanoporous carbide derived carbon (CDC) electrodes. In particular, with the latter it has been possible to explain the high capacitance exhibited by nanoporous carbon based supercapacitors [START_REF] Merlet | On the molecular origin of supercapacitance in nanoporous carbon electrodes[END_REF]. Complex nanoporous electrodes are more difficult to simulate, on the one hand because they are not easy to generate and on the other hand because it is computationally more demanding to work with disordered systems involving a large number of atoms. Furthermore, it is not as easy as with planar structures to measure their performances. To calculate the accurate capacitance of a nanoporous disordered electrode, it is necessary to use a planar electrode as a counter electrode with which the Poisson potential can be defined [START_REF] Pean | Single electrode capacitances of porous carbons in neat ionic liquid electrolyte at 100 • C: a combined experimental and modeling approach[END_REF]. Several studies on the equilibrium arrangement of the species inside supercapacitors electrodes have already been carried out in our group [START_REF] Merlet | Imidazolium Ionic Liquid Interfaces with Vapor and Graphite: Interfacial Tension and Capacitance from Coarse-Grained Molecular Simulations[END_REF][START_REF] Merlet | New Coarse-Grained Models of Imidazolium Ionic Liquids for Bulk and Interfacial Molecular Simulations[END_REF]. Furthermore, the effect of solvation on static properties (structure of the liquid, global capacitances, density profiles) have been studied on planar electrodes [START_REF] Merlet | Influence of solvation on the structural and capacitive properties of electrical double layer capacitors[END_REF]. Despite such significant progress, a real difficulty in the electrochemical study of supercapacitors remains to apprehend dynamical processes, either experimentally or by simulation. Kornyshev et al. investigated the transport properties inside slit pores [START_REF] Kondrat | Accelerating charging dynamics in subnanometre pores[END_REF][START_REF] He | Dynamic Charge Storage in Ionic Liquids-Filled Nanopores: Insight from a Computational Cyclic Voltammetry Study[END_REF]. Currently, only few experimental techniques give access to dynamical properties, e.g. EQCM [START_REF] Tsai | Electrochemical Quartz Crystal Microbalance (EQCM) Study of Ion Dynamics in Nanoporous Carbons[END_REF], infrared spectroelectrochemistry [START_REF] Richey | Ion dynamics in porous carbon electrodes in supercapacitors using in situ infrared spectroelectrochemistry[END_REF][START_REF] Richey | Ionic liquid dynamics in nanoporous carbon nanofibers in supercapacitors measured with in operando infrared spectroelectrochemistry[END_REF] and NMR [START_REF] Griffin | In situ nmr and quartz crystal microbalance techniques reveal the electric doublelayer structure of supercapacitor electrodes[END_REF].

In this work we use a realistic model to investigate by classical molecular simulations the dynamical behaviour of EDLCs: carbon electrodes are made of complex CDCs, in contact with either pure or solvated ionic liquid [BMI][PF 6 ]. The model allows for the polarization of the electrodes, and thus to work with a constant potential [START_REF] Merlet | Simulating Supercapacitors: Can We Model Electrodes As Constant Charge Surfaces[END_REF], as it is the case experimentally. We already showed that these model systems were reliable because they allow to calculate reasonable charging times for corresponding macroscopic devices [START_REF] Pean | On the Dynamics of Charging in Nanoporous Carbon-Based Supercapacitors[END_REF]. Here, we perform out-of-equilibrium simulations that simulate the charging and discharging of supercapacitors. Charging and discharging are initiated by a perturbation, namely a sudden change in the potential difference. The path to return to equilibrium is then analyzed. The evolution of characteristic quantities, such as the number of ions inside the electrodes, the distance of the electrolyte molecules to the carbon surface or the degree of confinement of the ions [START_REF] Merlet | Highly confined ions store charge more efficiently in supercapacitors[END_REF] can be monitored as a function of time. By comparing the behaviour of pure and solvated systems charged at different voltages, we are able to conclude on the effects of solvation and of high potentials on the dynamics of charging and discharging. Furthermore, we are able to identify a stepwise mechanism for the charging process.

Computational Details

The simulation cells consists of [BMI][PF 6 ] ionic liquid (600 ion pairs) or [BMI][PF 6 ] solvated in acetonitrile (230 ion pairs among 2146 acetonitrile molecules, corresponding to a concentration of 1.5 M that is classically used in experiments) surrounded by two identical and symmetric nanoporous carbon electrodes, obtained by quenched molecular dynamics [START_REF] Palmer | Modeling the structural evolution of carbide-derived carbons using quenched molecular dynamics[END_REF]. For the pure ionic liquid, three different CDC electrode materials are considered (see Table 1). Following our previous works, we use for the ionic liquid the coarsegrained model of Roy and Maroncelli, which yields the correct structural, thermodynamic and dynamic properties of [BMI][PF 6 ] [START_REF] Roy | Dynamics in an idealized ionic liquid model[END_REF][START_REF] Roy | An Improved Four-Site Ionic Liquid Model[END_REF]. Three and one sites describe the cation and the anion, respectively [START_REF] Merlet | Imidazolium Ionic Liquid Interfaces with Vapor and Graphite: Interfacial Tension and Capacitance from Coarse-Grained Molecular Simulations[END_REF][START_REF] Merlet | New Coarse-Grained Models of Imidazolium Ionic Liquids for Bulk and Interfacial Molecular Simulations[END_REF]. For acetonitrile, also following our previous work, the coarse-grained model is the one of Edwards et al. [START_REF] Edwards | A computer simulation study of the dielectric properties of a model of methyl cyanide[END_REF]. The size of the simulation cells in the z direction are chosen in order to reproduce the experimental densities of the bulk electrolyte.

Molecular dynamics simulations are conducted on the four different systems described in table 1, with a timestep of 2 fs which is consistent with the coarse-grained models used. The forces are calculated as the sum of site-site Lennard-Jones and coulombic interactions. Parameters for the ions and carbon atoms are the same as in our previous work [START_REF] Merlet | On the molecular origin of supercapacitance in nanoporous carbon electrodes[END_REF]. In all the simulations, the Coulomb interactions are calculated through a two-dimensional Ewald summation [START_REF] Reed | Electrochemical interface between an ionic liquid and a model metallic electrode[END_REF][START_REF] Gingrich | On the Ewald summation of Gaussian charges for the simulation of metallic surfaces[END_REF]. The lengths in the x and y directions are the same (two dimensional periodic boundary conditions are used, i.e. there is no periodicity in the z direction).

The procedure to simulate the charging process consists in two steps: The system is first equilibrated for a few nanoseconds, first with a constant charge of 0 e on all carbon atoms and then with a 0 V potential difference between the two electrodes. Then at t = 0, this potential difference 0 is suddenly set to 1, 2 or 4 V ( + = 0 /2 and -= -0 /2) and maintained constant. The discharging process is simulated the opposite way: starting from the end of the charging simulations, at t = 0, the difference 0 is suddenly set to 0 and maintained constant. The electrodes are held at constant potential using a method developed by Reed et al. [START_REF] Reed | Electrochemical interface between an ionic liquid and a model metallic electrode[END_REF] from an original proposal by Siepmann [START_REF] Siepmann | Influence of surface topology and electrostatic potential on water/electrode systems[END_REF]. Simulations are continued enough time to reach equilibrium, i.e. to exhibit a long plateau for the value of the total charge on carbon atoms. The constant potential approach is computationally expensive compared to the constant charge approach but is compulsory to study realistically dynamic processes [START_REF] Merlet | Simulating Supercapacitors: Can We Model Electrodes As Constant Charge Surfaces[END_REF].

The simulations are performed in the NVE ensemble. The pure ionic liquid and the solvated one are simulated at 400 K and 298 K respectively. These temperatures are chosen because of the very high viscosity of [BMI][PF 6 ] at room temperature (261.4 mPa.s for [BMI][PF6] [START_REF] Tokuda | Physicochemical Properties and Structures of Room Temperature Ionic Liquids. 1. Variation of Anionic Species[END_REF]) and because acetonitrile boils at 355 K. Furthermore, at these temperatures we have equivalent mobilities for the ions. Even if the simulated temperatures are different, qualitative comparisons can be done between solvated and pure ionic liquids. The qualitative conclusions on the charging and discharging mechanisms are still valid because we identify mechanisms that are found for different temperatures. Indeed, high potential simulations (4 V) exhibit temperature elevation along the simulation, explained by the Joule effect [START_REF] Merlet | Simulating Supercapacitors: Can We Model Electrodes As Constant Charge Surfaces[END_REF]. The data presented in this article show that mechanisms are qualitatively the same for low or high potential differences, thus proving the weak effect of temperature here.

Results and Discussion

Charges on carbon atoms

Fig. 1 shows the evolution of the total charges on the carbon atoms of either the positive or the negative electrode (both are opposite because of the electroneutrality of the global system) during charging processes at different voltages or discharging on CDC-1200. A higher potential difference allow to reach a larger plateau value at the end of the simulation. During the charging, the total charge of the positive electrode increases, first very fast, and then more slowly, for all potential differences, according to our previous results [START_REF] Pean | On the Dynamics of Charging in Nanoporous Carbon-Based Supercapacitors[END_REF]. In the latter work, we describe the use of an equivalent circuit model to interpret the evolution of the total charges on carbon atoms. The charging plots were fitted accurately by a biexponentional function. This work allowed to demonstrate that our simulations are performed on a time scale that is consistent with the size of our model systems. Here, in order to make comparisons between characteristic times for different quantities, we choose the simplest model of a monoexponential fit with the general expression for the charge during charging:

Q(t) = Q max . (1 -exp(-t/ )
). Indeed, if the principal aim is to compare different systems and extract their qualitative behaviour, the fit does not need to be extremely accurate. From the values of the characteristic times presented on the top of Fig. 1, we can conclude that the higher the potential is, the faster the stationary state is reached. Furthermore, the integral capacitance C int = Q tot only slightly depends on the applied potential difference. We also note that when the potential is high, the monoexponential function really suits the raw data.

Then we compare the values of the for the charging of CDC-1200 in contact with either pure or solvated ionic liquid (Fig. 1, bottom). In both cases, we have ≈ 970 ps, revealing the same kinetics for the charging process at this concentration of acetonitrile of 1.5 M and at these simulated temperatures (298 K for acetonitrile

Table 1

Lengths and number of molecules for each simulation cell. The lengths in the x and y directions are the same. Two dimensional periodic boundary conditions are used, i.e. there is no periodicity in the z direction. N ions corresponds to the number of ion pairs and N ACN to the number of acetonitrile molecules. and 400 K for pure ionic liquid). Here, the solvent does not have an effect on the charging dynamics at the electrode scale. The plateau values reached at the end of the simulation are the same, indicating that the capacitance are close for the two systems, in agreement with our previous results [START_REF] Merlet | Influence of solvation on the structural and capacitive properties of electrical double layer capacitors[END_REF] where the gravimetric capacitance was of 87 F/g for the system with pure ionic liquid and 80 F/g for the system with solvent. We can then compare the charging and the discharging (Fig. 2). At first glance, the two processes seem really symmetric, the charging ending at around 25 e and the crossing between charging and discharging occuring at 12.5 e. To analyze in further detail the two mechanisms, we separate each electrode into four slices of same volume, and sum slice per slice the values of the charges for both the charging and the discharging. The results of the sum consists in four globally constant plots. This methodology allows us to compare the kinetics of the charging and discharging phenomena in a trustful way because we do not assume any mathematical function to describe the evolution of the charges on atoms. The four plots obtained are constant, which reveals that the charging and the discharging happen with the same kinetics.

Systems

Number of atoms inside the electrodes

In the tables of Fig. 3, we examine the total number of ions inside the electrodes at the initial and final states of the simulation, in cases of pure ionic liquid. A non charged state is characterized by a small difference between cations and anions (|C -A| ≤ 5), whereas charged states are characterized by a larger difference which depends on the applied potential difference. The sum C + A is not particularly important with high potential differences. Furthermore, it fluctuates with no clear tendency (the standard deviation is around 5) between initial and final states, confirming the exchange of ions between the electrode and the bulk during the charging process, already highlighted in our previous work [START_REF] Merlet | On the molecular origin of supercapacitance in nanoporous carbon electrodes[END_REF].

For CDC-1200, during the charging process at 1 V, the number of counter-ions increases by 20 -35% for the anions, by 10 -20% for the cations, and the number of co-ions decreases by 20% for the anions, and slightly less than 10% for the cations. Cations are therefore less mobile than anions, according to their more complex geometry and their larger size. Symmetrically, for the discharging, we observe larger variations for anions (by 20 -25%) than for cations (a little more than 10%). The final state of the discharging is similar to the initial state of the charging, even though not identical, consistently with the fact that the latter has been obtained by the constant charges methodology [START_REF] Merlet | Simulating Supercapacitors: Can We Model Electrodes As Constant Charge Surfaces[END_REF].

We now consider the simulations with solvent (tables of Fig. 4). The number C -A is very similar in the cases with or without solvent, according to the values of the charges on electrodes (Fig. 1). The number C + A is smaller in the presence of solvent, because there are less ions in the simulation cells with acetonitrile. It is therefore possible to have the same value of the charge on the electrode with fewer ions, if we add solvent to the system. As the charging and discharging plots of Fig. 1 are very similar with or without acetonitrile, one hypothesis would have been to think that acetonitrile is not present in the electrodes. From the number of molecules in the tables of Fig. 4, we see that it is not the case at all. Inside the electrodes, 80% of the molecules are acetonitrile molecules, i.e. only marginally less than in the bulk (90%). Furthermore, we notice that, along the charging and discharging processes, the number of acetonitrile molecules stays almost constant, which was also recently evidenced from combined EQCM/NMR experiments [START_REF] Tsai | Electrochemical Quartz Crystal Microbalance (EQCM) Study of Ion Dynamics in Nanoporous Carbons[END_REF][START_REF] Griffin | In situ nmr and quartz crystal microbalance techniques reveal the electric doublelayer structure of supercapacitor electrodes[END_REF].

The evolution as a function of time of the number of anions and cations, for a charging process at a 4 V potential difference, is presented in Fig. 5. As for smaller potential differences, two regimes are observed: a first one where the number of ions changes quickly, and a second regime where it evolves slowly. As done above for the charges on the electrodes, it is possible to fit the evolution of the number of ions by a monoexponential function. The characteristic times are given on Fig. 5. They are all of the same order of magnitude (around 300 ps) and very close to that obtained for the charging at 4 V. The two phenomena (total charge on electrode and number of ions evolution) therefore happen with the same kinetics. This is not surprising since the total charge on one electrode directly compensates the charge of the difference C -A. If we now look accurately at the evolution of different characteristic times for the number of ions, we observe that the entrance of counter-ions occurs faster than the departure of co-ions.

The evolution of the number of ions as a function of time for charging simulations at a low potential difference (1 V) with and without solvent is presented in Fig. 6. When possible, the plots are fitted by a monoexponential exponential function and the associated characteristic times are given. In the case of pure ionic liquid, we find the same result as the one obtained at higher potential difference: the entrance of counter-ions occurs faster than the departure of co-ions. The exponential fit is of poorer quality with acetonitrile. It seems that the expulsion of co-ions is faster than counter-ions entrance in the presence of solvent ( ≈ 1 ns vs ≈ 100 ps), whereas it is the contrary without solvent, whatever the Fig. 6. Evolution of the number of ions inside the positive and negative electrodes during a charging process at 1 V in the cases with or without solvent. When possible, the plots are fitted by a monoexponential function (orange) and the characteristic times are given. The quality of the fit is poorer in the case of acetonitrile simulations, and also poorer than for simulations at high potential differences, as it was the case for the evolution of the total charges. voltage. Other simulations at various potential differences on solvated systems should be carried out to confirm this last result.

Density profiles

After analyzing the global quantities at the electrode scale, it is possible to perform more detailed analysis at the pore scale. It is easy to calculate ionic densities close from a planar surface [START_REF] Kirchner | Electrical double layer in ionic liquids: Structural transitions from multilayer to monolayer structure at the interface[END_REF]. In the case of the irregular surfaces of CDCs, this calculation becomes more delicate. A local accessible surface is defined by the center of a probe sphere of argon traveling on a grid along the electrode (experimentally, it is also argon that is commonly used to measure specific surface areas in the case of porous materials [START_REF] Brunauer | Adsorption of gases in multimolecular layers[END_REF]). Locally, another planar surface, parallel to the previous one, can be defined by the contact between the radii of the carbon spheres and the probe sphere: it is the so called Connolly surface [START_REF] Connoly | Analaytical molecular-surface calculation[END_REF][START_REF] Connolly | Computation of molecular volume[END_REF][START_REF] Gelb | Pore size distributions in porous glasses: a computer simulation study[END_REF]. This methodology suffers from not distinguishing the surfaces that are actually accessible to the electrolyte (open pores) from the ones that are not (closed pores). Nevertheless, there are few closed pores in the investigated electrodes here, and not taking them into account just adds a small fixed shift on the distance of the ions to the surface. As we use data to make comparisons and not to conclude on absolute values, this does not influences our results.

In Fig. 7 we show the density profiles of anions and cations as a function of the distance to the carbon atoms surface, for different chosen times. Looking at the y-axis, we note that counter-ions are more numerous than co-ions, according to the polarization of the electrode. During the charging process, the counter-ions density increases keeping the same distance to the carbon. At 1 V potential differences (data not shown here), the co-ions density also decreases keeping the same distance to the carbon. The density maxima of counter-ions and co-ions are located at almost the same distance to the carbon (the cation seems in some cases to get closer from the surface than the anion, according to our previous results [START_REF] Merlet | On the molecular origin of supercapacitance in nanoporous carbon electrodes[END_REF], but this could be an artefact caused by the coarse-grained model). Looking at the kinetics of evolution of the density peaks, we notice that most of the changes occurs during the first nanosecond, as it is the case with the number of ions evolution. After ≈ 6 ns, the peaks do not evolve anymore. Again, we find the two regimes during the dynamical processes, one very fast, and the other slower. The same behaviour is obtained for the negative electrode. The analysis of the density profiles during the discharging processes gives the same results in line with the symmetric evolution of the total charges on electrodes.

For the charging at 1 V, the density of ions evolves with the peak staying at the same distance to the carbon. We also observe that behaviour for the counter-ions density at 4 V. However, at 4 V, the co-ions density not only decreases (the integral of the plot at 2900 ps is inferior to the one at 100 ps) but also moves away from the carbon surface. A similar smaller shift is observed in simulations at 2 V. At high potential differences, co-ions have the possibility to move away from the carbon surface, which facilitates the access of the counter-ions and thus their adsorption on the polarized electrode [START_REF] Pean | Confinement, desolvation, and electrosorption effects on the diffusion of ions in nanoporous carbon electrodes[END_REF]. There is then a different organization of the ionic liquid at high potential differences. At 0 V, the charge on the electrodes vanishes on average, counter-ions and co-ions can coexist close to the carbon without any preferential organization besides the one caused by steric effects between charged molecules. With a relatively high potential, co-ions move away from carbon to facilitate the approach of counter-ions and because the polarization of the electrode is unfavourable to them. Nevertheless, we do not observe two distinct layers of counter-ions and co-ions because their density peaks at 2900 ps are located at 0.62 nm and 0.76 nm respectively, leading to a 0.14 nm difference, which is smaller than the sum of the radii of anions and cations (equal to 0.56 nm). The same behaviour is obtained for the co-ions in the negative electrode. These conclusions on density profiles can also be found for CDC-950, which has the same average pore size of CDC-1200 (0.9 nm [START_REF] Palmer | Modeling the structural evolution of carbide-derived carbons using quenched molecular dynamics[END_REF]), but a more irregular distribution of pore sizes (it presents a large pore in the middle of the electrode). This suggests that the average pore size is the structural parameter of main importance. In CDC-800, at 4 V, we observe the removal of co-ions, but to a lesser extent than the other two electrode structures. Indeed, CDC-800 has a smaller average pore size (0.75 nm [START_REF] Palmer | Modeling the structural evolution of carbide-derived carbons using quenched molecular dynamics[END_REF]), thus preventing species to really move away or get closer to the carbon surfaces. Again, the analyses of the discharges give the same results, confirming the symmetric mechanism for the discharging process.

For the charging with solvent at 1 V, we observe a slight shift of the density peak of co-ions (Fig. 8). This effect is more visible on the positive electrode, which mainly attracts anions, less voluminous and more mobile than cations. The density profiles for acetonitrile do not vary as a function of time, proving again that acetonitrile as a neutral molecule is not really affected by the polarization of the systems. Overall, adding solvent to pure electrolytes facilitates the charging, avoiding to apply a high voltage, because it allows the co-ions to move away from the surface.

Looking at the pore scale

In order to learn more about this removal of species from the carbon surfaces, we use our previous definition of the degree of confinement [START_REF] Merlet | Highly confined ions store charge more efficiently in supercapacitors[END_REF] distinguishing four different types of pores, from less confined to more confined: edges, planes, hollows, pockets. Fig. 9 presents the evolution of anions or cations in each type of pore during the charging of CDC-1200 at 4 V with pure electrolyte. Counter-ions move towards high confinement sites, whilst co-ions move towards low confinement sites. The latter can never access pockets, for example. Furthermore, we can note the following result: co-ions go mostly from planes to edges, and counter-ions go mostly from edges to hollows and pockets. The fact that it is possible to identify distinct co-ion and a counter-ion behaviours suggests that these conclusions could be generalizable to other Fig. 9. Percentage of anions and cations in the four different confinement sites [START_REF] Merlet | Highly confined ions store charge more efficiently in supercapacitors[END_REF] in the case of charging at 4 V for pure ionic liquid, and percentage of acetonitrile molecules in the case of charging at 1 V for solvated electrolyte systems. When possible, the plots are fitted by a monoexponential function (dashed lines).

ionic liquids, at least in the cases where there is a relative morphology difference between anions and cations.

When it is possible, the plots of Fig. 9 are fitted with a monoexponential function (indeed, not all observed quantities behave in a monoexponential way). From the corresponding characteristic times, we conclude that the reorganization of counter-ions is faster than that of the co-ions. Comparing these characteristic times to the ones already obtained in all this study, we can also conclude that the reorganization of counter-ions process is the fastest event. The analyses of the reorganization for the discharging processes leads to the same conclusions. For low potential differences, the same behaviours are observed, but the processes are slower, as for the charges or the number of ions inside the electrodes. Furthermore, the same behaviours are obtained for the three studied CDCs. For solvated systems, the plots are noisier, but we can also identify that the reorganization of co-ions occurs between edges and planes and that the one of counter-ions occurs mainly between planes and hollows. The pockets are almost never occupied by the ionic liquid because they contain the solvent, which is smaller (r ≈ 0.242 nm whereas for the cation, r ≈ 0.31 nm, and for the anion, r ≈ 0.253 nm). Finally, it is also interesting to note that for acetonitrile, there is no preference for a certain type of site. The proportion of acetonitrile molecules in each site is constant in time for either charging or discharging processes, according to what we already concluded for neutral molecules.

Conclusions

The results presented in this article allow us to conclude that charging and discharging are symmetrical processes that happen with the same kinetics. There are two regimes of response for the different quantities observed: first, a very fast regime and then, a slower second one. Furthermore, a high potential difference between electrodes accelerates the dynamical processes. The addition of solvent changes slightly the charging mechanism at the pore scale even if no major difference is seen at the electrode scale. A parallel can be drawn between the effects of increasing potential and adding solvent: the co-ions displacements is more evident and in particular their distance to the carbon electrodes can change more easily. In addition, we characterized the kinetics of the phenomena happening either at the pore scale or the electrode scale. Characteristic times have been identified, allowing to establish a classification for the charging mechanism, confirmed by the study of three different structures of carbon electrodes at various potentials:

-the inter-site reorganization of counter-ions is faster than that of co-ions -the counter-ions entrance into the electrode is faster than the exit of co-ions, with pure ionic liquids. Other simulations with various solvent concentrations submitted to various potential differences would be necessary to confirm that the mechanism is opposite with solvent. -the total charge on the electrode and the total number of ions inside the electrode evolves with the same characteristic time because the total charge on the electrode exactly compensates that of the ions inside the electrode, and reciprocally

This work provides knowledge on the behaviour of [BMI][PF 6 ] ionic liquid based supercapacitors, and as emphasized in the article, some co-ions and counter-ions behaviours suggest that the results could be the extendable to other ionic liquids. Understanding precisely the functioning of supercapacitors is a major step for improving performances of these devices.

Fig. 1 .

 1 Fig. 1. Top: evolution of the charge on the positive electrode of CDC-1200 in contact with pure [BMI][PF6] for charging processes at three different potential differences: 1 (black), 2 (red) and 4 V (blue). Bottom: evolution of the charge on the positive electrode of CDC-1200 during charging or discharging processes in the case of pure [BMI][PF6] (black) or solvated in acetonitrile (red). The monoexponential fits are in orange.

Fig. 2 .

 2 Fig. 2. Sum slice to slice (slice 1 correspond to the deepest slice in the electrode, slice 4 is the slice at the interface) of the values of the charges on the carbon atoms on the positive electrode during charging and discharging processes. The four plots obtained are globally constant, which suggests that the two processes are symmetric.

Fig. 3 .

 3 Fig. 3. Number of ions (A: anions; C: cations), their difference and their sum, inside positive and negative electrodes, at the initial and final states of the simulations of charging at 1 V, 2 V or 4 V and the corresponding discharging, for pure ionic liquid, and for different CDCs. Results are averaged over 1 ps of simulation. The standard deviation is around 5.

Fig. 4 .

 4 Fig. 4. Number of molecules (A: anions; C: cations; ACN: acetonitrile), their difference and their sum, inside positive and negative electrodes, at the initial and final states of the simulations of charging at 1 V and the corresponding discharging, for solvated ionic liquid. Results are averaged over 1 ps of simulation. The standard deviation is around 5.

Fig. 5 .

 5 Fig. 5. Evolution of the number of ions inside the positive and negative electrodes during a charging process at 4 V for the pure ionic liquid. The plots are fitted by a monoexponential function (orange) and the characteristic times are given.

Fig. 7 .

 7 Fig. 7. Density profiles of anions (left panel) and cations (right panel) inside the positive electrode (top) and the negative electrode (bottom) as a function of the distance to the carbon surface, for charging simulations at 4 V, in the case of pure ionic liquid, at 6 chosen times. Each plot is obtained by averaging the density over 100 picoseconds. The black plot corresponds to the beginning of the simulation and that in light blue corresponds to the end of the simulation.

Fig. 8 .

 8 Fig. 8. Density profiles of anions (left panel) and cations (right panel) inside the positive electrode (top) and the negative electrode (bottom) as a function of the distance to the carbon surface, for charging simulations at 1 V, in the case of solvated ionic liquid, at 5 chosen times. Each plot is obtained by averaging the density over 100 picoseconds. The black plot corresponds to the beginning of the simulation and that in orange corresponds to the end of the simulation.

Acknowledgments

We acknowledge the support from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC grant Agreement n.102539 (Advanced grant, Ionaces project). We are grateful for the computing resources on OCCIGEN (CINES, French National HPC) and CURIE (TGCC, French National HPC) obtained through the project x2015096728.