Reuven Segevt

Amit Ailon

A Geometrical Setting for the Newtonian Mechanics of Robots

A geometrical setting for the Newtonian mechanics of mechanical manipulators is presented. The configuration space of the mechanical system is modelled by a differentiable manifold. The kinematics of the system is formulated on the tangent and double tangent bundles of the corifiguration space, and forces are defined as elements of the cotangent bundle. The dynamical properties of the system are introduced by specifying a Riemannian metric on the configuration space. The metric is used in order to generate the generalized momenta and the kinetic energy from the generalized velocities, and the connection it induces makes it possible to formulate a generalization of Newton's second law relating generalized forces and generalized accelerations.

I. Introduction

,

Local properties of mechanical systems such as the equations of motion can be described by the classical treatment of mechanics. However, these formulations cannot provide a framework for the discussion of global questions such as the existence and uniqueness of the solutions of the equations of motion. While the classical control theory deals with differential equations in Rn, the motion of a mechanical system takes place (in the holonomic case) on differentiable manifolds.

Comprehensive formulations of Lagrangian and Hamiltonian mechanics on differentiable manifolds have been available for several years, see, for example [START_REF] Abraham | Foundation of Mechanics[END_REF][START_REF] Arnold | Mathematic al Methods of Classical Mechanics[END_REF]. These general formulations, which are based on the notion of a symplectic manifold, cannot be utilized directly in engineering applications because, as their names indicate, they assume that the motion is given in terms of a Lagrangian function or a Hamiltonian function, and the Newtonian notion of a force is missing.

In this paper, we present a global geometrical formulation of the dynamics of multi-degrees of freedom mechanical systems which incorporates the notion of a force. Such a formulation is suitable for the analysis of the dynamics of robots where (generalized) forces serve as control parameters. In the suggested formulation the equations of motion are written directly on the basis of a generalized invariant version of Newton's second law. For the sake of completeness and easy reference, t Fellow of the B. De Rothchild Foundation for the Advancement of Science in Israel.

~ we present in an informal manner all the necessary mathematical terms and we define all the mathematical notions in the context of mechanics.

Section II introduces the notion of a differentiable manifold as the configuration space of a mechanical system and Section III introduces the tangent bundle and the double tangent bundle through the notions of motion and generalized velocity. Section IV introduces the cotangent bundle of the configuration space and presents the statics of a mechanical system as an example. In Section V we introduce the Riemannian metric (kinetic energy), and the structure it induces on the configuration space. In Section VI we use the structure provided by the kinetic energy in order to relate the forces and accelerations by the equation of motion. Geometrical results concerning the case of a free motion of the system are stated in Section VII. It turns out that these results allow us to write the equations of motion in an uncoupled linearized form if the gravity forces are neglected.

II. The Configuration Space

The basic notion of any analytical study of a mechanical system is the configuration space. The configuration space of a mechanical system is the collection of all possible configurations or states of the system. Clearly, the configuration space describes the kinematical properties of the mechanical system completely. The configuration spaces of mechanical systems which are composed of material particles and rigid bodies can be associated with geometric objects. For example, the configuration space of the double planar pendulum can be identified with the surface of a torus. Given a mechanical system, we denote its configuration space byQ.

In order to a~ply mathematical analysis to the configuration space, generalized coordinates are introduced. The number of generalized coordinates needed to specify any configuration is the number of degrees of freedom of the system. The following assumptions allow us to assign coordinates to the various configurations in a meaningful way (i) We assume that for every configuration q 0 in Q there exists a collection of neighboring configurations U and a function cp that assigns to any configuration q in U an ntuple of real numbers cp(q) = (ql(q), q2(q), ... , qn(q)), such that the image of cp is an open subset of Rn, and if q and q' are two distinct configurations in U, then cp(q) =1= cp(q').

The subset U is called a coordinate neighborhood and we say that cp is a chart on U and that qi(q), i = I, ... , n, are the coordinates of q under cp.

(ii) It might happen that two coordinate neighborhoods U and U' with charts cp and cp' intersect and we can have the transformation of coordinates cp' ocp-I on cp(U n U'), which gives the coordinates qi'in the chart cp' in terms of the coordinates qi in the chart cp in the form of n functions qi' = {(qi). We assume that both cp' 0 cp -I and its inverse transformation can be continuously differentiated as many times as we wish and it follows that the matrix oqi'joqi is nowhere singular.

The collection of coordinate neighborhoods and the charts defined on them are called an atlas on Q. Together with such an atlas Q has the geomctric sti-ucture of an n-dimension al differentiabl e manifold. It can be shown that the coordinate neighborho ods form a basis for a topology on Q. It can also be shown that if Q is compact in that topology it cannot be covered by a single coordinate neighborho od and it follows that the torus and the sphere, for example, cannot be covered by a single coordinate neighborho od. Let h: Q -+ Q' be a continuous mapping from a manifold Q into a manifold Q'. Then, by the definition of the topology on a manifold, for every coordinate neighborho od Von Q', there exists a coordinate neighborho od U on Q such that h( U) c V. However, denoting the charts on the two coordinate neighborho ods by <p and t/J, the unique relation between points on manifolds and their coordinates implies that the mapping t/Joho<p-1 associates the coordinates si of h(q) with the coordinates qi of any point q in U, in the form sf = sf(qi). The mapping t/Joho<p-I is called a local representati ve of h. We say that a mapping h: Q -+ Q' is differentiab le if every local representati ve of h is differentiab le. Seemingly, this definition depends on the charts used to construct the local representati ves. However, the differentiability of the coordinate transformat ions guarantees that if a local representative of a mapping with respect to some charts on Q and Q' is differentiab le, then, any other local representati ve of h with the same domain on Q will be differentiab le so that the notion of differentiab ility is well defined.

III. The Tangent Bundle

