Guillaume Jeanneau
email: guillaume.jeanneau@ls2n.fr

Vincent B Égoc
email: vincent.begoc@icam.fr

S Ébastien Briot
email: sebastien.briot@ls2n.fr

Experimental Safety Analysis of R-Min, an underactuated parallel Robot

The R-Min robot is an intrinsically safe parallel manipulator dedicated to pick-and-place operations. The proposed architecture is based on a five-bar mechanism, with additional passive joints in order to obtain a planar sevenbar mechanism with two degrees of underactuation, allowing the robot to reconfigure in case of a collision. A preload bar is added between the base and the endeffector to constrain the additional degrees of freedom. This article presents an analysis of the workspace and of the safety performances of the R-Min robot, and it compares them with those of the five-bar mechanism, in order to evaluate the benefits of introducing underactuation in a parallel architecture to obtain intrinsically safer robots. The geometrico-static model of the R-Min robot is formulated as an optimisation problem. The direct and inverse kinemato-static models are derived from the geometricostatic model and they allow to express the singularity conditions of the R-Min robot. An analysis of the singularity loci is carried out among the robot's workspace. A controller based on the dynamic model is proposed and experimentally validated on a prototype of the R-Min robot. Finally, the safety performances of the R-Min robot are evaluated experimentally and they are compared with that of an equivalent five-bar mechanism, using the maximum impact force as a safety criteria in accordance with recent international standards.

INTRODUCTION 1.Context and state of the art

The introduction, in the industry, of collaborative robotic cells able to safely share their workspace with humans requires finding solutions to reduce the risk to an acceptable level. Improving the safety of robots can be achieved by sensing or control strategies (reviews can be found in [START_REF] Robla-Gómez | Working together: A review on safe human-robot collaboration in industrial environments[END_REF][START_REF] Vicentini | Collaborative Robotics: A Survey[END_REF]). However, those strategies are limited by the bandwidth and the response time of the controller, making it hardly suitable to face with the highly dynamics and unpredictable nature of human displacements. Therefore, strategies based on the mechanical design to obtain intrinsically safe robots are worth being deeply investigated.

A first strategy consists in reducing the mass of the robot links in order to reduce the reflected mass at the impact [START_REF] Bischoff | The kuka-dlr lightweight robot arm -a new reference platform for robotics research and manufacturing[END_REF][START_REF] Kim | Anthropomorphic low-inertia high-stiffness manipulator for high-speed safe interaction[END_REF][START_REF] Tsumaki | A 7-dof wire-driven lightweight arm with wide wrist motion range[END_REF] and thus to mitigate the impact severity. However, when applied to a serial architecture, this strategy is limited by the rigidity of the resulting architecture and the accuracy of the end-effector's position. Another strategy consists in using intrinsically compliant actuators such as SEA [START_REF] Pratt | Series elastic actuators[END_REF] (Series Elastic Actuators), VSA [START_REF] Bicchi | Fast and" softarm" tactics [robot arm design[END_REF] (Variable Stiffness Actuators) or the DM 2 (Distributed Macro-Mini) actuation approach [START_REF] Zinn | A new actuation approach for human friendly robot design[END_REF]. This principle allows to decouple the inertia of the links from that of the motor, thus reducing the reflected mass at the impact. If the use of SEA leads to concepts suffering from poor accuracy, the VSA and DM 2 lead to complex and costly actuators. Therefore, other interesting concepts have been proposed based on the use of preloaded compliant mechanisms [START_REF] Park | Safe link mechanism based on nonlinear stiffness for collision safety[END_REF][START_REF] Hyun | Variable stiffness mechanism for human-friendly robots[END_REF][START_REF] López-Martínez | Design and analysis of a flexible linkage for robot safe operation in collaborative scenarios[END_REF][START_REF] Zhang | On the design of mechanically safe robots based on spatial isotropic force modules and torque limiters[END_REF][START_REF] Seriani | Development of n-DoF Preloaded Structures for Impact Mitigation in Cobots[END_REF] or clutch mechanisms [START_REF] Lauzier | A Comparison of the Effectiveness of Design Approaches for Human-Friendly Robots[END_REF]. Such solutions have been used in multi-link serial architectures, however the experimental evaluation of their safety are often conducted on a single link mechanism.

The use of a parallel architecture to design safe robots has been little investigated in the literature [START_REF] Campa | A 2 dof continuum parallel robot for pick & place collaborative tasks[END_REF], although it allows to reduce the mass of the moving links. In previous works [START_REF] Jeanneau | R-min: a fast collaborative underactuated parallel robot for pick-and-place operations[END_REF], we introduced the R-Min robot, an underactuated parallel robot, whose design is based on the architecture of the five-bar mechanism in which two supplementary passive joints are added. The obtained robot can self-configure in case of a collision with the environment allowing to mitigate the severity of an impact with a human being. However, no experimental proof of the safety performances of this concept has been provided so far.

Assessing the safety of a collaborative robotic application requires the use of safety-related criteria. The automotive industry widely used the HIC [START_REF] Versace | A review of the severity index[END_REF] (Head Injury Criteria) to evaluate the severity of a crash accident for human beings. However, it appears that this index is not well suited to analyze the severity of an impact with a robot [START_REF] Haddadin | Safety evaluation of physical human-robot interaction via crash-testing[END_REF][START_REF] Haddadin | The role of the robot mass and velocity in physical human-robot interaction-part i: Nonconstrained blunt impacts[END_REF], since the velocity of the robot tip is much lower than that of a car. More recently, biomechanical thresholds have been furnished in the standard ISO/TS 15066 [START_REF]Robots and robotic devices -Collaborative robots[END_REF] for 29 different body parts, based on the experimental analysis of the pain tolerance of 100 healthy subjects [START_REF] Melia | Measuring mechanical pain: the refinement and standardization of pressure pain threshold measurements[END_REF], following the early works of [START_REF] Yamada | Fail-safe human/robot contact in the safety space[END_REF]. These thresholds define the maximum quasi-static (resp. transient) force and pressure which can be continuously (resp. temporary) applied by the robot on a human body part.

Multiple experimental works investigated the evaluation of the safety performances of new concepts of safe robot. Early works were carried out by the DLR institute in a crash-test centre dedicated to car certification. They used a dummy sitting in a car seat, and mounted a force sensor at the end of the DLR LWR III robot [START_REF] Haddadin | Safety evaluation of physical human-robot interaction via crash-testing[END_REF]. Such an experimental setup is obviously complex and costly to reproduce. Recently, the consortium of the European project COVR, dedicated to safety-related issues, published protocols to help end-users and system integrators to measure the force and pressure during a constrained [START_REF] Covr | Test robot arm for collision with fixed object (measurement of pressure over time)[END_REF] or unconstrained impact [START_REF] Covr | Test robot arm for collision with a movable object (measurement of pressure over time)[END_REF] (an impact is said to be constrained when the impacted human body region is clamped between the robot and the wall or any other rigid and fixed part of the workspace). This protocol involves the use of appropriate force and pressure sensors but also a dedicated mechanical system which imitates the bio-mechanical properties of the human body, including a spring to mimic the stiffness of the impacted human body part and a foam to mimic the non-linear properties of the skin. The characteristics of the spring (stiffness) and the foam (thickness and hardness) are specified in a former study [START_REF] Dguv | Kollaboriende Roboter, Planung von Anlagen mit der Funktion "Leistungs-und Kraftbegrenzung[END_REF] of the DGUV (Deutsche Gesetzliche Unfallversicherung) institute.

Contribution highlights

The aim of this paper is to demonstrate the interest of introducing underactuation and compliance in the mechanical architecture of a parallel robot in order to improve the safety of the robot. The performances of the R-Min robot are thus compared with those of the five-bar mechanism from which it is inspired. First, the workspace of the R-Min robot is compared with that of the five-bar mechanism. Therefore, the direct and inverse geometrico-static model and kinemato-static model of the R-Min robot are introduced in a forward and unified manner. The geometrico-static model is formulated as an optimization problem minimizing the potential energy of the robot, and it conducts, through a numerical analysis, to the definition of the robot's workspace. Singularity conditions are derived from the kinemato-static model allowing to determine the singularity loci in the robot's workspace. Then, the safety performances of the R-Min robot are experimentally compared with those of the five-bar mechanism. Therefore, a prototype of the R-Min robot was designed and manufactured. This prototype can be turned into a five-bar mechanism with equivalent inertia properties allowing a fair comparison. Following recent international standards, a measurement device is used to measure the maximum impact force and evaluate the severity of an impact. A set of experiments is defined in order to analyze the effect of multiple parameters on safety such as the robot's configuration at the time of impact, the position of the contact on the robot, the velocity of the robot, the stiffness of the preload bar.

The paper is constructed as follows: Section 2 presents the R-Min robot and the manufactured prototype. In Section 3, the geometrico-static and kinemato-static models of the R-Min robot are introduced and a singularity analysis is presented. The Cartesian workspace of the R-Min robot is drawn and compared with that of the five-bar mechanism. The effect of the robot's stiffness on the size of the workspace is analyzed. In Section 4, the dynamic model of the robot is presented as well as the associated computed torque controller. The tracking performances of the controller are then evaluated experimentally, on a prototype of the R-Min robot. Section 5 presents the experimental setup and the experimental protocol designed to allow the comparison of the safety performances of the R-Min robot with respect to that of the five-bar mechanism based on the measurement of the impact force. The effect of the stiffness of the preload bar on safety is analyzed. Conclusions are drawn in Section 6.

PRESENTATION OF THE R-MIN ROBOT

In this Section, the kinematic architecture of the robot is introduced. Its behavior when faced with different operating scenarios (in particular during a collision with an operator) is illustrated. The corresponding prototype is introduced for further numerical and experimental analysis.

Description of the R-Min kinematic architecture

The R-Min robot is an underactuated parallel mechanism designed as a proof of concept for pick-and-place operations in a collaborative cell. Its architecture presented on Fig. 1a is inspired from the five-bar parallel mechanism depicted on Fig. 1b. The traditional fivebar parallel mechanism is composed of two parallel legs made with four links in totality, which are linked through three passive revolute joints at points O 12 , O 22 and P and actuated by two motors located on revolute joints at points O i1 (i = {1, 2} the index of the leg). The resulting fivebar mechanism has two dofs (degrees of freedom) and is thus fully actuated and able to position its end-effector with a high rigidity.

Even if such a parallel architecture allows to reduce the moving masses in comparison with a serial architecture, a five-bar mechanism may still be harmful in case of a collision with a human because of its high stiffness. Therefore, the R-Min robot was designed so as to keep the interesting inertia properties of the five-bar mechanism while adding underactuation in the architecture in order to bring reconfiguration possibilities in case of impact. The resulting design, deeper described below, is composed of an underactuated seven-bar mechanism and a preload system to constrain the robot's configuration in the normal operating mode (Fig. 1a).

The seven-bar parallel mechanism is actuated by two motors located on revolute joints at points O i1 , and composed of passive revolute joints at points O i2 , O i3 and P (i = {1, 2}). All joint axes are aligned along the y-axis normal to the plane P 0 : (A, x, z). This mechanism has four dofs (two more than the five-bar mechanism) making it underactuated with two unconstrained dofs. At this stage, the resulting chain O 12 O 13 P O 23 O 22 would have no rigidity making such robot of little practicability. We therefore introduce a preload system whose kinematic chain is composed of two passive revolute joints at points A and P (axes aligned along the y axis), a passive prismatic joint aligned along -→ AP (thus lying in P 0 ) and a compression spring exerting a force between points A and P . This preload system plays the same role as elastic elements introduced in the design of underactuated hands [START_REF] Birglen | Underactuated robotic hands[END_REF], i.e. it drives the unconstrained dofs of the robot and kinematically constrains its configuration. This permits to tense the two distal legs O 12 O 13 P and O 22 O 23 P in order to obtain a robot stiff enough to be controlled, but less stiff than a traditional five-bar mechanism for safety purpose. In addition, since the preload system is located between the two legs, it cannot collide with an operator. This practical solution avoids adding masses on external legs (such as installing a torsion spring on passive joints O i3 of the seven-bar mechanism), providing a safer design.

In Figure 2, the configuration of the R-Min robot is illustrated in two different situations: during a pick-andplace trajectory, when no contact occurs (Fig. 2a) and Fig. 2: Representation of the configuration of the R-Min robot during a pick-and-place trajectory: (a) without contact, the R-Min robot lies in a configuration close to the five-bar mechanism due to the force exerted by the preload bar so that the end-effector is precisely positioned, (b) when a collision occurs on the distal link, the robot self-configures to mitigate the impact force. when a collision occurs on distal links (Fig. 2b). When the robot is subject to gravity effect only (including a payload at the end-effector), points O i2 , O i3 and P are close to be aligned thanks to the preload system (Fig. 2a). As a result, the configuration of the R-Min robot is close to that of the five-bar mechanism and the position of the endeffector can be precisely determined, making this solution practicable for pick-and-place operations. On the other hand, when loadings likely to appear during an impact with a human are applied on the robot (on a distal link), it encounters large internal reconfigurations (Fig. 2b). This behavior should avoid transmitting a large part of the robot's kinetic energy during impact, since a large part of the collision energy will be dissipated through the reconfiguration of the robot. This behaviour makes the R-Min robot a good candidate for safe physical interactions during pick-and-place operations, as it will be confirmed experimentally in what follows.

Presentation of the R-Min prototype

A prototype of the R-Min robot was manufactured based on the kinematic architecture previously defined, and is presented in Fig. 3. The key features of the design of this prototype are described below:

-All links are made in aluminium hollow tubes of circular cross-section, with (a) an external diameter of 30 mm, a thickness of 5 mm for proximal links 11 and 21 (Fig. 1a), and (b) an external diameter of 20 mm, a thickness of 2 mm for all other links. -Passive joints are made with bronze plain bearings instead of ball bearings, to lighten the mechanism. -The compression spring of the preload system is replaced by a set of two traction springs in order to avoid the self-locking effect of the translational guidance. -Two direct drive motors SIMOTICS S-1FL6 are used (maximum torque of 23.9 Nm when using motor torque control).

The geometric and inertia properties of the R-Min robot are provided in Table 1. Two interchangeable springs of the preload system are indicated. In the following, the robot with the default spring (k t = 111.2 N/m) is referred to as R-Min, while the robot with the stiffer spring (k t = 239.2 N/m) is referred to as R-Min + .

The R-Min prototype can be turned into a five-bar mechanism, by removing the preload system and locking the two passive joints located at O i3 . This solution allows to obtain two robots with comparable dimensions, inertia and contact geometry, and thus to obtain a fair comparison of the robots performances. The geometric and inertia properties of the obtained five-bar mechanism are

WORKSPACE ANALYSIS OF THE R-MIN ROBOT

In this section, the computation of the geometricostatic model of the R-Min robot is presented (in a more unified way than in [START_REF] Jeanneau | Geometrico-static analysis of a new collaborative parallel robot for safe physical interaction[END_REF]). The kinemato-static model1 is derived from the geometrico-static model and allows to express the singularity conditions. The workspace of the R-Min robot is computed numerically and compared with that of the five-bar mechanism.

Geometrico-and kinemato-static models of R-Min robot 3.1.1 Geometrico-static problem

The R-Min robot being underactuated, its equilibrium configuration is obtained at a minimum of the potential energy U . The robot configuration is parameterized with the joint angles q ij (with i = {1, 2} the index of the leg and j = {1, 2, 3} the index of the link) defining the angle between the links (i, j), q 4 , the direction of the preload system and ρ its length. In what follows, we will denote as q a = [q 11 q 21 ] T ∈ R 2 the vector of the active joint coordinates, q d = [q 12 q 13 q 22 q 23 ] T ∈ R 4 the vector of the passive joint coordinates (excluding the preload system) and p = [x z] T ∈ R 2 the end-effector position vector expressed in the frame F 0 .

In order to solve the geometrico-static problem, for a given value of the motor angles q a (considered to be blocked), it is necessary to find the configuration q dp = [q T d p T ] T ∈ R 6 which minimizes the robot potential energy [START_REF] Khalil | Chapter 14motion control[END_REF], under geometric constraints [START_REF] Robla-Gómez | Working together: A review on safe human-robot collaboration in industrial environments[END_REF]. Indeed, due to the closed-chain nature of the robot, the variables q a , q d and p are not independent. Four constraint equations grouped in the vector

ϕ ϕ ϕ = [ϕ ϕ ϕ T 1 ϕ ϕ ϕ T 2 ]
T can be defined:

ϕ ϕ ϕ i = 0 = ---→ OO i1 + ℓ i1 u i1 + ℓ i2 u i2 + ℓ i3 u i3 -p (1)
with i = {1, 2} the index of the leg and u ij = cos j k=1 q ik sin j k=1 q ik T the direction vector of the link (i, j).

Mathematically speaking, the geometrico-static problem can be formulated as:

q dp = argmin U (q a , q d , p) subject to ϕ ϕ ϕ(q a , q d , p) = 0 (2)
The total potential energy of the robot U can be expressed as a function of the coordinates q a , q d and p. For the sake of conciseness, details of the computation are given in Appendix A. This optimization problem can be reformulated using the augmented Lagrangian function defined as:

L = U (q a , q d , p) + ϕ ϕ ϕ T (q a , q d , p)λ (3)
where λ ∈ R 4 is a vector of Lagrange multipliers. Solutions to (2) can be found by solving the first-order stability conditions [START_REF] Nocedal | Numerical optimization[END_REF] ∇ q dp ,λ L = 0, which are formulated as follows:

h(q a , q d , p, λ) = 0 ϕ ϕ ϕ(q a , q d , p) = 0 (4)
with h(q a , q d , p, λ) = ∇ q dp U + ∇ q dp ϕ ϕ ϕλ, the gradient of the Lagrangian with respect to q dp . The system of equations ( 4) is composed of 10 equations, and contains 12 unknowns (q a , q d , p, λ). Therefore, it can be solved by fixing two of the unknowns to a desired value, in order to get either the forward geometrico-static model (FGSM) or the inverse geometrico-static model (IGSM). Since the solutions to Eq. ( 4) cannot be obtained analytically, a numerical solver for nonlinear systems of equations must be used. In this work, we used the Matlab implementation of the Levenberg-Marquardt algorithm implemented in the fsolve function.

Forward geometrico-static model By fixing the value of q a to a given value q * a , it is possible to compute a solution (q * d , p * , λ * ) ∈ R 10 to Eq. ( 4).

Inverse geometrico-static model In the same manner, by fixing the value of p to a given value p * , a solution (q * a , q * d , λ * ) ∈ R 10 to Eq. ( 4) can be found.

Second-order stability condition The computed configuration (q * a , q * d , p * , λ * ) is not necessarily a stable configuration, since it is not necessarily a minimizer of the potential energy U , but could be a maximizer or a saddle point. In order to assess the stability of a configuration satisfying the first-order condition given in Eq. ( 4), second-order conditions must be checked as detailed in [START_REF] Nocedal | Numerical optimization[END_REF]. These conditions can be formulated using the following two matrices:

• The Hessian matrix of the Lagrangian function with respect to q dp :

H(q * a , q * d , p * , λ * ) = ∇ 2 q dp q dp L(q * a , q * d , p * , λ * ) = D q dp h(q * a , q * d , p * , λ * ) (5)
• A matrix Z ∈ R 6×2 spanning the null space of the transposed gradient of the constraints:

Z ∈ Ker ∇ q dp ϕ ϕ ϕ T (q * a , q * d , p * , λ * ) (6)
Then, the configuration

(q * a , q * d , p * , λ * ) is stable if the projected (or reduced) Hessian H z = Z T HZ is positive definite, i.e. if H z ≻ 0 [29].

Kinemato-static modelling of the R-Min robot

The kinemato-static model of the robot can be obtained by differentiating the geometrico-static model ( 4) with respect to time:

D qa h qa + D p h ṗ + D q d h qd + D λ h λ = 0 D qa ϕ ϕ ϕ qa + D p ϕ ϕ ϕ ṗ + D q d ϕ ϕ ϕ qd = 0 (7)
Equation ( 7) can be simplified by left-multiplying its first line by matrix Z T defined in [START_REF] Pratt | Series elastic actuators[END_REF], noticing that D λ h = ∇ q dp ϕ ϕ ϕ. We then obtain:

A qa + P ṗ + D qd = 0 (8) with A = Z T D qa h D qa ϕ ϕ ϕ ∈ R 6×2 , P = Z T D p h D p ϕ ϕ ϕ ∈ R 6×2 and D = Z T D q d h D q d ϕ ϕ ϕ ∈ R 6×4
The system of equations ( 8) is made of 6 equations, and contains 8 unknowns ( qa , qd , ṗ). Therefore, by fixing two of the unknowns to a desired value, Eq. ( 8) can be solved, in order to get either the forward kinemato-static model (FKSM) or the inverse kinematostatic model (IKSM).

Forward kinemato-static model It allows to express the end-effector velocity ṗ and the passive joint velocities qd as a function of the motor velocities qa :

ṗ qd = -G PD -1 A qa (9)
where G PD = P D ∈ R 6×6 . The conditions of degeneracy of the matrix G PD are discussed later.

Inverse kinemato-static model It allows to express the motor velocities qa and the passive joint velocities qd as a function of the end-effector velocity ṗ:

qa qd = -G AD -1 P ṗ ( 10
)
where

G AD = A D ∈ R 6×6 .
The conditions of degeneracy of the matrix G AD are discussed in the next paragraph.

Singularity analysis

The system of equations ( 8) is analogous to that of a continuum parallel robot [START_REF] Black | Parallel continuum robots: Modeling, analysis, and actuation-based force sensing[END_REF], as detailed in [START_REF] Briot | Singularity conditions for continuum parallel robots[END_REF]. Based on this paper, several types of singularities may appear for this system of equations, among which the main ones are:

-Type 1 singularities: appear when matrix G AD is rank deficient, while the matrices A and D are full rank. These are singularities of the IKSM: In them, a motion of the motors do not lead to a motion of the end-effector in some given direction. They define the workspace boundaries. -Type 2 singularities: appear when matrix G PD is rank deficient, while the matrices P and D are full rank. These are singularities of the FKSM: For a fixed position of the motors, a small displacement of the end-effector is allowed. As shown in [START_REF] Briot | Singularity conditions for continuum parallel robots[END_REF], they define the boundaries of the stable configuration domains and lead to a loss of stiffness of the robot. -Type 3 singularities: appear when matrices G AD and G PD are rank deficient at the same time, while the matrices A, P and D are all full rank. These are a combination of both other types of singularities.