By a motion of the mechanical system having a configuratio n space Q, we mean a differentiab le mapping m: R -+ Q, where R is the time axis. Mathematic ally., such a motion is called a curve iJil Q. Let q be a configuratio n of the mechanical system and let m and m' be two curves such that m(O) = m'(O) = q. If <p is a chart in some neighborho od of q then the two motions induce the motions <pom = (ql(t), . .. , qn(t)) and <pom' = (ql'(t), . .. , if'(t)) in a neighborho od of <p(q) in Rn. We say that m and m' are tangent at q if d/dt (<pom) 11= 0 = d/dt (<pom') 11= 0, i.e. their representati ves are tangent in Rn at the coordinates of q. Again, it can be shown that if two motions are tangent when we use the chart <p they are tangent with respect to any other chart.

Let v denote a class of all motions that are tangent to a certain motion m at q. We call v a tangent vector at q or a generalized velocity at q and we say that m or any other motion in v represents v. The tangent space to Q at q is the set of all velocities at q and it is denoted by TqQ.

We show now that TqQ has the structure of a vector space. Let v and u be two velocities at Q, let <p be a chart in a neighborho od of q and let m and m' be motions representing v and u. Then, for a and b in R, we define av + bu to be the tangent vector whose representing motion is <p-I(<p(q)+ t(av+bu)), where

V= d/dt(<p0m)lt =o, u = d/dt (<pom') It=o, i.e. Vi = d/dt(qi(t))lt= o, d = d/dt ({(t)) It=o.
One can easily verify that this definition satisfies all the axioms of a vector space and that it is independen t of the chart <po

The tangent bundle to the configuratio n space is defined as TQ = UqeQTqQ. The tangent bundle projection 't"Q: TQ -+ Q is the mapping which assigns to every tangent -= vector v the configuration q where it is tangent. the tangent -bundle TQ can be given the structure of a differentiable manifold in the following way. Given any tangent vector v, let cp be a chart in a neighborhood U Qf q = LQ(V). We can define the chart <1>: LQ I (U) -+ Rn x Rn by

<I>(v) = (cp(q),d/dt(cpom)/t=o) = (ql, ..

. ,qn,d/dt(ql(t)), ... ,d/dt(qnct)))/t=o,

where m is any motion representing v and qi(t) are the components of cpom. This chart on TQ induced by the chart cp on Q is called the natural chart and it clearly satisfies the first axiom of manifolds. In addition, by the definition of the linear structure on TqQ, <l>q: TqQ -+ Rn is an isomorphism of vector spaces. Assuming that cp' is another chart in a neighborhood of q and that <1>' is the chart it induces on the tangent bundle, we relate the coordinates (l, vi) of v with respect to <I> with the coordinates (l', 1/) of v with respect to <1>' . Let qi(t) be the components of cpom, then, by definition, qi = qi(O), vi = d/dt (qi(t)) / t= o. The components of cp' am will be {(t) = qi'(qi(t)), and using the chain rule and the summation convention, the coordinates of v will be Hence, denoting the partial differentiation O{/oqi by q:;, we have vi' = q:;v i and we note that q:; is a vector space isomorphism.

Given a differentiable mapping h: Q -+ Q' there is an induced mapping Th: TQ -+ TQ' which is defined as follows. Let v be any tangent vector, let LQ(V) = q and let m be a motion representing v. Then, we define Th(v) to be the tangent vector to h(q) which is represented by the motion hom in Q'. If c((q) are the coordinates in a neighborhood of q and siCs) are the coordinates in some chart in the neighborhood of h(q), the coordinates (sr, Uk) of Th(v) in the natural chart are given by

s = s(c(), if = d/dt (s"(qi(t))) / t= 0 = s",jd,
where (q", d) are the coordinates of v in the natural chart on TQ. Thus, for a fixed q, the restriction T;: TqQ -+ Tf(q)Q' of Th to the tangent space at q is a linear mapping and given charts in neighborhoods of q and h(q), it is locally represented by the derivative of the local representative of h. The mapping Th is called the tangent mapping or the derivative of h.

The process of generating the tangent bundle can be repeated and we can construct the tangent bundle to the tangent bundle T(TQ), or T2Q. Clearly, elements of T2Q have coordinates in R 4 n, and if the local representative of a curve in TQ is of the form (qi(t), viet)), the coordinates in the natural chart of the element of T2Q that the curve represents are (qi, vi, d/dt (qk), d/dt (V S

)) / t= o. The tangent bundle projection for this second tangent bundle is LTQ: T2Q -+ TQ, and locally, it is given by (qi, vi, uk, w) -+ (qi,V i ). Given a motion m: R -+ Q there is an induced curve dm in TQ which assigns to any time to the tangent vector (velocity) The curve dm is called the lifting of m and d 2 m is called the second lifting of m. Just as dm(t) represents the velocity at time t, d 2 m(t) represents the (generalized ) acceleration of the system undergoing the motion m at the time t.

As a special case of the derivative of a mapping between two manifolds we have I the mapping TrQ: T2Q ~ TQ, which is the tangent mapping of the tangent bundle projection. Since the local representati ves of rQ are of the form (qi, vi) ~ (qi), the local representati ves of TrQ are of the form (qi, vi, uk, W S ) ~ (qi, u'). We note that if w is in the image of d 2 m for some motion m, then TrQ(w) = rTQ(w),

A vector field on Q is a mapping X: Q ~ TQ such that rQoX(q) = q. We note that the derivative of a vector field is a mapping TX: TQ ~ T2Q. If the vector field is given locally in the form (qi, Xi(qi», then the local representati ve of TXis of the form (qi, Vl) ~ (q, Xk(qi), vi, XS,,(qi)V'). A vector field on Q induces a first-order differential equation on Q in the following way. Let qo be a configuratio n. We say that the motion m on Q, with m(O) = qo, is the solution of the differential equation induced by X if its lifting dm is the restriction of X to the image of m. It follows that if (qi, Xi(qi» are the coordinates of X, the coordinates of the solution curve qi(t) satisfy the equation d/dt (ql)(t) = Xi(qi(t», which justifies the term. Similarly, a first-order equation on TQ is induced by a vector field Y: TQ ~ T2Q on TQ. A vector field Y on TQ is a second-orde r differential equation if r TQ ° Y = Tr QO Y, or equivalently , if the coordinates of Yare of the form (q, vi, vk, YS(qi». If (q(t),vi(t» are the coordinates of the solution curve on TQ to a second-orde r differential equation, they satisfy the relation d/dt (qi) = vi, d/dt (vi)(t) = yi(qi(t», so that the use of the term is justified.

IV. The Cotangent Bundle and Forces in Statics

Let qeQ and let T*Qq denote the dual vector space to TqQ. We call T*Qq the cotangent space to Q at q and we call T*Q = UqeQT*Qq the cotangent bundle of Q. An element of T*Q will be referred to as a covector, and nQ: T*Q ~ Q will denote the mapping that assigns to any covector f, the configuratio n q such that / e T*Qq. We endow T*Q with the structure of a differentiab le manifold as follows. Let q> be a chart in a neighborho od U of q, then the chart <1>*: nQ I(U) ~ Rn x Rn is defined by <1>*(/) = (qi, jj) = (qi, fee)~, where ei is the tangent vector whose natural coordinates are (0, ... , 1, ... ,0), with 1 in the jth place. In other words, if (qi, jj) are the coordinates off, thenf (v) = !jv i , where vi are the natural coordinates of v. One can easily show that if a coordinate transformat ion qi' = qi' (qi) is given in a neighborho od of q, then, j;, = j;q:j' , i.e. the coordinates of f are related by the transpose of the matrix of the partial derivatives.

Let h: Q ~ Q' be a differentiab le mapping. Note that for q e Q we have a mapping T*hh(q):T*Q~(q) ~ T*Qq which is the adjoint of Tqh and th~~efore its local representative is the transpose of the matrix of T ;.

Convectors represent mathematically the generalized/orees and the meaning of the real number obtained by the action of a force on a generalized velocity is simply the virtual work that the force performs on the virtual displacement or the power supplied by the force. Indeed, in order that the virtual work will be independent of the coordinates by which we choose to represent the velocities, the transformation of coordinates for forces should be the one we obtained for covectors. For example, let B be a rigid body and let Q be its configuration space.

Choosing a Cartesian coordinate system which is fixed to the body and a coordinate system which is fixed in space, we can find for each configuration qo a neighborhood of qo in which the rectangular components of the origin of the body frame and the rotations of the body axes about the space axes can serve as coordinates. In this coordinate system the components of the generalized velocity will be composed of the three components of the linear velocity of the origin of the body frame and the three components of the angular velocity of the body frame. The coordinates of a generalized force will be composed of the three components of the force acting on the body and the three components of the torque (by torque we mean the components of the generalized force that correspond to the rotations). Let A be a point in B. Using the space frame, the configuration space of A can be identified with the Euclidean space R 3 and the coordinates of a generalized force are simply the components of the force acting on A. We can write the mapping h: Q ~ R3 which assigns to any configuration q of the body the corresponding position of the point A in space. Then, for any configuration q, the linear mapping T*h: T*Rl(q) ~ T*Qq will assign a force acting on the body to any force acting on A. It turns out that the torque 6n B due to a force f acting on A is given by the vector product r x f where r is the position vector of point A in the body frame.

V. The Riemannian Metric and the Kinetic Energy

Let g: TQ ~ T*Q be a differentiable mapping such that for each qE Q the restriction of g to q is an isomorphism TqQ ~ T*Qq. The mapping g induces a nonsingular bilinear form K on TQ by K(v,u) = }g(v) (u). If the induced bilinear form is symmetric and positive definite g is called a Riemannian metric on Q. The dynamical properties of a mechanical system are represented by a Riemannian metric g in the following way. For every velocity v, the kinetic energy of the system at that velocity is K(v, v) = }g(v)(v) and the generalized momentum is the element

g(v) of T*Q.
Since the restriction of g to TqQ for any q E Q is a symmetric linear transformation, it is represented locally by a symmetric tensor gks called the metric tensor such that if / = g(v), the coordinates jj of / are given by gJsif, where if are the coordinates of v. The local expression for K is given by the same tensor in the form K(v, u) = }gijViu J • Since K is a symmetric bilinear form it can be given the interpretation of a scalar product so that 2K(u, v) is the scalar product between u and v, and (2K(u, u» 1/2 is the length of the vector u. One can easily verify that if qi' = l'(q) is a coordinate transformat ion, then the coordinatet r;~sfor7ria tion for the metric tensor is gi'j' = 9 ij if)' q:i"

As we saw earlier, the derivative of a vector field is a mapping from TQ into T2Q. Similarly, the acceleration is an element of T2Q. We will see in the following that one can use the metric in order that the acceleration , and any other member of T2Q, can be identified with an element of TQ. In such a case, the derivative of a vector field will be a linear mapping from TQ into TQ and the same will hold for a second-ord er differential equation.

We recall that although all the tangent spaces are isomorphic to Rn, these isomorphis ms are induced by the charts and therefore they are not unique. If however, we were given isomorphism s of the tangent spaces in a certain neighborhood of a configuratio n q to a given vector space, we could associate with any vector v in TqQ a unique vector v' E Tq,Q. Given such isomorphis ms this operation is called a parallel translation.

Let m be a motion on TQ and let w denote the element of T2Q that this motion represents. Using the parallel translation we can associate with w the element c(w) of Tqq, q = 'rQ(m(O)), which is given by c(w) = dldt (m'(t))lt=o, where m'(t) is the result of the parallel translation of m(t) to TqQ. (We note that the differentiati on makes sense because we calculate the derivative ofa curve in the vector space TqQ.) Such a mapping c: T2Q ~ TQ is called the connection mapping of the parallel translation. Given a motion m on TQ, we can use the parallel translation to define the curve Vm: R ~ TQ by Vm(t o ) = dldt (u')(t o ) = c(dm(to)), where u' is the parallel translation of m(t) from Tq,Q, q' = 'rQ(m(t)), to TqQ, q = 'rQ(m(to)).

Given a vector field X whose coordinates are (qi, Xi) and parallel translations in neighborho ods of all points in Q, we can define the covariant derivative V X: TQ ~ TQ as follows. The restriction of V X to q is a linear mapping TqQ ~ TqQ whose representing , matrix is Xij = oloqj(X'i), where Xli are the coordinates of the parallel translation of X(q') to TqQ. Given a vector field Y whose coordinates are (qi, yj) we denote the evaluation of the covariant derivative of X on Y by V yX. The coordinates of the vector field VyX are Xijyi.

Assuming that a Riemannian metric is given on Q it can be shown that there is a unique parallel translation on Q that satisfies the following properties:

(i) parallel translation preserves the kinetic energy, i.e. Kq,(u,v) = Kq(u',v'), where u' and v' are the parallel translations of u and v from Tq,Q to TqQ; (ii) if e i is the vector field in the neighborho od of q whose coordinates satisfy (Xl) = (0, ... ,1, ... ,0), where the 1 is in the ith place, then o~ej = Veiei.

Let (qi, vi, uk, w') be the coordinates in the natural chart of the element w of T2Q. It can be shown that with the unique parallel translation that has the two aforementio ned properties, the coordinates of c(w) are (qi, w' + rjkviuk), where rjk = r~j = ifl(gsk,j+gis ,k-gjk,s) (ffi denote the component s of the inverse matrix of gjk)' The coefficients rjk are called the Christoffel symbols.

Let (qi(t), viet)) be the local expression for the curve m in TQ, so that the coordinate representat ion of dm in the standard chart in T2Q is given by (qi(t), viet), dldt (t)(t), dldt (VS)(t)). It follows that the coordinate representat ion of Vm is (qi(t), dldt (VS)(t) + nkvk dldt (ql)(t)). In the particular case where m = dmo, where 0= mo is a curve on Q whose coordinates are (q), the local represent~tive ofVm = V dmo is (ql(t), d2/dt2(qS)(t)+n~i d/dt(q") d/dt(qi)(t)). Similarly, any second-order ,differential equation w: TQ ~ T2Q whose local representative is (qi, vi, v\ W(qi, vi)), induces a mapping cow: TQ ~ TQ, which is given locally by (qi, W+rjkvivk). Ifmo is the solution of the second-order differential equation w, it satisfies the equation V dmo = cod2mo = cowomo, so that second-order differential equations can be written on TQ. In addition, it can be shown that with this parallel translation we have X lj = X: j + rjkX". In conclusion, we note that for a given q E Q the mapping c: T 2 Qlq ~ TqQ x TqQ x TqQ, given by w ~ (LTQ(W), TLQ(W), c(w)) is one to one and onto and it can be used to identify any element of T2Q with a triplet of elements ofTQ.

VI. The Equations of Motion

Consider the tangent bundle to the cotangent bundle T(T*Q), equipped with the projections LT*Q: T(T*Q) ~ T*Q and TnQ: T(T*Q) ~ TQ which are given locally in the form (qi,jj,V\gs) ~ (qi,jj) and (qi,jj,V\gs) ~ (qi, v k ), respectively. Letp be a member of (T*Q)q for some configuration q and consider curves m in T*Q with the property met) = q for each t. Since these curves do not leave (T*Q)q they are of the formp+f(t) where f(t)E(T*Q)q for all t, and each such curve is tangent to the line p + if (0). It follows that the collection of tangents to these curves form the vector space Tp(T*Q) and that the elements of Tp(T*Q) have natural coordinates of the form (qi,pj' 0, fk) where (qi, fk) are the coordinates of f (0). Hence, given two elements p and f of (T*Q)q, an element of Tp(T*Q) is induced by the line p + if. Given pin T*Q, we denote by Ip the mapping (T*Q)q ~ Tp(T*Q), q = nQ(p), given, by f ~ d/dt (p = lif).This map is clearly an isomorphism and it follows that any force f can be identified with ~/f), where Oq is the zero element of (T*Q)q, and conversely, any element of ToiT*Q) (the coordinates of such an element are in the form (qi. 0, O,jj)),induces a unique force whose natural coordinates are (qi, jj).

In general, there is no natural mapping which acts as I;; I to associate two members of T*Q with an element of T(T*Q). However, given a Riemannian metric, the connection mapping c induces a connection mapping c*: T(T*Q) ~ T*Q given by c* = gocoTg-1 so that the following diagram is commutative.

2 Tg T Q ------+) T(T*Q) e1
Ie.

TQ _-.:::..g ~l T*Q Let (qi,pj' v k , Ps) be the coordinates of an element p of T(T*Q), then the coordinates of Tg-l(p) are given by (qi,gi'pj,v\g~;vlp,+gS'p,) so that the coordinates of C(Tg-l(p)) are (qi, r~kgj'pjVk+g~;vlp,+gS'p,). Thus, the coordinates of c*(p) are (qi, g/sr:kgi'pjVk + g/sg~;vlp, + g/sgS' p,) = (qi, g/sCkgi'PjV k -if' g/S,IV1p, + pD.

Here again c* induces an identification C*:

T(T*Q) ~ T*Q x TQ x T*Q which is given by p ~ (LT*Q(P), Tn(p), c*(p)).
Let v E TqQ and let p = g(v) be the momentum associated' ;ith th;' generalized velocity v. Then, given a force / E (T*Q)g, the effective force acting on the system due to / at the velocity v is C*l(p, v,/) E T(T*Q),q' Assuming that (qi, vj) and

(ql,i1J are the coordinates of v and f, respectively , the coordinates of the effective force due to / at the velocity v are (qi, girvr, vk, i1-glsCtVrVt + gls,tVV). Using the symmetry of the product vtv s and the relation between r: t and gj" the last expression reduces to (qi, gjrVr, v\ h+ glsr:tvrvl A/orcefiel dfis a mapping TQ -+ T*Q such that 'lrQ(f(v» = LQ(V). Given a force field f and a Riemannian metric g on the configuratio n space Q, a second-orde r differential equation p: TQ -+ T2Q is induced on Q by

p(v) = Tg-1oC*-I( g(V), v, f(v».
Assuming that the local representati ves off are given in the form (qi, jj(qk, VS», one can calculate the components of Tg-1oC*-I( g(V), v, f(v» to obtain (qi, vj, vk,gsrj,_ CtvrV t ) so that p is indeed a second-ord er differential equation. We say that p is the acceleration field induced by the force fieldf, and we note that if the motion m is a solution of this second-orde r differential equation, its local representati ves From this last expression it follows that if m is a solution of the acceleration field induced by f, then m satisfies g(Vdm(t» = f(m(t» which is the required generalizati on of Newton's second low.

VII. Free Motion and Norl{1la/ Coordinates

By a free motion or a geodesic we mean a motion of the mechanical system under the zero force field. It follows that a free notion satisfies the equation Vdm = 0, and in coordinates

d2/dt2(qi)(t )+r~s d/dt(qr)d/d t(if)(t) = O.
The following properties of free motion can be proved [see [START_REF] Arnold | Mathematic al Methods of Classical Mechanics[END_REF][START_REF] Sternberg | Lectures on Differential Geometry[END_REF]] for the case of a compact configuratio n space. a. For every configuratio n q and every velocity VE TqQ there exists a unique free motion met) such that m(O) = q and dm(O) = v.

b. Let a: R -+ R be a change of scale of the time domain in the form aCt) = at+j3.Then , if m:R -+ Q is a free motion, i.e. a solution of Vdm = 0, satisfying the initial condition m(O) = v, the motion moa is also a free motion satisfying the same initial condition. In addition, the free motion m' satisfying the initial condition yv can be obtained by a reparametri zation of m and it satisfies m'(t) = m(yt), dm'(t) = ydm(yt). c. The kinetic energy is conserved along a free motion, i.e.

d/dt (K(dm(t), dm(t») = O.
d. Consider the mapping <1>: R x TQ -+ TQ, defined by <1>(t, v) = dm(t), where m is the unique free motion whose initial condition is ~. This ni1.pping, called the geodesic flow, can be shown to be differentiable.

The differentiable mapping exp: TQ ~ Q defined by exp (v) = m(1) = LQ(<I>(1,V)) is called the exponential mapping. The restriction eXPq: TqQ ~ Q of the exponential mapping to TqQ induces a one to one and onto differentiable mapping from a neighborhood V of the zero vector in TqQ onto a neighborhood U of q in Q. Thus, choosing a basis in TqQ we can construct a chart on U in the following way. For every point q' in U let its coordinates Vi be the components of the vector v in V such that q' = exp (v) = m(l) (m is the free motion with initial condition v). This chart in the neighborhood U of q is called the normal chart. In particular, we can choose an orthonormal basis in TqQ so that gij = bij (the Kronecker b) in the normal chart. It can be shown that in the normal chart we have rjk = o. Thus, using the compactness assumption, it follows that Q can be covered by a finite number of normal charts in which the equations of motion are of the form d 2 q i dt 2 = j;. e. The configuration space Q is geodesic ally complete, i.e. for any two configurations q and q' there exists a free motion which starts at q and ends at q'. f. Given the configurations q and q', the functional L(m) = .r' (K(dm(t), dm(t)) 1/2 dt which can be interpreted as the length of the curve joining the two configurations, and the functional E(m) = r (K(dm(t), dm(t))dt are at a local minimum if m is a free motion between q and q'.

VIII. Conclusions

In the preceding sections, we presented the Newtonian mechanics of mechanical manipulators from a geometric point of view. Generalized forces that serve as the control parameters in the control of robots were introduced as elements of the cotangent bundle of the configuration space. The kinetic energy was modelled by a Riemannian metric and the isomorphism it induces between the tangent and the cotangent spaces was used in order to obtain the generalized momentum corresponding to a given generalized velocity. The generalized accelerations, which are elements of the second tangent bundle, were related to the forces by using the connection induced by the Riemannian metric (kinetic energy). To the best of our knowledge this is the first geometric invariant formulation of the equations of motion for a multi-degrees of freedom mechanical system in the form of Newton's second law. In addition, the standard form of the equation~ ;;[ motion found in the literature, i.e. M(q)q+B(q ,q)q = f which is written in a coordinate non-invaria nt form, can be interpreted now as follows: M is the matrix of the metric tensor (kinetic energy) and B is composed of the covariant Christoffel symbols r ijk = girrSk of the connection induced by the metric.

The introductio n of the normal coordinates provides us with a procedure for the linearizatio n and decoupling of the equations of motion in the case where we do not consider generalized forces due to gravity [ef [START_REF] Koditschek | Robot kinematics and coordinate transformati ons[END_REF][START_REF] Hunt | Design for multi-input nonlinear systems[END_REF]]. We note that the configuratio n spaces of most mechanical manipulato rs are compact differentiab le manifolds. For example, the configuratio n spaces of mechanical systems composed of rigid bodies that are joined together so that one of the bodies rotates about a fixed point and the rest of the degrees of freedom of the system are rotations, are compact manifolds. It follows that the results of the previous section regarding compact manifolds are all applicable.

  to m at time to. If (qi(t» are local representati ves of this motion, the local representati ves of dm are (qi(t), d/dt (qi(t»). Similarly, m induces a curve d 2 m on T2Q which assigns to any time to the acceleration at that time, i.e. if (qi(t» are the local representati ves of m, the local representati ves of d 2 m are (qi(t), d/dt (qi(t», d/dt (qi(t», d 2/dt(qi(t»).


Flavien Vucko

Cédric Bosch

Asdin Aoufi

David Delafosse

Research report # ENSMSE-SMS-2014-01 Palladium coating on quenched-tempered martensitic steel for hydrogen electrochemical permeation tests

Keywords: Palladium coating, High strength steel, Hydrogen permeation

A coating process is proposed to obtain a thin palladium layer on high-strength steel. The quality of the layer is investigated by SEM observation.Potentiokinetic curves are performed to ensure no hydrogen formation during the coating process. Coated samples are tested by hydrogen permeation test and compared to non-coated samples. Permeation data analysis and boundary conditions are discussed based on permeation results and 1D simulations.

I. Introduction

Hydrogen mobility and trapping in metals can be study by electrochemical permeation (EP) of hydrogen. This method was firstly developed by Devanathan and Stachurski in 1962 [1] and it has been widely used on several metals [2][START_REF] Manolatos | The electrochemical permeation of hydrogen in steels without palladium coating. Part II: Study of the influence of microstructure on hydrogen diffusion[END_REF][START_REF] Manolatos | Necessity of a palladium coating to ensure hydrogen oxidation during electrochemical permeation measurements on iron[END_REF][START_REF] Manolatos | The electrochemical permeation of hydrogen in palladium: Boundary conditions during a galvanostatic charging under low charging current densities[END_REF][START_REF] Zakroczymski | Electrochemical determination of hydrogen in metals[END_REF][START_REF] Zakroczymski | Adaptation of the electrochemical permeation technique for studying entry, transport and trapping of hydrogen in metals[END_REF][START_REF] Hadam | Absorption of hydrogen in tensile strained iron and high-carbon steel studied by electrochemical permeation and desorption techniques[END_REF][START_REF] Frappart | Study of the hydrogen diffusion and segregation into Fe-C-Mo martensitic HSLA steel using electrochemical permeation test[END_REF][START_REF] Frappart | Hydrogen solubility, diffusivity and trapping in a tempered Fe-C-Cr martensitic steel under various mechanical stress states[END_REF][START_REF] Parvathavarthini | Influence of microstructure on the hydrogen permeability of 9% Cr-1% Mo ferritic steel[END_REF][START_REF] Ramunni | Interaction of hydrogen with the microstructure of low-carbon steel[END_REF][START_REF] Serna | Electrochemical technique applied to evaluate the hydrogen permeability in microalloyed steels[END_REF][START_REF] Hurtado Noreña | Effect of microstructure on hydrogen diffusion and trapping in a modified 9%Cr-1%Mo steel[END_REF][START_REF] Rivera | Hydrogen trapping in an API 5L X60 steel[END_REF].It is very accurate [START_REF] Coudreuse | Fragilisation des aciers par l'hydrogène: étude et prevention[END_REF] as long as boundary conditions are perfectly known. Indeed, hydrogen entry and detection are modified by surface processes. For steels, the oxide layer on the detection side acts as a barrier and limits the stabilised hydrogen flux[2, [START_REF] Manolatos | The electrochemical permeation of hydrogen in steels without palladium coating. Part II: Study of the influence of microstructure on hydrogen diffusion[END_REF][START_REF] Casanova | The influence of an oxide layer on hydrogenpermeationthroughsteel[END_REF]. But the stability of this oxide layer is a proof that hydrogen oxidation takes place at the oxide/metal interface [START_REF] Casanova | The influence of an oxide layer on hydrogenpermeationthroughsteel[END_REF][START_REF] Ruetschi | Cation Vacancies in Mn02 and their influence on electrochemical reactivity[END_REF]. So there is a modification of the diffusion of hydrogen in the metal by the hydrogen ion in the oxide layer [START_REF] Song | Hydrogen Permeation Through the Passivation Film on Iron by Time-Lag Method[END_REF]. In addition to this, complete oxidation of the hydrogen flux at the detection side is not fulfilled on iron or steel surfaces [START_REF] Manolatos | Necessity of a palladium coating to ensure hydrogen oxidation during electrochemical permeation measurements on iron[END_REF].

Moreover, during the desorption transient, due to the oxide barrier at the detection side, hydrogen desorption by the charging side increases, and detected hydrogen decreases. On the charging side, the oxide layer spontaneously formed before EP test is reduced during hydrogen charging [START_REF] Casanova | The influence of an oxide layer on hydrogenpermeationthroughsteel[END_REF]. As hydrogen absorption rate is strongly related to the coverage of hydrogen adsorbed on the surface, any modifications of this surface change hydrogen flux.

A solution to ensure boundary conditions during EP test is to use a palladium coating [START_REF] Manolatos | A thin palladium coating on iron for hydrogen permeation studies[END_REF]. This metal is totally stable in solutions of permeation cells (NaOH 0.1M) and the diffusion coefficient of hydrogen is close the one of steels [START_REF] Coudreuse | Fragilisation des aciers par l'hydrogène: étude et prevention[END_REF].Many coating methods exist, but sometimes require a heating step during the process. As the microstructure of the steel is metastable, this kind of heat treatment is prohibited, to keep the original microstructure, i.e. the same hydrogen traps sites distribution.

An electrochemical process to obtain a thin palladium coating on high-strength steel is presented. The best coating conditions are determined thanks to Hull cell tests and potentiokinetic curves. The quality of the coating is investigated by scanning electron microscopy (SEM) observations. Finally, EP tests are performed to compare the results with or without the coating on the detection side. A diffusion model is used to determine the hydrogen concentration profile in both cases.

II. Material and experimental procedure II.1 Material

The material of the study is a high strength steel S690QL (EN 10137-2). It exhibits a tempered martensitic microstructure consisting of laths of about 200 nm width. Laths are gathered in blocks and packets which are contained in prior austenite grain with a cluster size ranging from 10 to 20 µm. The steel chemical composition is given in table 1. The tempering consisted in two hours at 550°C after austenitization at 920°C and water quenching. No retained austenite was revealed by X-ray analysis.

II.2 Coating process

An electroplating bath for plating a coating of palladium is prepared with PdCl 2 salt (5g.L -1 ), ammonia (4.2M) and water. Some millilitres of commercial surface active agent are added to improve the wettability of the solution on the steel surface and quality of the coating. Final solution pH is 12.1. This high pH is needed to displace the water reduction reaction (hydrogen formation) to more negative potentials, necessary condition to ensure no hydrogen formation during the process.

For the coating process, a Platinumelectrode and saturated calomel electrode (SCE) are used as the counter and the reference electrodes.The working electrode is a cylindrical sample of 8mm of diameter encased in polyester. The contact surface is about 0.5 cm².

Surface preparation consists in mechanical polishing with SiC paper up to 4000 grade, degrease in ethanol and in some cases anodic pickling in H 2 SO 4 (30%) solution.

II.3 Hull cell tests

II.4 Potentiokinetic curves

In addition to tests in Hull cell, potentiokinetic curves are made to determine the potential of hydrogen formation. Indeed, absorption of hydrogen during the coating process could modify the permeation curve and lead to a misinterpretation of the results.

Potentiokinetic curves are performed in the electroplating bath with the same sample geometry and electrodes as for the coating process. A scan rate of 5mV.s -1 is chosen, starting from the free potential to more negative potentials.

II.5 Electrochemical permeation test

EP tests are performed in a two-compartmental set-up consisting of a charging cell and an oxidation (or detection) cell.Both cells are filled with 0.1M sodium chloride solution deoxygenated by nitrogen.

The steel sample is a thin cylindrical membrane with a thicknessl s of about 0.6 mm and a diameter of 40 mm. The same surface preparation as for coating process is performed (mechanical polishing up to 4000 grade SiC paper).The membrane is mounted in a holder which reduces thecontact surface to a circular area of 1cm² on both faces. This holder can be independently attached to each cell. Platinum counter electrode and saturated calomel electrode are used in both cells and the steel sample is the common working electrode. A floating ground galvanostat is plugged to the charging cell and a potentiostat to the detection side.

The detection side is continuously polarized at -400mV/SCE for non-coated samples and -100mV/SCE for palladium coated surfaces. The potential is lower without the coating to limit the oxide growth at the detection side. In the first case, stabilization of the current at about 0.1µA.cm -² is attended after almost 12 hours. When the steady-state is established, instantaneously with palladium coating, the charging cell is connected to the holder and the charging side of the sample is subjected to a cathodic current (-5 to -60 mA/cm²). The recorded anodic current on the detection side is the permeation rate of hydrogen through the membrane.

II.6 Modelling

The mathematical modelling of hydrogen diffusion through a layer composed of steel with thickness l s = 0.6 mm and a thin layer based on the 1 st and 2 nd Fick's laws is analysed. Two kinds of layer are investigated. The first one is palladium coating with thickness of 230 nm and a diffusion coefficient 7.10 -11 m².s -1 [START_REF] Manolatos | The electrochemical permeation of hydrogen in palladium: Boundary conditions during a galvanostatic charging under low charging current densities[END_REF], the second one is an oxide layer of 2.6 nm thick with a diffusion coefficient of 1.4.10 -22 m².s -1 [START_REF] Song | Hydrogen Permeation Through the Passivation Film on Iron by Time-Lag Method[END_REF].

Hydrogen diffusion in steel as proposed by Sofronis [START_REF] Sofronis | Numerical analysis of hydrogen transport near a blunting crack tip[END_REF] is expressed by a nonlinear diffusion equation (Eq.1)presented below, where Ω is the lattice (steel) domain.

∀ Ω : 1 + , , 1 -, , = . , Eq.1 C T and C L are respectively the trapped and lattice hydrogen concentration. N T is the number of trapping sites. The apparent lattice diffusion coefficient D eff is calculated from Eq.2:

, , , = For a current density of 2mA.cm -², 3 minutes are required to cover the entire surface. On the contrary, at 8mA.cm -² the surface is entirely coated after 1 minute. But cracks appear as the thickness of the coating grows and the current density increases.

A potentiokinematic curve, presented on figure 3, shows the formation of hydrogen by water reduction at a potential of about -950mV/SCE. It corresponds to a current density of 5mA.cm -².

To conclude, the best electrochemical conditions to coat a thin palladium layer are from 2 to 4mA.cm -² during 1 to 4 minutes.

III.2 Analysis of the coating

Steel samples were coated with palladium at a current density 4mA.cm -² during 2 minutes. Both crosssection and surface observationswere performed by SEM (see fig. 4).

The thickness was determined by XRF and image analysis. The average value determined from both techniques is identical and equal to 230nm in these coating conditions.Adherence is tested with an adhesive tape. It could be improve for thicker coating with a heat treatment, but it has been excluded for this steel. The coating is homogeneous, dense and adhesive, so perfectly adapted for hydrogen permeation applications.

IV. Electrochemicalpermeation tests

Permeation tests were performed on coated and non-coated specimens. The first rise transients and last decay transientsare compared on figures5 and 6 respectively.

The current density at steady-state is quite more important with the palladium coating (0.53A.m -²) than for uncoated sample (0.024A.m -²). In the one hand, the oxide layer seems to act as a barrier [START_REF] Song | Hydrogen Permeation Through the Passivation Film on Iron by Time-Lag Method[END_REF]. In the other hand, the hydrogen flux oxidation seems to be more complete with palladium for two main reasons: the hydrogen oxidation reaction is easily on palladium than on iron surface [START_REF] Manolatos | Necessity of a palladium coating to ensure hydrogen oxidation during electrochemical permeation measurements on iron[END_REF], and an higher oxidation potential is applied on palladium coating due to its higher stability than the raw steel.

The rise transient (fig. 5) differs from the Fick model in both cases. Indeed, stable charging conditions are not fulfilled at the beginning of the test. Native oxide on charging side is progressively reduced by cathodic current and there is a direct relation between the permeation current and this oxidation state [START_REF] Casanova | The influence of an oxide layer on hydrogenpermeationthroughsteel[END_REF]. Palladium coating on the charging side was tested with the same deposition process to remove this phenomenon. But the coating is deteriorated quite rapidly by hydrogen bubbling. A thicker coating is necessary by probably another coating process. The concentration of reversibly trapped hydrogenis determined as proposed by Zakroczymski [START_REF] Zakroczymski | Adaptation of the electrochemical permeation technique for studying entry, transport and trapping of hydrogen in metals[END_REF](see fig. 6). The apparent lattice diffusion coefficient, used for the Fick model, is determined by partial decays and build-ups (see fig. 7).The applied current density varies from 2 to 60 mA.cm -². Results are similar for all investigated transients.

The apparent lattice diffusion coefficient is quite similar with or without the palladium coating ≈ 5 × 10 -10 2 . -1 and is consistent with literature for various ferriticsteels [START_REF] Zakroczymski | Adaptation of the electrochemical permeation technique for studying entry, transport and trapping of hydrogen in metals[END_REF][START_REF] Hadam | Absorption of hydrogen in tensile strained iron and high-carbon steel studied by electrochemical permeation and desorption techniques[END_REF][START_REF] Frappart | Hydrogen solubility, diffusivity and trapping in a tempered Fe-C-Cr martensitic steel under various mechanical stress states[END_REF][START_REF] Serna | Electrochemical technique applied to evaluate the hydrogen permeability in microalloyed steels[END_REF]. As a result, the concentration of hydrogen reversibly trapped in the uncoated and palladium coated specimens are 0.115wppm and 0.34 wppmrespectively, with the assumption that no desorption occurs on the charging side. This condition is probably not fulfilled here, and desorption on both faces should be more consistent. But to evaluate the ratio of desorbed hydrogen through each side, the hydrogen concentration profile before desorption have to be known. If hydrogen is considered to desorb similarly from entry and exit side during desorption, only 1/3 of the hydrogen amount will escape from the exit side with the palladium coating, while 1/2 will desorb for the uncoated specimen. So, the concentration of hydrogen reversibly trapped in the uncoated and palladium coated specimens are respectively 0.23 wppm and 1.02 wppm. Ratios of hydrogen reversibly trapped over lattice hydrogen are 5.2 and 2.5 respectively for uncoated and palladium plated specimens.

Over estimation of reversibly trapped hydrogen is attended with an oxide layer.

V. Conclusions

Conditions for an electrochemical palladium coating process are determined to plate a low carbon steel with a thin layer of palladium. The quality of the coating is investigated by SEM observations. Optimal conditions are determined to obtainhomogeneous, dense and adhesive layer, without producing any hydrogen during the process. Electrochemical permeation tests are compared with uncoated and palladium plated specimens. The ratio of reversibly trapped hydrogen over lattice hydrogen is evaluated, using simulation. This ratio looks overestimated without palladium coating.

VI. Bibliography

  the coating process can be investigated with Hull cell tests. The cell, presented on figure 1, has a specific geometry to obtain a variation of the current density along the sample. It depends on the applied current and is given by an abacus. An applied current of 200 mA is chosen to obtain a current density from 2 to 8 mA/cm². The time of deposition varies from 1 to 3 minutes. X-Ray Fluorescence (XRF) is used to evaluate the palladium thickness along the sample considering homogeneous palladium coating on iron-based substrate. SEM observations are performed to analyze the quality of the coating.

  the true lattice diffusion coefficient without effects of trapping.The diffusion coefficient is constant in the coating and no trapping occurs.The boundary condition at entry side is constant hydrogen concentration C 0 =15 mol.m -3 (1.9 wppm) and C=0 at exit side. The trapping is simulated with a trap density N T = 50 mol.m -3 and a trapping energy E T = 20 kJ.mol -1 . Fully implicit 1D conservative finite volume numerical simulation is performed over 2000 seconds (50 increments) with our numerical software written from the ground in C++.At each time-step the algebraic nonlinear system that takes into account the coupling between Version 1.0.a 2014.02.04 6 hydrogen concentration in the lattice and coating layer is iteratively solved by a fixed point method which converges in less than 5 iterations. Influence of the mesh size and time-step over the computed solution is investigated in order to ensure that the retained numerical solution is grid independent. 1000 elements are finally used for the steel domain and 100 for the layer.III. Palladium coating III.1 Electrochemical conditions for coatingTwo examples of SEM observations of sample tested in Hull cell are presented on figures 2a and 2b.

  evaluated in both cases by 1D finite elements simulation. Results are presented in figure8. As expected, with an oxide layer, the concentration is quite homogeneous at steady-state, whereas the profile is linear with the palladium coating. The lattice concentration for the uncoated specimen is equal to the concentration at the entry side, while it is the half for the palladium coated one. The concentration atentry side subsurface can be evaluated with permeation test thanks to the equation[START_REF] Zakroczymski | Adaptation of the electrochemical permeation technique for studying entry, transport and trapping of hydrogen in metals[END_REF]: ,0 = ∞ . . with i ∞ the steady-state permeation current density, l s the thickness of the steel membrane, F the Faraday constant and D L the apparent lattice diffusion coefficient. As a result, the lattice hydrogen concentration in the uncoated and palladium plated specimens are respectively about 0.04wppm and 0.41wppm.

[ 1 ]

 1 Devanathan, MAV, Stachurski, Z, "The Adsorption and Diffusion of Electrolytic Hydrogen in Palladium". Proc. Royal Society A: Mathematical, Physical and Engineering Sciences, 270, 90-102 (1962) [2]Manolatos, P, Jerome, M, Duret-Thual, C, Le Coze, J,"The electrochemical permeation of hydrogen in steels without palladium coating. Part I: Interpretation difficulties". Cor. Sci., 37, 1773-1783 (1995)

Table 1 :

 1 Chemical composition of steel S690QL (wt. %)

	Elements	C	Si	Mn	Mo	Nb	Ni	Ti	S
	wt.%	0,16	0,33	1,22	0,25	≤0,10	0,11	≤0,05	≤0,001
	Elements	P	V	Zr	N	B	Cr	Cu	Fe
	wt.%	≤0,02	≤0,1	≤0,1	≤0,015	≤0,005	0,29	≤0,1	Bal.
	This low-alloyed steel is susceptible to general corrosion in aqueous or acid environment. It exhibits a
	passivity state in alkaline media.


H Édi Soula
email: hedi.soula@insa-lyon.fr

Bertrand Car

Guillaume Beslon

Hugues Berry
email: hugues.berry@inria.fr

Comments to the Editor Reply to the Comment by V.P. Shkilev on "Anomalous versus slowed-down Brownian diffusion in the ligand-binding equilibrium"

In his comment in this issue of the Biophysical Journal, V.P. Shkilev [START_REF] Shkilev | Comments to the Editor. Comment on "Anomalous versus slowed-down Brownian diffusion in the ligand-binding equilibrium[END_REF] disagrees with some of the findings of our article on the influence of diffusion on the ligand-binding equilibrium (2). V.P. Shkilev makes two main claims in his comment: (i) when diffusion is uniformly Brownian (with a space-independent diffusion coefficient), "the equilibrium fraction of bound receptors [...] should depend on the diffusion coefficient" and (ii) when molecules undergo transient subdiffusion due to a Continuous-Time Random Walk (CTRW), the equilibrium fraction of bound receptors can only be equal or larger (but not smaller) than the value expected for uniform Brownian motion with comparable diffusion coefficient. While the former is a frequent misconception related to the ligand-binding equilibrium, the latter relates to an important issue in the field of (anomalous) subdiffusion-reaction couplings.

Diffusion-dependence of the ligand binding equilibrium.

Even when the movement of each molecule is a simple Brownian motion with constant diffusion coefficient, the influence of diffusion on the ligand-binding equilibrium, eq. ( 1), is non-trivial:

L + R kon -⇀ ↽- k off C (1)
where k on and k off are the global forward and backward reaction rates, respectively. An overview of the impact of Brownian diffusion on this reaction can be found in the book by J. Linderman and D. Lauffenburger (3) (see section 4.2 in their book). We recall below their major results. The central point is that the classical one step process eq. ( 1) is actually not a good way to represent this reaction-diffusion process. The ligand-binding equilibrium is in essence a two-step process, that requires first the diffusive transport of individual molecules L and R (with rate constant k + ) until their distance is small enough to allow the intrinsic binding step to take place:

L + R k+ ←→ [LR] k b -⇀ ↽- ku C (2)
This two-step process thus defines an "encounter complex"

[LR] in which L and R are close enough to bind but have not yet bound.

[LR] then can undergo an intrinsic reversible binding reaction with intrinsic reaction rates k b and k u , that do not depend on the transport process. In agreement with the celebrated Smoluchovski's rate law, the transport rate constant (or diffusion-limited rate constant) reads: [START_REF] Shoup | Role of diffusion in ligand binding to macromolecules and cell-bound receptors[END_REF], where D L and D R are the diffusion coefficients of L and R, respectively, and a is the capture radius. Now, the global forward reaction rate k on depends on the diffusion coefficient through k + according to [START_REF] Shoup | Role of diffusion in ligand binding to macromolecules and cell-bound receptors[END_REF]. Therefore, the global forward constant rate is expected to increase with the diffusion coefficients D L + D R . This can for instance be seen in figure 2A of our article [START_REF] Soula | Anomalous versus slowed-down Brownian diffusion in the ligand-binding equilibrium[END_REF] where the time needed to reach equilibrium increases when the diffusion coefficient (thus k on ) decreases. Now, when applied to the global backward constat rate k off , this approach leads to the remarkable conclusion that the global backward step is not a local process but also depends on diffusivity. This can be illustrated with the following toy example. Consider the case where the diffusion coefficients vanish for all species (D L = D R = D C = 0) and the reaction is initiated with only C complexes (no free L nor R initially). In this case, dissociated species L and R cannot appear; the dynamics is restricted to transitions between "real" complexes C and "encounter" complexes [LR]. It follows that the global dissociation rate constant k off = 0 in this case, which illustrates the fact that k off too depends on the diffusion coefficient.

k + = 4π(D L + D R )a (3,
k on = k + k b / (k b + k + ) (3,
Theoretical approaches by e.g. Lauffenburger and Linderman (3), Shoup and Szabo (4) indeed lead to k off = λk u where λ = k + / (k b + k + ) is the escape probability, i.e. the probability that an encounter complex [LR] actually dissociates into free R and L molecules. Since the global associa-tion rate constant k on has the same dependence on D via k + , the dependence on the diffusion coefficient cancels out in the ratio K D = k off /k on . Therefore, when molecules move via a (uniform) Brownian motion, the equilibrium fraction of bound receptors C eq = L T R T (K D + L T ) (where L T and R T are the total concentrations of L and R, respectively) does not depend on the diffusion coefficient, in agreement with standard equilibrium thermodynamics.

Contrarily to the claim by V.P. Shkilev in his comment (1), even for D = 0.02 or 0.01, the equilibrium fraction of bound receptors (at long times) in our simulations with Brownian motion does not depend on the value of D. We agree that this was not obvious from Figure 2A where the reaction times were set so as to evidence the difference in the time needed to reach equilibrium but was not large enough to actually show the equilibrium for the smallest diffusion coefficients. However this result was directly displayed in our Figure 3B1 that shows the equilibrium constant K D (thus, indirectly, the equilibrium fraction of bound receptors) for Brownian motion with various values of the diffusion coefficient D. This figure clearly evidenced that the equilibrium fraction of bound receptors does not depend on D when molecules move via Brownian motion.

The ligand binding equilibrium with CTRW-based transient subdiffusion.

In our article, transport by CTRW-based transient subdiffusion is obtained when the distribution of waiting times (between two movements) is given by the power law φ(τ ) = ατ -(1+α) / (∆t -α -τ -α c ) where ∆t is the simulation time step. This distribution has two parameters: the anomalous exponent α and the cut-off time τ c . When τ c → ∞ or, equivalently τ c > t max (where t max is the maximal simulation time), diffusion is anomalous at all simulation times and equilibrium can never been reach in simulations.

To overcome this issue (and to improve the biological realism of our simulations), we introduced the cut-off τ c as an upper-limit of the waiting times. In this case diffusion is transiently anomalous (subdiffusive), until a crossover time (that increases with τ c ), after which diffusion gets back to a Brownian regime (see Figure 1B in our article). One major result of our article [START_REF] Soula | Anomalous versus slowed-down Brownian diffusion in the ligand-binding equilibrium[END_REF] is that with such CTRW-based transient subdiffusion, the equilibrium fraction of bound receptors appears to vary when we vary α (Figure 2C and3B3) or τ c (Figure 2D).

Yet, as rightly pointed out by V.P. Shkilev in his comment, these two parameters define the diffusion coefficient of this transport process according to D = (∆x) 2 4 τc ∆t τ φ(τ )dτ -1 where ∆x is the lattice spacing.

Therefore our simulations suggest that, contrarily to Brownian motion, the equilibrium fraction of bound receptors changes with the diffusion coefficient when molecules move via CTRW-based transient subdiffusion.

We moreover found that the equilibrium fraction of bound receptors with CTRW-based transient subdiffusion is actually lower than the Brownian case at least for small diffusion coefficients (see small values of α and large values of τ c in our Figures 2 and3). This result disagrees with the meanfield equations obtained by V.P. Shkilev in his comment (1) (eq. 14-19).

We see four major reasons that can explain this discrepancy. First, it has already been realised at several occasions that contrarily to Brownian motion, when CTRW-based subdiffusion is coupled to the reaction terms, the macroscopic behavior of the system might be very sensible to microscopic details of the system, for instance whether reactions can occur during molecule waiting times or only in association with molecule jumps [START_REF] Hornung | Morphogen gradient formation in a complex environment: an anomalous diffusion model[END_REF][START_REF] Yuste | Reaction-subdiffusion model of morphogen gradient formation[END_REF]. Such an explanation may be at work in our case.

Alternatively, the "random trap model" used by V.P. Shkilev in his comment (and his related article [START_REF] Shkilev | A model of anomalous transport[END_REF]) gives rise to a transient subdiffusive behavior with a very specific nonlinear increase of the mean-squared displacement r 2 (t) = 6D{t + (b -1)τ (1 -exp(-t/τ ))} (eq. 24 in Shkilev [START_REF] Shkilev | A model of anomalous transport[END_REF]). Like the mean-squared displacement observed in CTRWbased transient subdiffusion (our Figure 1B), this formula gives a Brownian (linear) dependence on time for t → 0 and t ≫ τ and a non-linear behavior in between. However, the general shape of the non-linear regime is quite different from what is obtained with CTRW-based transient subdiffusion, in particular the duration of the non-linear regime can be much larger with CTRW-based transient subdiffusion (several decades as in our Figure 1B) than with the random trap model. In a context where, as already note above, the macroscopic behavior of the system is known to critically depend on microscopic details, such a discrepancy might be important.

Another possible explanation can be found in the "weakergodicity breaking" property of CTRW-based subdiffusion, for which time-averages are not equivalent to ensemble averages [START_REF] Barkai | Strange kinetics of single molecules in living cells[END_REF]. Even though the subdiffusion regime is transient in our case, the effect of the weak ergodicity breaking could still be significant. In this case, the ensemble averages we use throughout our simulations might not be adequate ways to measure species concentration.

Finally, the mean-field equations developed by V.P. Shkilev in his comment are based on a one-step process for the ligand-binding equilibrium, i.e. eq. ( 1) above (note that, to the best of our understanding, his previously published work has not focused on the reversible bimolecular reaction). As such, his mean-field analysis of the Brownian case (eq.17-19 in his comment) fails to predict the fact that the equilibrium fraction of bound receptors in the Brownian case should not depend on the diffusion coefficient. We think that theoretical analysis of the ligand-binding equilibrium as a two stepprocess including transient subdiffusive transport similar to what was achieved e.g. in Lauffenburger and Linderman [START_REF] Lauffenburger | Receptors: Models for Binding, Trafficking, and Signaling[END_REF] for Brownian diffusion (i.e. eq. ( 2) above), could help resolve the contradiction between our simulation results and the mean-field analysis.

The theoretical and simulation analysis of the coupling between complex reaction schemes and subdiffusion is still in its infancy in spite of the effort of several groups. The results obtained hitherto have already produced several results that are far from trivial compared to their Brownian counterpart (see e.g. [START_REF] Henry | Anomalous diffusion with linear reaction dynamics: from continuous time random walks to fractional reaction-diffusion equations[END_REF][START_REF] Fedotov | Subdiffusive master equation with space-dependent anomalous exponent and structural instability[END_REF]). The dependence on diffusion coefficients of fundamental biological processes like the ligand-binding equilibrium might provide another example.