Other types of singularity may be found, as mentioned in [START_REF] Briot | Singularity conditions for continuum parallel robots[END_REF], but do not appear in the study of the prototype of the R-Min robot, therefore we do not recall them.

Robot workspace analysis

In what follows, the algorithm provided in [START_REF] Jeanneau | Geometrico-static analysis of a new collaborative parallel robot for safe physical interaction[END_REF] is used to draw the Cartesian workspace of the R-Min robot. Since, there exist multiple solutions to the geometricostatic problem of the R-Min robot (as detailed in [START_REF] Jeanneau | Geometrico-static analysis of a new collaborative parallel robot for safe physical interaction[END_REF]), the following inequalities are added in order to keep only the set of solutions corresponding to the configuration of the robot depicted in Fig. 2a:

-4π/3 < q 11 < -π/6 (11)
-5π/6 < q 21 < π/3 (12)

y < 0 (13) -π/2 < q i3 < π/2 (14) 0 < θ 1 < π (15) -π < θ 2 < 0 ( 16
)
with θ i being the angle between

----→ O i1 O i2 and ----→ O i2 O P .
The search of the workspace is performed only in half of the joint space since the robot is symmetric, therefore the additional constraint πq 21 < q 11 is considered.

The working mode depicted in Fig. 2a is suitable for multiple reasons: both distal arms can fold in case of a collision with a human, -the preload bar does not present any risk of injury since it lies between both distal arms, -it offers a large Cartesian space for the controlled point, -the weight of the payload has no (or little) effect on the position of the end-effector, as shown in [START_REF] Jeanneau | Geometrico-static analysis of a new collaborative parallel robot for safe physical interaction[END_REF].

The obtained workspace for the current prototype is represented on Fig. 4. It corresponds to the set of stable solutions to the geometrico-static problem, satisfying inequalities ( 11)- [START_REF] Jeanneau | R-min: a fast collaborative underactuated parallel robot for pick-and-place operations[END_REF]. The parameters of the prototype provided in Tab. 1 have been considered for the computation with the default setting of the preload system.

The values of the inverse condition numbers of the Figure . 5 decomposes the boundaries of the workspace of the R-Min robot: curves L 1 , L 4 and L 6 correspond to singularities (as it can be seen on Fig. 4), while L 2 and L 3 correspond to the mechanical limits of both motors (resp. q 21 < π/3 and q 11 < -π/6), L 5 corresponds to the symmetry of the robot. When considering unlimited motor torques, L 1 to L 6 thus represent the boundaries of the workspace for the R-Min robot. When considering the motor torque limitations of our prototype (23.9 Nm), a new boundary appears which depends on the stiffness of the preload system. This boundary can significantly reduce the workspace. Thus, the curve L 7 (resp. L 8 ) indicates the position for which one motor reaches the maximal value 23.9 Nm for the R-Min robot (resp. for R-Min + ).

L 1 L 2 L 3 L 4 L 5 L 6 x (m) z (m) 0 0.4 0.2 0 -0.3 -0.2 -0.1 -0.6 -0.5 -0.4 -0.7 L 7 L 8
The whole workspace of the R-Min robot is drawn 2 In the present work, we compute the condition number of matrices whose components have non-homogeneous units. This is valid as long as we intend to analyze the degeneracy of the corresponding matrices, and not to characterize the physical performance of the robot [START_REF] Merlet | Jacobian, manipulability, condition number, and accuracy of parallel robots[END_REF], which is not the case in this paper. on Fig. 6a. When no limitation of the motor torque is considered, the workspace of the R-Min robot is close to that of the five-bar mechanism corresponding to Tab. 2 (Fig. 6b). However, when considering the maximal torque that can be delivered by each motor (23.9 Nm for the prototype), the size of the workspace of the R-Min robot decreases. This workspace is represented for two interchangeable springs in the preload system, with a stiffness k t = 111.2 N/m and a free length ℓ 0 = 82.7 mm for the R-Min robot (Fig. 6c) and k t = 239.2 N/m, ℓ 0 = 102.5 mm for the R-Min + robot (Fig. 6d). Indeed, when increasing the spring stiffness, the size of the workspace decreases because motors cannot deliver enough torque to compensate the wrench applied by the preload system at the end-effector.

We may thus conclude that the spring stiffness of the preload system has a considerable impact on the size of the robot workspace: the lower the stiffness, the larger the workspace. In addition, a lower stiffness would intuitively be related to a safer robot as will be discussed later in section 5. However, since the R-Min robot is underactuated, there exist internal oscillations of the passive joints during a trajectory, whose amplitude is related to the robot stiffness. These oscillations are directly related to the positioning accuracy of the end-effector which is evaluated in the following section.

DYNAMIC MODELING AND CONTROL OF THE R-MIN ROBOT

In this section, the dynamic model of the R-Min robot is introduced. Based on this model, a modified version of a computed torque controller is proposed. The tracking performances of the implemented controller are then evaluated experimentally.

Dynamic model of the R-Min robot

The dynamic model of the R-Min robot is derived from the Lagrange equations. Since the robot has four dofs, it can be parameterized by four independent coordinates that will be denoted by the vector q g , defined as q g = [q T a q T u ] T ∈ R 4 , with q a = [q 11 q 21 ] T and q u = [q 12 q 22 ] T . The remaining coordinates defining the configuration of the robot are grouped into the vector q c = [q 13 q 23 p T ] T ∈ R 4 . The dynamic model is thus given by: with τ qa , τ qu and τ qc , the generalised torques associated to the coordinates q a , q u and q c , obtained using the equations of Lagrange

  τ 0 0   +   ∇ qa ϕ ϕ ϕ ∇ qu ϕ ϕ ϕ ∇ qc ϕ ϕ ϕ   ν =   τ qa τ qu τ qc   (17
τ q k = d dt ∂L ∂ qk T -∂L ∂q k T with k = a, u, c . L = E -U is the robot Lagrangian.
For the sake of conciseness, the expressions of the potential energy U and the kinetic energy E are detailed in appendices A and B. τ ∈ R 2 is the vector of the motor input efforts and ν ∈ R 4 is a vector of Lagrange multipliers (which should not be confused with vector λ introduced in the augmented Lagrangian formulation of the optimization problem in Section 3.1.1).

This model can also be rewritten under the form below, by removing the Lagrange multipliers ν, using the third line of eq. ( 17):

τ 0 = τ qa τ qu - ∇ qa ϕ ϕ ϕ ∇ qu ϕ ϕ ϕ (∇ qc ϕ ϕ ϕ) -1 τ qc ( 18
)
where ∇ qc ϕ ϕ ϕ ∈ R 4×4 is a square and invertible matrix. Finally, after some tedious undetailed mathematical simplifications, and by replacing the expressions of the variables in the vector q c and their time-derivatives by their expressions as a function of q g , qg , qg , the dynamic model can be obtained under the following standard form [START_REF] Khalil | Chapter 14motion control[END_REF]:

τ 0 = M 11 (q g ) M 12 (q g ) M 21 (q g ) M 22 (q g ) qa qu + c 1 (q g , qg ) c 2 (q g , qg ) (19)
where M kl ∈ R 2×2 (k, l = {1, 2}) are submatrices composing the robot inertia matrix, and c k ∈ R 2 are vectors containing the centrifugal, Coriolis, gravity and spring effects. It should be noted that friction terms in the joints could also be added in the vectors c k by using the methodology presented in [START_REF] Briot | Dynamics of Parallel Robots: From Rigid Bodies to Flexible Elements[END_REF].

Robot controller

A modified version of the Computed Torque Controller [START_REF] Khalil | Chapter 14motion control[END_REF] is applied for the control of the R-Min robot. For computing the motor torques, let us modify the equations of the dynamic model as follows. First, using the bottom of Eq. ( 19), we get:

qu = -M -1 22 M 21 qa -M -1 22 c 2 (20)
Then, introducing this expression of qu in the top of Eq. ( 19), the input torques τ can be obtained as follows:

τ = M 11 -M 12 M -1 22 M 21 qa +c 1 -M 12 M -1 22 c 2 (21)
In Equation ( 21), the terms M kl and c k depend obviously of q a , q u and their time derivatives. However, the R-Min prototype (presented in section 2) is equipped with only two sensors positioned on both motors, which is not enough to compute the full robot configuration. It was then necessary either to predict the behaviour of the robot or to make an assumption on the robot's configuration. In order to simplify the problem, and because the preload system tense the distal arms, angles q 13 , q 23 and their time derivatives are assumed to be null at any time, allowing to estimate q u and its time derivatives based on the second line of Eq. ( 7). This hypothesis is of course not true in practice, but this approximation will be compensated when controlling the input efforts τ thanks to the use of an integral term, as proposed in the following equation, derived from [START_REF] Melia | Measuring mechanical pain: the refinement and standardization of pressure pain threshold measurements[END_REF]:

τ = M 11 -M 12 M -1 22 M 21 q * a + k i (q * a -q a ) + k p (q * a -q a ) + k d ( q * a -qa ) + c 1 -M 12 M -1 22 c 2 (22)
with q * a , q * a and q * a , the desired position, velocity and acceleration, respectively, and k d , k p , k i three gains defined by [START_REF] Khalil | Chapter 14motion control[END_REF]: where, for our prototype, the damping coefficient ξ is fixed at 0.72 and the controller cut-off pulsation ω is experimentally set at 22 rad/s.

k d = (1 + 2ξ)ω (23)
k p = (1 + 2ξ)ω 2 (24)
k i = ω 3 (25)

Experimental evaluation of the controller performances

In this section, the tracking performances of the controller are evaluated experimentally on the prototype of the R-Min robot. A reference trajectory is defined in the Cartesian space. It consists of a linear, horizontal trajectory of 30 cm generated using a quintic polynomial motion law so as to reach a maximal linear velocity of the end-effector of 1.5 m/s (see Fig. 7). From this trajectory, the desired motor positions, velocities and accelerations q * a (t), q * a (t) and q * a (t) (the controller inputs) are computed using the dynamic model. It should be noted that the dynamic model was previously validated experimentally by a proper comparison of simulation and experimental results, whose details are not shown in this article.

The tracking errors on motor coordinates were mea-sured using the motor encoders. Results are presented in Fig. 8a, showing a maximal error of 3.2 deg. This value is relatively low regarding the simplicity of the used controller and seems low enough to achieve pick-and-place operations. The tracking error at the level of the endeffector was obtained from a camera 3 . Results are presented in Fig. 8b and show a maximal error of 17 mm along the x direction and 14 mm along the z direction. At the end of the motion, a residual error of 4 mm is probably due to the presence of dry friction into the passive joints and to encoder calibration errors which are not compensated. These errors appear to be reasonable for a pick-and-place application. Figure 8c shows the evolution of the passive angles q 13 and q 23 along the trajectory.

A maximal amplitude of 9.3 • is reached which indicates that the amplitude of these uncontrolled internal motions remains low enough not to represent an additional risk for human operators.

EXPERIMENTAL ANALYSIS OF THE SAFETY PERFORMANCES OF THE R-MIN ROBOT

In this section, we first introduce the experimental setup that will be used in order to analyze experimentally the safety performances of the R-Min robot and to compare them with those of the five-bar mechanism. Following the recommendations given in the standard ISO 15066 [START_REF]Robots and robotic devices -Collaborative robots[END_REF], the maximum impact force is chosen as the safety criteria. This impact force is measured using a dedicated Impact Force Measurement Device (IFMD), that will be presented below. This IFMD allows to mimic the bio-mechanical properties of different parts of the human body, in case of a constrained impact (i.e. when the human body is clamped between the robot and a wall or the ground). Multiple impact scenarios are conducted in order to evaluate the effects of different parameters on safety such as the configuration of the robot, the location of the impact on the robot (at the end-effector or on middle of the distal arm), the velocity of the end-effector at the time of impact, the impacted human body part or the stiffness of the preload spring.

Description of the experimental setup

Description of the compared robots

The experimental setup is constituted of two robots: a prototype of R-Min presented in Sec. 2.2 and a five-bar mechanism. The five-bar mechanism is obtained by lock- the positioning error at the end-effector and (c) the passive angles q 13 and q 23 obtained during the reference trajectory with the R-Min robot (k t = 111.2 N/m).

ing the passive joints at O 13 and O 23 of the R-Min robot and by removing the preload bar. This permits to obtain two robots with comparable geometric and inertia properties. The R-Min robot is controlled using the strategy given in Sec. 4.2, whereas the five-bar mechanism is controlled using a standard Computed Torque Control [START_REF] Khalil | Chapter 14motion control[END_REF]. Additionally, in the controller, an impact detection was settled to avoid the application of a continuous force on the impacted body after the collision. Once an impact is detected, the controller generates motor torques that compensate gravity effects only. As mentioned in the previ-ous section, a high-frequency camera allows to measure the instant configuration of the robot and obtain its timederivatives numerically.

Impact Force Measurement Device (IFMD)

The impact force is measured using a dedicated Impact Force Measurement Device (IFMD) which allows to mimic the bio-mechanical properties of different body parts such as the stiffness of the impacted body part and the non-linear properties of the skin, in case of a constrained impact. In accordance with the protocols given in [START_REF] Covr | Test robot arm for collision with a movable object (measurement of pressure over time)[END_REF][START_REF] Covr | Test robot arm for collision with fixed object (measurement of pressure over time)[END_REF], this system is composed of a translating impactor, a changeable spring, a changeable soft material cover (foam) and a force sensor. The stiffness of the spring and the hardness and thickness of the foam are specified for 29 body parts in a document of the DGUV [START_REF] Dguv | Kollaboriende Roboter, Planung von Anlagen mit der Funktion "Leistungs-und Kraftbegrenzung[END_REF]. These protocols [START_REF] Covr | Test robot arm for collision with a movable object (measurement of pressure over time)[END_REF][START_REF] Covr | Test robot arm for collision with fixed object (measurement of pressure over time)[END_REF] also recommend measuring the pressure during the collision, however, since the current paper focuses on the effect of the robot's architecture on safety and not on the local geometry of the robot where the contact occurs, the pressure was not measured. The clamping force, which is another important safety criteria [START_REF]Robots and robotic devices -Collaborative robots[END_REF] when considering a constrained impact, is not analyzed either in our evaluation. Indeed, thanks to the impact detection implemented in the controller and the subsequent reaction applied on motor torques, the clamping force is systematically null.

For the needs of the experiments, we designed our own IFMD (see Figs. 9a and9b) following the recommendations given in [START_REF] Dguv | Kollaboriende Roboter, Planung von Anlagen mit der Funktion "Leistungs-und Kraftbegrenzung[END_REF][START_REF] Covr | Test robot arm for collision with a movable object (measurement of pressure over time)[END_REF][START_REF] Covr | Test robot arm for collision with fixed object (measurement of pressure over time)[END_REF]. It is constituted of the following parts:

an impactor of mass m p linked to the ground through a prismatic joint, -a set of three foams with different Shore A hardness and thickness ℓ f (this set was provided by the Pilz company), -a set of three springs with different stiffnesses k h , -a force sensor Delta IP 68, with a measurement range of 1980 N along the impact direction, a resolution of 0.5 N, a data acquisition frequency of 7 kHz. The signal is filtered with a 4 th order low-pass butterworth filter with a 3 dB frequency of 1 kHz, -a stiff positioning table to position the force sensor with respect to the chassis of the robot.

The foam and the spring were chosen to match the properties of three different body parts: the belly, the leg and the head, as displayed on Tab. 3. These three body parts were selected because they represent the range of body stiffnesses given in the standards ISO/TS15066 [START_REF]Robots and robotic devices -Collaborative robots[END_REF] from 10 to 150 N/mm and use the three possible foams. As shown in Tab. 3, the stiffnesses of the springs selected within the manufacturer's catalogue differ from the recommended values in [START_REF]Robots and robotic devices -Collaborative robots[END_REF] for multiple reasons related to the mechanical design. However, it is not our objective in this article to assess the safety of the R-Min robot, but only to compare its safety performances with respect to the five-bar mechanism.

Experimental protocol

In the following, we explain how we designed our experiment by detailing the impact scenarios that were Table 3: Parameters of the IFMD used in the experiments, k h is the stiffness given by the manufacturer, whereas k rec h is the stiffness recommended in the standard [START_REF]Robots and robotic devices -Collaborative robots[END_REF]. selected for testing. We therefore introduce four reference impact configurations. From these impact configurations, multiple impact scenarios are defined in order to evaluate the effect of the following parameters on impact forces: the end-effector velocity, the impacted body part and the stiffness of the R-Min robot.

Body

Definition of tested Impact Configurations (IC)

We introduce here four reference impact configurations defined by the configuration of the robot at the time of impact, the direction of the end-effector's velocity and the location of the impact on the robot's arm. These impact configurations are presented on Fig. 10 and can be summarized as follows:

-Configuration "Co" corresponds to an impact on the distal arm (in the middle of link (2, 2) for the five-bar mechanism, on O 23 for R-Min), while the robot is in a centred configuration (Fig. 10a), -Configuration "Cp" corresponds to an impact on the end-effector P , while the robot is in a centred configuration (Fig. 10b), -Configuration "So" corresponds to an impact on the distal arm (in the middle of link (2, 2) for the five-bar mechanism, on O 23 for R-Min), while the robot is in a side-configuration (Fig. 10c), -Configuration "Sp" corresponds to an impact on the end-effector P , while the robot is in a sideconfiguration (Fig. 10d).

We did not consider any impact with the proximal arm (link (i, 1)), since we assume that a proper protective housing could be designed to avoid any contact with the proximal arm of the robot.

In the collision protocols furnished by the COVR consortium [START_REF] Covr | Test robot arm for collision with a movable object (measurement of pressure over time)[END_REF][START_REF] Covr | Test robot arm for collision with fixed object (measurement of pressure over time)[END_REF], it is recommended to position the IFMD such that the foam surface is normal to the direction of the velocity of the contact point attached to the robot v i . However, this protocol is well-adapted for serial robots, when impact occurs on the end of the arm. In our case, due to the parallel architecture of the robot and Table 4: Set of parameters defining the four reference Impact Configurations (IC). The orientation of the impact velocity θ i is computed from the quasi-static configuration and may vary slightly for the R-Min robot due to internal dynamics of passive joints. its bulky nature, the IFMD was positioned normal to the distal arm at the time of impact (see Fig. 10), so that the contact is punctual and takes place on the foam surface and not on its sides. Table 4 regroups the main parameters defining each four impact configurations, i.e. the end-effector's coordinates at the time of impact, the orientation of the endeffector's velocity vector θ e , the orientation of the IFMD θ c and the orientation θ i of the velocity vector of the contact point. This latter angle is not measured but computed using the geometrico-static model given in Sec. 3.

IC x(t c ) z(t c ) θ e θ c θ i m m

Impact scenarios

The design of experiments is based on multiple impact scenarios defined by the selected impacting robot (R-Min, R-Min + or five-bar mechanism), the impacted body part (belly, leg or skull head) and the selected impact configuration (as defined in Tab 4 and Fig. 10). Impact scenarios have been selected in order to allow the comparison of the R-Min robot with the five-bar mechanism from a safety perspective and to evaluate the effect of the following parameters on impact forces: the configuration of the robot, the location of the impact on the robot's arm, the impacted body part (belly, leg or head), the stiffness of the preload bar and the end-effector's velocity. Therefore, we propose the following design of experiments:

-Eight impact scenarios are defined to analyze the effect of the impact configuration on impact forces. Both robots (R-Min and the five-bar mechanism) are tested impacting the IFMD set to mimic the head. The four impact configurations are considered: Co, Cp, So, Sp. -Eight additional impact scenarios are defined to analyze the effect of the impacted body part on impact forces. Both robots are tested impacting the IFMD in two additional settings (belly and leg). Two impact configurations are considered: Co, Cp. -Two additional impact scenarios are defined to analyze the effect of the stiffness of the spring of the preload bar on impact forces. We therefore con-sider the R-Min + version of the R-Min robot which uses a set of springs with a larger stiffness (k t = 231.9 N/m). This robot is thus tested impacting the IFMD in two additional settings (belly and leg). Two impact configurations are considered: So, Sp. Impact configurations Co and Cp could be not tested for R-Min + , because the corresponding robot configurations are not inside the workspace due to the motor torque limitations. -The 18 previously defined impact scenarios are realized imposing three different end-effector's velocities at the time of impact v e (t c ) = {1, 1.5, 2} m/s in order to analyze its effect on impact forces.

In order to ensure the repeatability of the obtained measures, each impact scenario was repeated 5 times for a given end-effector's velocity. Experimental results are obtained with a good reproducibility, indeed we computed the standard deviation of the maximal impact force for each impact scenario and obtained an average standard deviation of 1.7 % and a maximum one of 4.0 % among all impact scenarios.

Experimental results and analysis

In this section, we present the experimental results obtained from the previously presented design of experiments. The time profile of the impact force was acquired for each of the 18 impact scenarios and each of the three end-effector's velocities. These results allow to analyze the effect of various parameters on safety as detailed in the following.

Comparison of R-Min with the five-bar mechanism

In this section, we present the impact forces obtained using each of the two robots impacting with the IFMD set to mimic the head with an end-effector velocity v e (t c ) = 1 m/s. All four impact configurations are considered (Co, Cp, So, Sp).

As can be seen on Fig. 11a, the maximum impact force obtained with the R-Min robot (56 N) is much lower than that obtained with the five-bar mechanism (226 N), when impact occurs on the distal arm with the robot being in a centred configuration (impact configuration Co). In this case, the maximum impact force is divided by a mean ratio 4 of r = 4.0. It should be noted that the impact duration is also significantly lower for the R-Min robot, showing that the mechanical energy transmitted to the human body part through the impact is much lower when using the R-Min robot.

However, when considering an impact at the level of the end-effector (Cp), the maximum impact forces are very similar for R-Min and the five-bar mechanism (280 N and 272 N respectively). In this case, the obtained mean ratio is r = 0.94. Similar results are obtained when considering the robot in a side-configuration (So, Sp) as depicted on Fig. 11b. We then obtain a mean ratio r = 4.6 when impact occurs on the distal arm (So) and r = 1.01 when impact occurs on the end-effector (Sp).

These results show that the underactuated architecture of the R-Min robot allows to significantly reduce the impact force when impact occurs on the middle of the distal arm at O 23 , whereas the impact force can increase slightly when impact occurs on the end-effector. Indeed, as shown on Fig. 12a, the robot stiffness k r 5 in the direction of impact of the R-Min robot is significantly lower than that of the five-bar mechanism, when considering an impact at point O 23 for both centred and sided configurations. Interestingly, in this case, the stiffness of the R-Min robot is not isotropic but is particularly low in the direction of a probable impact (normal to the distal links 22 and 23), in contrast with the five-bar mechanism. Conversely, the R-Min stiffness is rather isotropic at the endeffector but larger than the stiffness of the five-bar mech-anism (Fig. 12b).

This behaviour of the R-Min robot can be explained by the fact that, at the instant before collision, the distal links (i, 2) and (i, 3) are nearly aligned due to the tension applied in the preload bar, making the robot as stiff as the five-bar mechanism at the end-effector, independently of the stiffness of the preload bar. Figure 13b shows that the distal links of the R-Min robot remain aligned after an impact at the end-effector, in contrast with an impact at O 23 as shown in Fig. 13a. For both robots, the stiffness at the end-effector mainly relies on the controller stiffness. Indeed, the stiffness of the robot mechanics at the end-effector is infinite when considering rigid links and joints, for both the five-bar mechanism and R-Min (considering aligned distal links in this latter case). The controller stiffness is given by the proportional gain k p defined in Eq. ( 22), which depends on the inertia of the robot. This explains why the R-Min robot, whose inertia is larger than that of the five-bar mechanism, has a larger stiffness ellipsoid.

The introduction of a mechanical stop between links i2 and i3 on both arms, hindering the alignment of these two links, could help reducing the stiffness and the impact force at the end-effector, but, conversely, would increase the stiffness and the impact force at point O i3 . Authors believe a compromise could be found in the definition of this mechanical stop in order to obtain a safe behaviour at both locations, but did not investigate this idea in the present article.

In the following two sections, we check that these results remain true for different impacted body parts and different values of the end-effector velocity.

Effect of the impacted body on safety

We now investigate the effect of the impacted body part on the impact force. Therefore, three different settings of the IFMD are considered, in order to respectively mimic a constrained impact with the belly, the leg and the head (see Tab. 3). Both robots are considered in a centred configuration, impacting on distal arm (Co), with an end-effector velocity v e (t c ) = 1 m/s.

The corresponding impact force time profiles are presented on Fig. 14 and show that the maximum impact force is significantly lower for the R-Min robot than for the five-bar mechanism for each of the three body parts. A decreasing stiffness of the IFMD induces a decreasing value of the maximum impact force for both R-Min and the five-bar mechanism. A ratio r = 5.2 is obtained for an impact on the belly, r = 4.2 on the leg and r = 4.0 on the head, showing that the R-Min has a very positive effect on safety in each of these three cases. Considering the impact configuration Cp, a ratio slightly lower than one is obtained in each case, as observed in the previous section.

To compare with, the admissible transient impact forces given in the standard ISO/TS 15066 [START_REF]Robots and robotic devices -Collaborative robots[END_REF] are, respectively, 260 N and 220 N for the leg and the belly. No value is given for the head since it is a critical zone with whom any contact must be avoided. Results depicted on Fig. 14 show that the R-Min allows to obtain a maximum or the five-bar mechanism (dashed lines), with different body parts (belly, leg and head). The impact configuration Co is considered with v e (t c ) = 1 m/s. impact force that is significantly lower that the allowable thresholds for the leg and the belly and, thus, would allow operating at higher velocity, whereas the five-bar mechanism is already close to the maximum value, especially for the leg. However, this conclusion is not valid when impact occurs on the end-effector.

Effect of the end-effector velocity on safety

As mentioned in multiple works [START_REF] Lauzier | A Comparison of the Effectiveness of Design Approaches for Human-Friendly Robots[END_REF][START_REF] Herbster | A new conversion method to evaluate the hazard potential of collaborative robots in free collisions[END_REF], there exists a linear relationship between the maximum impact force and the velocity of the impact point at the time of impact v i (t c ). In what follows, we show that this linear relation is verified for each of the 18 impact scenarios.

Indeed, a linear regression going through the origin was performed for each impact scenario, considering the three end-effector's velocities v e (t c ) = {1, 1.5, 2} m/s (see Fig. 15). It is recalled that each impact scenario is led five times for each end-effector's velocity, providing 15 points for the linear regression. We obtain a maximum error on the maximum impact force of 12.0 % using the linear regression model, among all impact scenarios, with a mean standard deviation of 3.7 %. Once this linear relation has been determined for each scenario, it would be possible to determine, by extrapolation, the end-effector velocity, knowing the maximum admissible impact force defined in the ISO/TS 15066 standard [START_REF]Robots and robotic devices -Collaborative robots[END_REF].

Effect of the stiffness of the preload bar on safety

It was shown in Sec 3.2 that the stiffness of the preload bar influences the size of the workspace when considering motor torque limitations (Fig. 6). In this section, we analyze how it affects the safety performances of the R-Min robot. The impact force is measured for two different settings of the preload bar, for R-Min k t = 111.2 N/m and ℓ 0 = 82.7 mm (see Fig. 6c) for R-Min + k t = 239.2 N/m, ℓ 0 = 102.5 mm (see Fig. 6d).

Results are presented on Fig. 16 considering impact configurations So and Sp, v e (t c ) = 1 m/s and an impact with the head. These results show that the stiffness of the preload system has a very negligible effect on the impact force, whatever the location of the impact on the robot. When impact occurs on O i3 , a ratio r = 5.4 (resp. r = 4.6) is obtained with R-Min + (resp. R-Min), showing that the maximum impact force is even decreased when using a stiffer spring in the preload bar.

CONCLUSION

This paper analyzed the safety performances of the R-Min robot, a planar underactuated parallel robot designed for intrinsically safe collaborations with human. After having introduced the mechanical architecture and presented its prototype, the workspace of the R-Min robot was computed numerically based on the geometricostatic, the kinemato-static models and a singularity analysis. The obtained workspace is very close to that of an equivalent five-bar mechanism. However, the R-Min workspace decreases as the stiffness of the spring in the preload bar increases due to practical motor torque limitations. The dynamic model of the R-Min robot was then introduced as well as the modified computed torque controller. An experimental analysis showed good tracking performances despite the underactuated nature of this robot.

In order to compare the safety performances of the R-Min robot with those of the five-bar mechanism, we set up a dedicated experimental setup, i.e. a five-bar mechanism with geometric and inertia properties equivalent to that of the R-Min robot, and a bio-fidelic Impact Force Measurement Device (IFMD). Following a detailed protocol, experimental results show that the underactuated architecture of the R-Min robot allows to reduce by a ratio of 4 to 5 the impact force with respect to the five-bar mechanism when impact occurs on the middle of the distal arm. In the meantime, when impact occurs on the end-effector, the impact force is similar to that obtained with the five-bar mechanism. The effect of the stiffness of the preload bar was shown to have a negligible effect on the impact force. We also showed that there exists a linear relation between the impact force and the end-effector velocity each scenario, allowing to determine, by extrapolation, the maximum end-effector velocity from the allowable impact force defined in the ISO 15066 standard.

The proof of concept R-Min allowed us to obtain encouraging results showing the interest of introducing underactuation in the mechanical architecture of a parallel robot to obtain an intrinsically safer robot. Therefore, further works will investigate strategies to reduce the impact force when impact occurs both on the distal arm or on the end-effector. We will also extend this planar architecture to the spatial case. Simulations tools will be used to optimize the design of such robot with respect to safety criteria. In parallel, we will investigate the use of flexible links in the mechanical design of a parallel robot to obtain a highly underactuated system.

A EXPRESSION OF THE POTENTIAL ENERGY

The potential energy of the link (i, j) due to the gravity field, denoted as U i,j , depends on the links masses m ij (respectively m 4 and m 5 for links (4) and ( 5)), their centre of mass position S ij (respectively S 4 and S 5 ) and their length ℓ ij . The length ℓ Sij is equal to ℓ Oij Sij . ℓ S4 is the length from P to S 4 , ℓ S5 is the length from A to S 5 .

It is then possible to express the potential energy due to the gravity for each link:

U i,1 = gm i1 ℓ Si1 sin(q i1 ), (26)
U i,2 = gm i2 (ℓ i1 sin(q i1 ) + ℓ Si2 sin(q i1 + q i2 )) (

U i,3 = gm i3 (ℓ i1 sin(q i1 ) + ℓ i2 sin(q i1 + q i2 ) +ℓ Si3 sin(q i1 + q i2 + q i3 ))

U 4 = gm 4 z (1 -ℓ S4 /ρ) (28)
U 5 = gm 5 ℓ S5 sin(q 4 ) = gm 5 ℓ S5 z/ρ [START_REF] Black | Parallel continuum robots: Modeling, analysis, and actuation-based force sensing[END_REF] where the angle q 4 = tan -1 (z/x) and the length ρ = √ x 2 + z 2 parameterize the configuration of the preload system.

The potential energy of the robot due to the gravity can then be summarized as: U ℓ (q a , q d , p) = 2 i=1 3 j=1 U i,j (q a , q d ) + U 4 (p) + U 5 (p) [START_REF] Briot | Singularity conditions for continuum parallel robots[END_REF] The potential energy of the spring included in the preload system can be expressed as follows:

U s = 1 2 k t (ℓ 4 -ρ -ℓ 0 ) 2 (32)
where k t and l 0 are respectively the stiffness and the free length of the spring and ℓ 4 is the length of the preload bar from the spring to the end-effector position. As a result, the total potential energy of the robot U can be expressed as a function of the coordinates q a , q d and p by: U (q a , q d , p) = U ℓ (q a , q d , p) + U s (p)p c f [START_REF] Khalil | Chapter 14motion control[END_REF] where f is an external force, assumed to be conservative, applied at a contact point p c of the robot.

B EXPRESSION OF THE KINETIC ENERGY

The kinetic energy of the link (i, j) is denoted E i,j and obtained as follows:

E i,1 = 1 2 J i1 + m i1 ℓ 2 Si1 q2 i1 (34)
E i,2 = 1 2 J i2 + m i2 ℓ 2 Si2 ( qi1 + qi2 ) 2 + m i2 ℓ 2 i1 q2 i1 + 2m i2 ℓ i1 ℓ Si2 qi1 ( qi1 + qi2 ) cos(q i2 ) (35)
E i,3 = 1 2 J i3 + m i3 ℓ 2 Si3 ( qi1 + qi2 + qi3 ) 2 + m i3 ℓ 2 i1 q2 i1 + m i3 ℓ 2 i2 ( qi1 + qi2 ) 2
+ 2m i3 ℓ i1 ℓ i2 qi1 ( qi1 + qi2 ) cos(q i2 ) + ℓ Si3 ( qi1 + qi2 + qi3 ) ℓ i1 qi1 cos(q i2 + q i3 ) + ℓ i2 ( qi1 + qi2 ) cos(q i3 ) (36)

E 4 = 1 2 J 4 + m 4 ℓ S4 2 q2 4 + m 4 ẋ2 + m 4 ẏ2
-2ℓ S4 q4 (cos(q 4 ) ẋ + sin(q 4 ) ẏ)

(37)

E 5 = 1 2 J 5 + m 5 ℓ S5 2 q2 4 ( 38
)
where J ij (J 4 , J 5 , resp.) is the moment of inertia around y of the body (i, j) (4, 5, resp.). As a result, the total robot kinetic energy E can be expressed as a function of the coordinates q a , q d and p, and their time derivatives, by: E(q a , q d , p, qa , qd , ṗ) = 2 i=1 3 j=1 E i,j (q a , q d , qa , qd ) + E 4 (p, ṗ) + E 5 (p, ṗ) (39)

C EXPRESSION OF THE ROBOT STIFFNESS

Let us recall the static equilibrium equation presented on Eq. ( 4) : ∇ q dp U (q a , q d , p) + ∇ q dp ϕ ϕ ϕ(q a , q d , p)λ = τ e ϕ ϕ ϕ(q a , q d , p) = 0 (40) where τ e represents an external force applied on our mechanism, considered constant.

From equation [START_REF] Yamada | Fail-safe human/robot contact in the safety space[END_REF], in static, we may consider that the controller acts like a spring system of stiffness matrix K p = (M 11 -M 12 M -1 22 M 21 )k p , with a potential energy: U a = 1 /2(q aq v ) T K p (q aq v ), q v being the motor desired reference. The equilibrium equation can be expressed using this additional spring:

              
∇ q dp U (q a , q d , p) ∇ qa U (q a , q d , p) + K p (q aq v ) + ∇ q dp ϕ ϕ ϕ(q a , q d , p) ∇ qa ϕ ϕ ϕ(q a , q d , p) λ = τ e 0 ϕ ϕ ϕ(q a , q d , p) = 0 (41) Differentiating Eq. (41), considering a small movement of the robot around this equilibrium configuration, due to an effort dτ applied on the mechanism, we obtain: Hdq t + ∇ qt ϕ ϕ ϕdλ = dτ ∇ qt ϕ ϕ ϕ(q a , q d , p)dq t = 0 (42) with H being the Hessian of the Lagrangian L = U (q a , q d , p) + U a (q a , q v ) + ϕ ϕ ϕ T (q a , q d , p)λ with respect to variables q t = q dp q a T . The second part of Eq. (42) signifies dq t lies in the kernel of ∇ϕ ϕ ϕ. Let us define Z ∈ R 8×2 a basis of the kernel of ∇ qt ϕ ϕ ϕ T , dq t can be written under the form dq t = Zv with v ∈ R 2 .

Left-multiplying the first line of Eq. ( 42) by Z T gives:

Z T HZv = Z T dτ (43)
In the absence of Type 2 singularities, the projected Hessian H r = Z T HZ is invertible. Left-multiplying Eq. ( 43) by ZH -1 r gives:

dq t = ZH -1 r Z T dτ = Sdτ (44)
with S, being the softness matrix of the manipulator. Considering a ponctual effort df c applied at point p c applied on the mechanism, we have:

dτ = J T c df c ( 45
)
dp c = J c q t ( 46
)
with J c the contact Jacobian matrix.

We finally obtain:

dp c = J c SJ T c df c ( 47
)
where JSJ T is the softness matrix at the point of contact and K c = JSJ T -1 , the stiffness matrix. Projecting Eq. (47) along a direction n, we obtain:

dp c = n T J c SJ T c ndf c = df c k r (48)
with k r , the robot stiffness along direction n.

Fig. 1 :

 1 Fig. 1: Kinematic chains of (a) the R-Min robot and (b) the five-bar mechanism.

Fig. 3 :

 3 Fig. 3: Prototype of the R-Min robot. This prototype can be turned into an equivalent five-bar mechanism by removing the preload system and locking the two passive joints located at point O i3 (i = {1, 2}).

3 3 • 10 - 3 1 • 10 - 2 3 • 10 - 2 11 • 10 - 3 3 • 10 - 3 1Fig. 4 :

 31031021021031034 Fig. 4: Inverse condition numbers of matrices (a) G P D and (b) G AD plotted in the Cartesian workspace of the R-Min robot.

Fig. 5 :

 5 Fig. 5: Decomposition of the workspace boundaries of the R-Min robot.

x

  R-Min workspace without motor torque limitations (b) Workspace of the five-bar mechanism 0 R-Min + workspace with motor torque limitations

Fig. 6 :

 6 Fig. 6: (a-b) Representation of the Cartesian workspace of the R-Min robot for two different values of the spring stiffness in the preload bar considering a maximum motor torque of 23.9 Nm, compared with (c) the Cartesian workspace of the five-bar mechanism. The robot parameters given Tab. 1 have been considered.

  Reference trajectory for the robot end-effector (in green) given in the R-Min workspace with initial and median configuration obtained from the inverse geometrico-static model.

Fig. 7 :

 7 Fig. 7: Trajectory used for the evaluation of the controller performances.

Fig. 8 :

 8 Fig.8: Evolution of (a) the controller tracking errors, (b) the positioning error at the end-effector and (c) the passive angles q 13 and q 23 obtained during the reference trajectory with the R-Min robot (k t = 111.2 N/m).

Fig. 9 :

 9 Fig. 9: Picture of the designed Impact Force Measurement Device (IFMD) mounted on the positioning table (a) and its schematic (b).

Fig. 10 :

 10 Fig. 10: Representation of the four reference impact configurations (a) Co, (b) So, (c) Cp and (d) Sp. The trajectory is represented in green, the robot is drawn in the initial configuration (grey) and at the time of impact (black). The workspace is drawn for two different cases: (a-b) R-Min and (c-d) R-Min+. The velocity of the end-effector and that of the impact point, as well as the location of the sensor are shown for each of the four reference impact configurations.

Fig. 11 :

 11 Fig. 11: Representation of the impact force time profile due to an impact between the R-Min robot or the five-bar mechanism and the head, with v e (t c ) = 1 m/s. The four impact configurations are considered: (a) Co and Cp, (b) So and Sp.

Fig. 12 :

 12 Fig. 12: Representation of the robot stiffness ellipsoid of the R-Min robot and the five-bar mechanism, considering (a) an impact at the end-effector and (b) an impact at point O 23 .

Fig. 13 :Fig. 14 :

 1314 Fig. 13: Captures of the R-Min robot 100 ms after the collision, considering impact configurations (a) Co and (b) Cp, with v e (t c ) = 1 m/s.

Fig. 15 :

 15 Fig.15: Representation of the maximum impact force as a function of the velocity of the impact point considering the R-Min robot impacting with the head in all reference impact configurations. A linear regression is performed for each impact configuration and shows a linear behaviour of the impact force with respect to the velocity of the impact point.

Fig. 16 :

 16 Fig. 16: Representation of the impact force time profile for two different stiffnesses of the preload bar: (a) R-Min and (b) R-Min + . Scenarios So and Sp are considered, with v e (t c ) = 1 m/s and an impact with the head.

Table 1 :

 1 Geometric, inertia properties and design parameters of the prototype of the R-Min robot. Center of masses S ij (Fig.1a) are located on the lines O ij O ij+1 .x i,j is the distance from S ij to O ij .

	Link	Length Mass Inertia a	COM b
	(i; j)	ℓ i,j	m i,j	J i,j	x i,j
		(m)	(kg) (kg.m 2 )	(m)
	(1&2;1)	0.28	6.92	0.0348	0.0148
	(1&2;2)	0.2	0.302	0.007	0.126
	(1&2;3)	0.2	0.250 0.0037	0.086
	(-;4)	0.87	0.54	0.173	-0.429
	(-;5)	-	0.81	0.004	0
	ℓ O11O21 = 0.25 m
	Parameters of the preload bar spring:
	Default setting: k t = 111.2 N/m, ℓ 0 = 82.7 mm
	R-Min + setting: k t = 239.2 N/m, ℓ 0 = 102.5 mm

a Moments of inertia around y are provided at the COM. b

Table 2 :

 2 Geometric, inertia properties and design parameters of the prototype of the equivalent five-bar mechanism.

	Link	Length Mass Inertia COM
	(i; j)	ℓ i,j	m i,j	J i,j	x i,j
		(m)	(kg) (kg.m 2 )	(m)
	(1&2;1)	0.28	6.92	0.0348 0.0148
	(1&2;2)	0.40	0.59	0.0348	0.213
	ℓ O11O21 = 0.25 m
	provided in Tab. 2.

Following[START_REF] Quennouelle | Kinematostatic modeling of compliant parallel mechanisms[END_REF], we prefer to replace the word kinetostatic by the word kinemato-static. Indeed, the former is an assembly of the words kinetics and statics, and is not related with our present interest in kinematics, i.e. with the study of the motion.

A high-frequency camera is used to obtain, a posteriori, the joint and Cartesian coordinates of the robot. The camera used is a Mikrotron MC4082, images were acquired at a frequency of 500 Hz with a resolution of 1500 dpi. The camera is not used in the control of the robot.

The ratio r is the division of the maximum impact force obtained with five-bar mechanism divided by the maximum impact force obtained with the R-Min robot. The ratio values given in this article correspond to the mean value obtained from five repetitions of an impact scenario.

Appendix C details how the robot stiffness is computed.


Marco Grimaldo

Felix Roosen-Runge

Marcus Hennig

Fabio Zanini

Fajun Zhang

Michaela Zamponi

Niina Jalarvo

Frank Schreiber

Tilo Seydel
email: seydel@ill.eu.

Salt-Induced Universal Slowing Down of the Short-Time Self-Diffusion of a Globular Protein in Aqueous Solution

come

L'archive ouverte

D iffusion is a ubiquituous process in a broad variety of condensed matter systems. [START_REF]Diffusion in Condensed Matter: Methods, Materials, Models[END_REF][START_REF] Crank | The Mathematics of Diffusion[END_REF][START_REF] Dhont | An Introduction to Dynamics of Colloids[END_REF] In the particular case of the biological cell, protein self-diffusion (i.e., the diffusion of single "tagged" protein macromolecules) governs, for example, transport processes, reaction kinetics, and information exchange in the biological cell. Importantly, self-diffusion is affected by single-particle properties as well as by environmental factors. One well-known environmental control parameter for selfdiffusion is macromolecular crowding, that is, the dense packing fraction at which macromolecules such as proteins are usually found in the aqueous intracellular environment (refs 4-6 and references therein). The motion of proteins is affected by two kinds of interactions with surrounding macromolecules: direct interactions such as excluded volume and Coulomb interaction influence the dynamics on time scales of roughly microseconds and above. On shorter time scales, for example, on the order of several nanoseconds, proteins interact only through quasiinstantaneous hydrodynamic interactions. [START_REF] Dhont | An Introduction to Dynamics of Colloids[END_REF] The presence of salts is another important control parameter that influences both structure and dynamics of protein solutions [START_REF] Roosen-Runge | Protein Diffusion in Crowded Electrolyte Solutions[END_REF][START_REF] Heinen | Viscosity and Diffusion:Crowding and Salt Effects in Protein Solutions[END_REF] as well as, for example, DNA solutions, [START_REF] Grosberg | Colloquium: The Physics of Charge Inversion in Chemical and Biological Systems[END_REF] and may lead to the formation of protein aggregates. [START_REF] Stradner | Equilibrium Cluster Formation in Concentrated Protein Solutions and Colloids[END_REF][START_REF] Soraruf | Protein Cluster Formation in Aqueous Solution in the Presence of Multivalent Metal Ions-a Light Scattering Study[END_REF] The interplay of macromolecular concentration and salts, physiological ones, and others to control condensation phenomena in protein solutions represents an important question with applications in medicine and biotechnology. Protein condensation is observed in several diseases, [START_REF] Gunton | Protein Condensation: Kinetic Pathways to Crystallization and Disease[END_REF] such as some forms of cataracts [START_REF] Benedek | Cataract as a Protein Condensation Disease: the Proctor Lecture[END_REF][START_REF] Bloemendal | Ageing and Vision: Structure, Stability and Function of Lens Crystallins[END_REF] and neurodegenerative diseases. [START_REF] Ross | Protein Aggregation and Neurodegenerative Disease[END_REF] In Parkinson's disease, raised levels of multivalent iron ions have been observed in sick neurons, [START_REF] Oakley | Individual Dopaminergic Neurons Show Raised Iron Levels in Parkinson Disease[END_REF] pointing to a yet-to-be-explored role of salt. A deep understanding of the aggregation process of proteins is also of high interest for the tunable production of nanostructures through self-assembly. [START_REF] Whitesides | Self-Assembly at all Scales[END_REF] The formation of clusters in dependence on protein concentration has been previously studied using small-angle scattering (SAS), [START_REF] Stradner | Equilibrium Cluster Formation in Concentrated Protein Solutions and Colloids[END_REF][START_REF] Stradner | A Small-angle Scattering Study on Equilibrium Clusters in Lysozyme Solutions[END_REF][START_REF] Kowalczyk | Equilibrium Clusters in Concentrated Lysozyme Protein Solutions[END_REF] nuclear magnetic resonance (NMR), [START_REF] Barhoum | NMR Detection of an Equilibrium Phase Consisting of Monomers and Clusters in Concentrated Lysozyme Solutions[END_REF] dynamic light scattering (DLS), [START_REF] Johnston | Concentrated Dispersions of Equilibrium Protein Nanoclusters that Reversibly Dissociate into Active Monomers[END_REF] and neutron spin-echo spectroscopy (NSE). [START_REF] Porcar | Formation of the Dynamic Clusters in Concentrated Lysozyme Protein Solutions[END_REF] Although studied since the pioneering work of Hofmeister more than a century ago, [START_REF] Hofmeister | Zur Lehre von der Wirkung der Salze[END_REF] the effect of salt ions on the phase behavior of protein solutions presents still numerous questions. For the particular case of acidic proteins, trivalent salts such as FeCl 3 and YCl 3 induce complex reentrant phase diagrams depending on the protein concentration, c p , and salt concentration, c s , [START_REF] Zhang | Reentrant Condensation of Proteins in Solution Induced by Multivalent Counterions[END_REF][START_REF] Zhang | Reentrant Condensation, Liquid-Liquid Phase Separation and Crystallization in Protein Solutions Induced by Multivalent Metal Ions[END_REF] including macroscopic aggregation, liquid-liquid phase separation (LLPS), and protein crystallization, [START_REF] Zhang | Reentrant Condensation, Liquid-Liquid Phase Separation and Crystallization in Protein Solutions Induced by Multivalent Metal Ions[END_REF][START_REF] Zhang | Novel Approach to Controlled Protein Crystalization Through Ligandation of Yttrium Cations[END_REF][START_REF] Zhang | Charge-Controlled Metastable Liquid-Liquid Phase Separation in Protein Solutions as a Universal Pathway Towards Crystallization[END_REF][START_REF] Zhang | The Role of Cluster Formation and Metastable Liquid-Liquid Phase Separation in Protein Crystallization[END_REF][START_REF] Sauter | Real-Time Observation of Nonclassical Protein Crystallization Kinetics[END_REF] and offer promising ways to manipulate the behavior of proteins. On the basis of the complex phase behavior, it is tempting to speculate that dynamic precursor processes of protein aggregation and crystallization occur already at lower trivalent salt concentrations. Such precursor processes may be the formation of small static or transient protein clusters. Indeed, our group has also explored the formation of clusters due to the presence of multivalent salts by SAS [START_REF] Zhang | The Role of Cluster Formation and Metastable Liquid-Liquid Phase Separation in Protein Crystallization[END_REF] and DLS. [START_REF] Soraruf | Protein Cluster Formation in Aqueous Solution in the Presence of Multivalent Metal Ions-a Light Scattering Study[END_REF] Three factors render short-time self-diffusion a promising quantity to probe clustering of proteins as controlled by salt. First, the short-time self-diffusion of proteins can be described remarkably well by the theory of diffusion for effective hard spheres, [START_REF] Roosen-Runge | Protein Self-Diffusion in Crowded Solutions[END_REF][START_REF] Grimaldo | Diffusion and Dynamics of Γ-Globulin in Crowded Aqueous Solutions[END_REF] which links experimental results to a physical picture only based on the diffusional tensor of the single particle and an effective viscosity of the environment. In particular, the shift of structural correlation peaks due to a varying cluster size does not influence the interpretation. Second, the addition of salt efficiently alters direct interactions such as Coulomb repulsion, while hydrodynamic interactions remain comparable, which allows us to study the nontrivial interplay of protein and salt concentration. Third, the short observation time scale much smaller than the collision time τ i ≈ 200 ns ensures that the effects of particle collisions and specific interactions with obstacles do not affect the observed diffusion. This short-time scale is routinely accessible by neutron spectroscopy as opposed to NMR, DLS, and fluorescence correlation (FCS) techniques.

In this context, incoherent quasi-elastic neutron backscattering (QENS) is an essential tool because it unambiguously accesses the self-diffusion of globular proteins on a nanometer length -scale and a nanosecond time scale, even for crowded and possibly opaque samples. [START_REF] Roosen-Runge | Protein Self-Diffusion in Crowded Solutions[END_REF] Employing QENS to systematically investigate the protein short-time self-diffusion can hence provide crucial information on the presence or absence of clusters with a lifetime down to nanoseconds. Spectra of protein solutions as a function of the energy transfer, ℏω, and the scattering vector, q, were recorded at the backscattering spectrometer BASIS [START_REF] Mamontov | A Time-of-flight Backscattering Spectrometer at the Spallation Neutron Source[END_REF] (SNS, ORNL) and fitted with the model described in detail elsewhere [START_REF] Roosen-Runge | Protein Self-Diffusion in Crowded Solutions[END_REF][START_REF] Stadler | Hemoglobin Dynamics in Red Blood Cells: Correlation to Body Temperature[END_REF][START_REF] Grimaldo | High-Resolution Neutron Spectroscopy on Protein Solution Samples[END_REF]

ω ω β γ ω ω β γ ω = ⊗ • + - Γ + S q q q A q q A q q q q ( , ) ( , ) { ( ) [ ( ) ( ( ), ) (1 ( )) ( ( ), ) ] ( ) ( ( ), )} 0 0 D O D O 2 2 (1)
Therein, ω q ( , ) represents the energy resolution function, and ω • ( , ) stands for Lorentzian functions with a line width related to dynamical processes in the sample. γ(q) describes the global diffusion of the proteins, consisting of both translational and rotational contributions. Γ(q) represents the internal relaxations accessed by the dynamic window of the spectrometer and convoluted with the global diffusion. β(q) and A 0 (q) are scalars. The contribution of the solvent is modeled by the fixed term

β γ ω ( , ) D O D O 2 2
, as explained elsewhere. [START_REF] Grimaldo | High-Resolution Neutron Spectroscopy on Protein Solution Samples[END_REF] To increase the contrast between the hydrogenated proteins and the surrounding solvent, we use D 2 O rather than H 2 O, as commonly done in neutron scattering by virtue of the large difference of the incoherent neutron scattering cross sections of the two hydrogen isotopes.

An example spectrum is shown in the inset of Figure 1. The red circles represent the recorded spectrum, the blue dotteddashed line depicts γ ω ( , ), and the green dashed line represents ω Γ ( , ) . Finally, the turquoise solid line super-imposed on the data is the result of the fit using eq 1. The widths γ(q) obtained from the fit function in Figure 1 are plotted in the main part of Figure 1. The solid red line denotes the fit with the Fickian law γ(q) = Dq 2 , where D is the global apparent self-diffusion coefficient of the proteins. [START_REF] Roosen-Runge | Protein Self-Diffusion in Crowded Solutions[END_REF][START_REF] Grimaldo | Diffusion and Dynamics of Γ-Globulin in Crowded Aqueous Solutions[END_REF] In this study, we explore the global diffusion coefficient D in solutions of the globular protein bovine serum albumin (BSA) depending on its concentration c p and on the concentration c s of the salt YCl 3 (cf. dots in Figure 2b). For several fixed protein volume fractions of approximately 7-20%, we increase c s while remaining on an area of the phase diagram where no macroscopic aggregation or phase transition occurs. [START_REF] Zhang | Reentrant Condensation of Proteins in Solution Induced by Multivalent Counterions[END_REF] Strikingly, when plotting the reduced self-diffusion coefficient D(c s ,c p )/D(c s = 0, c p ) as a function of the number of Y 3+ ions per protein N s = c s /c p in Figure 2a, we obtain a universal scaling for all measured protein concentrations. The master curve decreases continuously to 60% of the value with no added salt, when in the solution there are approximately six Y 3+ ions per protein. Essentially, this curve characterizes the diffusion of BSA in the presence of YCl 3 , in the entire explored range of c p and c s .

The dashed line in Figure 2a represents an empirical function for the master curve

= + + ⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ ⎛ ⎝ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ g c c a c c a c c ( / ) 1 s p 2 s p 2 4 s p 4 (2)
with the parameters a 2 = (-17.32 ± 2.76) × 10 -3 and a 4 = (18.86 ± 8.93) × 10 -5 .

Having established this rather remarkable universal behavior, we now turn to its possible origin. There may be several rather complex ways to model the salt dependence. To keep the model simple and limit the number of parameters, we consider a salt-induced cluster formation, inspired by the theory of ionactivated patchy particles, [START_REF] Wertheim | Fluids with Highly Directional Attractive Forces. I. Statistical Thermodynamics[END_REF][START_REF] Roosen-Runge | Ionactivated Attractive Patches as a Mechanism for Controlled Protein Interactions[END_REF] which was successfully used to semiquantitatively describe the phase diagram of the system of BSA and YCl 3 . In this physical picture, Y 3+ ions bind to patches on the particle surface and thereby activate them for subsequent bridging between protein molecules. Thus, increasing the ions

ρ * = - - -! - + ! ! - - p p p m m n n m n n n (1 ) [ (1 ) ] ( ) ( 2 2
)

n m m n b b b 2 1 ( 3
)
where m is the maximum number of ion-bridges per particle and n is the number of particles in the n-cluster. We chose m = 4 consistent with ref 35, inspired by the crystal structure displaying four ion bridges per monomer. [START_REF] Zhang | Novel Approach to Controlled Protein Crystalization Through Ligandation of Yttrium Cations[END_REF] To compare our data with this theory we assume that at the short times probed with BASIS the clusters are in good approximation rigid. Thus, movements of proteins within the clusters are neglected. Given the assumed strong bonds, these motions might occur on similar time scales as, for example, protein interdomain motions [START_REF] Biehl | Direct Observation of Correlated Interdomain Motion in Alcohol Dehydrogenase[END_REF] beyond the observation time scale of our experiment. Moreover, we assume that each ncluster diffuses according to the theory of colloidal hard-sphere suspensions as if they were in monodisperse suspensions at the actual volume fraction of the protein solution, that is, for every cluster size n

ϕ ϕ ϕ = = ν - D D n f ( ) ( 0) ( ) n t ( ) t (1) t (4)
ϕ ϕ ϕ = = α - D D n f ( ) ( 0) ( ) n r ( ) r (1) r (5)
where D t (n) and D r (n) are the translational and rotational diffusion coefficients of the n-cluster, D t [START_REF]Diffusion in Condensed Matter: Methods, Materials, Models[END_REF] (ϕ = 0) = k B T/ (6πηR p ), D r [START_REF]Diffusion in Condensed Matter: Methods, Materials, Models[END_REF] (ϕ = 0) = k B T/(8πηR p 3 ), and R p denotes, respectively, the dilute limit diffusion coefficients and the effective hydrodynamic radius of the monomers. η represents the viscosity of the solvent and f t,r (ϕ) is the theoretical reduced diffusion coefficient for hard-sphere suspensions. [START_REF] Roosen-Runge | Protein Self-Diffusion in Crowded Solutions[END_REF][START_REF] Tokuyama | Dynamics of Hard-Sphere Suspensions[END_REF][START_REF] Banchio | Short-Time Transport Properties in Dense Suspensions: From Neutral to Charge-Stabilized Colloidal Spheres[END_REF] Lacking a detailed description of the microscopic structure of the clusters, we choose ν = 1/3 and α = 1, valid in the limit of low fractality. Thus, the hydrodynamic radius of an n-cluster scales as that of a hard sphere of volume nV p , where V p is the volume of a monomer.

The dynamic structure factor of the global diffusion is modeled for the cluster solution by

∑ ω ρ ω = * γ S q s p nS q ( , ) ( ) ( , ) n n n b (6)
Here s is a scalar, and the well-known scattering function for a single n-cluster S n (q,ω) is provided in eq 7. Fitting of eq 6 to the experimental γ ω ( , ) obtained from eq 1 returns the two free parameters s and p b . Interestingly, p b for different protein concentrations, as derived from the master curve of the dynamics, follows itself a master curve (green line in the inset of Figure 2b). Thus, the model reflects the damping of the diffusion, provided that the ion-bridge probability p b depends only on c s /c p and not on ϕ. In other words, the observed universality of the dynamics would suggest that the contour lines of p b in the explored area on the c p -c s phase-diagram are straight lines with intercept 0 and increasing slope (cf. Figure 2b). A possible reason for the universality could be the nearly quantitative binding of cations to the protein due to a dominating ion binding energy. Such a result is important for a further refined picture along the model of ion-activated attractive patches in ref 35, also in connection with the generalized law of correspondent states for patchy interactions. [START_REF] Foffi | On the Possibility of Extending the Noro-Frenkel Generalized Law of Correspondent States to Nonisotropic Patchy Interactions[END_REF] Importantly, this cluster picture could also explain previous results from light scattering in BSA solutions with lower volume fractions. [START_REF] Soraruf | Protein Cluster Formation in Aqueous Solution in the Presence of Multivalent Metal Ions-a Light Scattering Study[END_REF] In these experiments, both reduced long-time collective diffusion and reduced isothermal compressibility showed a universal scaling with respect to c s /c p . With increasing c s , the diffusion of both a monomer and a cluster fraction showed a substantial slowing down and the compressibility diverged. Using the ion-bridge probability derived from QENS, the DLS observations could be explained by the occurrence of increasingly large clusters with a size distribution universal in protein concentration because the dynamics slows down and the total scattering power for light diverges once the clusters are large enough (cf. figure 3 in ref 35).

We emphasize that our experiment does not allow us to draw a conclusion of whether the observed clusters are dynamic or transient (i.e., have a finite lifetime) or static.

On the nanosecond time scale explored by our experiment, our model explains the measured short-time self-diffusion by at least temporarily "rigid" clusters, similar to NSE results on the formation of lysosyme clusters induced by the increase in the protein concentration (i.e., not by additional salt), where the The Journal of Physical Chemistry Letters lifetime of the observed clusters is longer than a few tens of nanoseconds. [START_REF] Porcar | Formation of the Dynamic Clusters in Concentrated Lysozyme Protein Solutions[END_REF][START_REF] Cardinaux | Cluster-driven Dynamical Arrest in Concentrated Lysozyme Solutions[END_REF] As an alternative scenario, the increase in attractive interactions between proteins due to bound Y 3+ ions on the protein surface could enhance fluctuations of the local volume fraction φ of monomers and thus decrease the averaged apparent self-diffusion coefficient without the formation of bonds between particles because most proteins would experience a denser packing. Under these assumptions, we would be able to describe our data assuming a well-separated bimodal distribution of φ; however, two points seem not to be consistent with a physical explanation along this scenario: (i) There is no obvious physical reason for a well-separated bimodal distribution G(φ) away from a phase separation. (ii) There is no obvious physical reason why density fluctuations should produce a universal slowing down because this implies a very specific relation between the fluctuation amplitude and the overall protein and salt concentration. (See Figure SI1 in the SI.)

In conclusion, we experimentally establish the existence of a remarkably universal slowing down of the short-time selfdiffusion of the protein BSA as a function of the number of multivalent Yttrium ions per protein. The observation can be connected to previous results 11 on a universal slowing down of the long-time collective diffusion of proteins in the same system (BSA + YCl 3 in solution) at lower protein concentration. In combination, the two results suggest that dynamics and thermodynamics in protein solutions can be tuned in a relatively universal manner for a broad range of protein concentration by the addition of multivalent cations.

Salt-induced cluster formation via ion bridges renders a clear physical explanation for the universality and strength of the slowing down, provided that the ion-bridge probability depends only on the number of ions per protein, as, for example, expected for nearly quantitative binding of cations to the protein surface, embedding well into a general qualitative understanding of dynamics and thermodynamics in protein solutions with multivalent cations by assuming ion-activated attractive patches. [START_REF] Roosen-Runge | Ionactivated Attractive Patches as a Mechanism for Controlled Protein Interactions[END_REF] Such a result is promising in several respects. First, similar results might also be observed in other systems with strong cross-linking. Second, if the observations were indeed to be explained with the formation of clusters distributed according to a remarkably universal cluster distribution, the finding would be extremely relevant for a rational choice of pathways for selfassembly.

■ EXPERIMENTAL METHODS

Sample Preparation. BSA and YCl 3 were purchased from Sigma-Aldrich with a purity of 98 (A3059) and 99.99% (451363), respectively. D 2 O was purchased from Acros Organics with an indicated purity of 99.8%. Samples were prepared as described in ref 6. After some hours on a roller mixer, all BSA powder was completely dissolved and the solutions appeared clear. For every protein concentration, a series of samples with increasing c s were prepared. The solutions were filled into double-walled aluminum cylinders (23 mm outer diameter, 0.15 mm gap, that is, difference between inner and outer radius), sealed against vacuum, and kept at T = 295 K for the measurements.

Quasielastic Neutron Backscattering. We used the backscattering spectrometer BASIS 31 at the ORNL's Spallation Neutron Source in Oak Ridge, Tennessee in the standard configuration with unpolished Si(111) analyzer crystals, which set the selected neutron wavelength to 6.27 Å. The energy resolution defined by the fwhm of the resolution function is ∼3.5 μeV. The accessible energy transfer range is -100 ≤ ΔE ⩽ +100 μeV. The resulting accessible time range is thus

τ ≤ = Δ ≤ h E 41 ps 1.2 ns
The q range 0.25 ≤ q ≤ 1.95 Å -1 sets the probed length scale l = 2π/q to ≥ ≥ l 25

3.22Å

The raw data were normalized to the incident neutron flux and relative detector efficiency. To remove the contribution of the sample holder from the spectra, we have used the Paalman-Pings coefficients, accounting for the q-dependent absorption of neutrons by the sample and the cell walls. The details can be found in the Supporting Information.

Determination of p b . The scattering function of an n-cluster is given by

∑ ω ω = + + + + + = ∞ S q B q D l l D q D l l D q ( , ) ( ) ( 1)
[ Therein, D t (n) and D r (n) are defined by the eqs 4 and 5, respectively. B l (q) is defined as where D t ≡ D t [START_REF]Diffusion in Condensed Matter: Methods, Materials, Models[END_REF] and D r ≡ D r [START_REF]Diffusion in Condensed Matter: Methods, Materials, Models[END_REF] , because at c s = 0 the protein solution is a monodisperse suspension of monomers (i.e., n = 1).

∫ ρ = + ∞ B q r r l j qr ( ) d ( )( 2 1) ( )

Figure 1 .

 1 Figure 1. Inset: Example backscattering spectrum I(q,ω) (red circles) recorded at BASIS for BSA and YCl 3 in D 2 O (c p = 150 mg/mL, c s = 8 mM, T = 295 K, individual detector at q = 0.65 Å -1 ). The blue dotteddashed line depictsγ ω ( , ) and the green dashed line pictures ω Γ ( , ) from eq 1. The turquoise solid line superimposed on the data is the result of the fit using eq 1. Main figure: Fitted γ (red circles) versus q 2 for the full q range of the example data. The fit γ = Dq 2 (solid red line) indicates a simple Brownian diffusive behavior with apparent self-diffusion coefficient D.

Figure 2 .

 2 Figure 2. (a) Reduced short-time apparent self-diffusion coefficient D(c s ,c p )/D(c s = 0, c p ) as a function of c s /c p (symbols). Independently from the protein concentration, all of the points lie on a master curve, which is parametrized by a polynomial function (dashed line, see text). (b) Contour plot of the bonding probability, p b , as a function of c p and c s . The dots denote the coordinates where experimental data were recorded. The area with the solid lines corresponds to the range of protein concentrations measured in this study. The dashed lines are extrapolations to lower and higher protein concentrations. Inset: Bonding probability, p b , as a function of the number of Y 3+ ions per protein c s /c p , as obtained from the experiment under the assumptions explained in the text, using eqs 2 and 6. (See the Experimental Methods section.)

  r) describes the radial density distribution of the hydrogen atoms in the molecule and j l (x) denotes the lth-order spherical Bessel function of first kind. Under the assumption of homogeneously distributed H atoms in effective hard spheres, ρ(r) = 4πr 2 . To obtain p b , one can fit eqs 6 and 7 to a Lorentzian function γ ω ( , ) with γ defined by c s /c p ) denotes the master curve and D(c s = 0, c p ) is the apparent diffusion coefficient of the solution at protein concentration c p and no added salt, which can be determined experimentally or calculated for an effective hard sphere by solving the equation[START_REF] Roosen-Runge | Protein Self-Diffusion in Crowded Solutions[END_REF]

Downloaded via INST LAUE-LANGEVIN on January 23, 2023 at 08:05:38 (UTC).See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

DOI: 10.1021/acs.jpclett.5b01073 J. Phys. Chem. Lett. 2015, 6, 2577-2582

■ ACKNOWLEDGMENTS

The research at Oak Ridge National Laboratory's Spallation Neutron Source was sponsored by the Scientific User Facilities Division, Office of Basic Energy Sciences, U.S. Department of Energy. The spectrometer BASIS is supported by the Julich Center for Neutron Science (JCNS), Germany, via the partner user programme and M.G. and F.Za. acknowledge travel funding by JCNS to conduct the experiments as well as a student grant by the Institut Laue-Langevin (ILL). M.H., M.G. and T.S. have obtained travel funding by ILL to conduct the experiments. We acknowledge E. Mamontov, M. Oettel, R. Roth, and H. Schober for fruitful discussions and R. Moody and V. Glenisson for assistance. We further acknowledge financial support by the DFG.

Present Address

# F.Za.: Max-Planck Institute for Developmental Biology, Spemannstraße 35, D-72076 Tubingen, Germany.

Notes

The authors declare no competing financial interest.