Jason Chevrie

Alexandre Krupa

Marie Babel

Real-time Teleoperation of Flexible Beveled-tip Needle Insertion using Haptic Force Feedback and 3D Ultrasound Guidance

Needle insertion procedures can greatly benefit from robotic systems to improve their accuracy and success rate. However, a fully automated system is usually not desirable and the clinicians need to be included in the control loop. In this paper we present a teleoperation framework for beveledtip flexible needle steering that enables the user to directly and intuitively control the trajectory of the needle tip via a haptic interface. The 6 degrees of freedom of the needle base are used to perform several automatic safety and targeting tasks in addition to the one controlled by the user. Real-time visual feedback is provided by a 3D ultrasound probe and used to track the 3D location of the needle and of a spherical target. Several haptic force feedback are compared as well as two different levels of mix between automated and user-controlled tasks. A validation of the framework is conducted in gelatin phantom and a mean targeting accuracy of 2.5 mm is achieved. The results show that providing an adequate haptic guidance to the user can reduce the risks of damage to the tissues while still letting the surgeon in control of the tip trajectory.

I. INTRODUCTION

Needle insertion is a widely performed medical procedure used for a large variety of applications, such as the treatment and diagnosis of cancers. Such procedures require a high accuracy in the placement of the needle tip to ensure the success of the operation and to avoid medical complications or misdiagnosis [START_REF] Abolhassani | Needle insertion into soft tissue: A survey[END_REF]. Many research work have focused on the modeling of the needle/tissue interaction to predict the deformation of the needle during the insertion [START_REF] Misra | Mechanics of flexible needles robotically steered through soft tissue[END_REF], [START_REF] Khadem | A two-body rigid/flexible model of needle steering dynamics in soft tissue[END_REF], [START_REF] Chevrie | Online prediction of needle shape deformation in moving soft tissues from visual feedback[END_REF]. Such models have then been used in robotic control frameworks to assist the surgeon and to improve the targeting accuracy [START_REF] Abayazid | Integrating deflection models and image feedback for real-time flexible needle steering[END_REF], [START_REF] Waine | Needle tracking and deflection prediction for robot-assisted needle insertion using 2d ultrasound images[END_REF], [START_REF] Adagolodjo | Inverse real-time finite element simulation for robotic control of flexible needle insertion in deformable tissues[END_REF]. Robots can indeed be of great help to perform repetitive gestures with consistency and accuracy.

Imaging modalities such as magnetic resonance imaging (MRI), computed tomography (CT) or ultrasound (US) are usually used to perform needle insertion procedures under visual guidance. In the case of MRI or CT, surgeons have a restricted access to the patient because of the bulky scanners which limit the size of the workspace or because of the ionizing radiations emitted during the CT image acquisition. On the contrary, with US the workspace is less limited and the insertion can be performed at the same time as the imaging process. However the accuracy of the needle insertion can be reduced by freehand probe holding since the surgeons have to coordinate the position of both the US probe and the needle [START_REF] Kaye | Robotic ultrasound and needle guidance for prostate cancer management: review of the contemporary literature[END_REF]. In each case, teleoperated robotic systems can offer better operating conditions for the J. Chevrie is affiliated with Univ Rennes, Inria, CNRS and IRISA, France.

jason.chevrie@gmail.com.

M. Babel is affiliated with Univ Rennes, INSA, CNRS, Inria and IRISA, France. marie.babel@irisa.fr.

A. Krupa is affiliated with Univ Rennes, Inria, CNRS and IRISA, France.

alexandre.krupa@inria.fr. surgeons, who do not need to be directly near the patient anymore to operate, but can be positioned in an ergonomic way with an optimized access to the different available intraoperative feedback modalities.

A. Related work

Many teleoperated robotic systems have been designed for needle insertion procedures [START_REF] Piccin | A force feedback teleoperated needle insertion device for percutaneous procedures[END_REF], [START_REF] Moreira | The miriam robot: A novel robotic system for mr-guided needle insertion in the prostate[END_REF]. However, the surgeon is still often required to directly control the different degrees of freedom (DOF) of the robot, which can have various and complex effects on the actual tip trajectory. In addition, teleoperated robots that are designed to work in the limited workspace of MRI or CT are often controlled in a discontinuous manner due to the non real-time nature of the imaging process. The insertion is then performed step by step and alternates between image acquisition, correction of the needle base placement and insertion steps [START_REF] Goldenberg | Robotic system for closedbore mri-guided prostatic interventions[END_REF].

In the case of a continuous insertion whose velocity control is left to the surgeon, the robot motion is often limited to the insertion and rotation around the needle axis [START_REF] Pacchierotti | Teleoperation of steerable flexible needles by combining kinesthetic and vibratory feedback[END_REF], [START_REF] Seifabadi | Toward teleoperated needle steering under continuous mri guidance for prostate percutaneous interventions[END_REF], [START_REF] Moreira | The miriam robot: A novel robotic system for mr-guided needle insertion in the prostate[END_REF]. Continuous needle insertion usually requires the use of US imaging to obtain a fast visual feedback and automatic needle insertion under 2D or 3D US guidance have been the focus of many work [START_REF] Abayazid | Integrating deflection models and image feedback for real-time flexible needle steering[END_REF], [START_REF] Adebar | 3d ultrasound-guided robotic needle steering in biological tissue[END_REF], [START_REF] Chatelain | 3d ultrasound-guided robotic steering of a flexible needle via visual servoing[END_REF], [START_REF] Mignon | Beveled-tip needle-steering using 3d ultrasound, mechanical-based kalman filter and curvilinear roi prediction[END_REF]. However there have been only few work on continuous needle steering under US guidance using all available DOF of the needle base [START_REF] Neubach | Ultrasound-guided robot for flexible needle steering[END_REF], [START_REF] Chevrie | Flexible needle steering in moving biological tissue with motion compensation using ultrasound and force feedback[END_REF] and, to the best of our knowledge, none of them included the surgeon in the control loop for an robot teleoperation.

B. Contributions

In this paper we propose a preliminary study of a control framework that enables the real-time semi-automatic teleoperation of the needle tip 3D trajectory during its insertion in soft tissues. The framework uses all 6 DOF of the needle base and intuitively gives the surgeon the control of the tip trajectory, while automatically performing an automatic orientation of the bevel and computing the optimal needle base velocity that eventually needs to be applied by the robot to achieve the desired tip motion. The insertion is performed using 3D US feedback in order to get a real-time tracking of both the needle and the target, which are considered as inputs for the controller, and to provide visual information to the surgeon. A haptic interface is used to get the input from the surgeon and to give an informative haptic force feedback. Different methods of computing the force feedback are here compared in terms of targeting performance and safety of the procedure.

This paper is divided as follows: the different constitutive parts of the control framework are presented in section II. The different strategies for the control and for the computation of the haptic feedback are also detailed. In section III we describe the experiments performed to validate the teleoperation framework and we discuss the results obtained from a comparative study of the different kinds of haptic feedback.

Finally section IV provides conclusions on this work and directions for future work and improvement.

II. SEMI AUTONOMOUS INSERTION METHODS

This section describes the different elements of our semiautonomous teleoperation framework summarized in Fig. 1.

A. Needle insertion framework

We consider the case of a needle attached by its base to a robot arm. The control of the robot is performed using the task function framework [START_REF] Siciliano | A general framework for managing multiple tasks in highly redundant robotic systems[END_REF]. Let v = [v ω] ∈ se(3) be the velocity screw vector applied to the needle base by the robot, with v ∈ R 3 the translation velocity and ω ∈ R 3 the angular velocity. We compute v such that N objective tasks can be fulfilled at the same time. Each task is described by a vector function e i ∈ R ni of dimension n i . The task variation ėi is linked to v by a time-varying Jacobian matrix J i such that ėi = J i v.

(

In order to set the values of the tasks' variations to their desired values ėd i , the optimal velocity can be computed as

v = J + ėd , (2)
with

J =    J 1 . . . J N    and ėd =    ėd 1 . . . ėd N    , ( 3
)
where + is a pseudo-inverse operator. In practice we use the damped pseudo-inverse [START_REF] Deo | Overview of damped least-squares methods for inverse kinematics of robot manipulators[END_REF] in order to avoid large velocities that can appear when J becomes near singular. The behavior of the controller then depends on the design of the different tasks. We choose here to design three tasks (N = 3) to control different aspects of the needle insertion.

The first task e 1 is used to control the safety aspect of the insertion. It is defined as the angle θ, represented in Fig. 2, between the needle base axis and the insertion point p in at the surface of the tissue, such that

θ = -atan2 x 2 in + y 2 in , z in , (4)
where x in , y in and z in are the coordinates of the insertion point p in expressed in the needle base frame (see Fig. 2).

Using the geometric differenciation of θ with respect to v, the Jacobian matrix J in ∈ R 1×6 associated to the task e 1 by (1) can be computed according to

J in = - x in z in ld 2 - y in z in ld 2 l d 2 y in l - x in l 0 , ( 5
)
where d = x 2 in + y 2 in + z 2
in is the distance between the needle base and the insertion point and l = x 2 in + y 2 in is the lateral distance between the needle base axis and the insertion point, as illustrated in Fig. 2.

This task ensures that the needle base automatically stays aligned with the insertion point, hence limiting the tissue deformation near the surface. The desired variation ėd 1 of the task is set to regulate θ toward zero according to

ėd 1 = -λ θ θ, (6)
where λ θ is a positive gain that tunes the decrease rate of θ. The second task e 2 ∈ R 3 is used to control the trajectory of the needle tip and is defined as the position of the tip. The Jacobian matrix J tip ∈ R 3×6 associated to the task e 2 is defined such that

ė2 = v tip = J tip v, (7)
where v tip ∈ R 3 is the velocity of the needle tip expressed in the tip frame. Note that J tip depends on the properties of the needle and of the tissues, such that a specific method needs to be used for its computation, as will be detailed in section II-B. This task is used to include the surgeon in the control loop via the haptic interface. In this study we will compare two different ways to compute the desired velocity v d tip ∈ R 3 , to see if it is better to let the user control the 3 DOF of the tip or only the insertion speed.

In the first case, full control over the tip trajectory is given to the user and v d tip is computed according to

ėd 2 = v d tip = α tip R h p h -p h,0 , (8)
where α is a positive scaling factor used to convert a distance into a velocity, tip R h ∈ SO(3) is the rotation matrix from the tip frame to the frame of the haptic device (see Fig. 3), p h is the current position of the handle of the haptic device and p h,0 is a default position at the center of the haptic device workspace, which will be used as a reference position for both the control law and the haptic feedback (see section II-D.2). Using this kind of control law, the more the user moves the handle of the haptic device away from the default position, the higher is the velocity applied to the tip. The user can thus freely control the trajectory of the tip and they are in charge of ensuring that it is effectively driven toward the target.

In the second case, only the insertion velocity is controlled by the user, while the direction of the tip velocity is maintained oriented toward the target. v d tip is then computed as

ėd 2 = v d tip = V t p tar p tar , (9)
with

V t = α tip R h p h -p h,0 .d tip , ( 10
)
where V t is the amplitude of the tip velocity, p tar ∈ R 3 is the position of the target and d tip ∈ R 3 is the direction of the tip axis (i.e. z t in Fig. 2), both expressed in the tip frame. The user can thus choose the insertion velocity V t , while the system automatically adapts the tip trajectory to reach the target. Note that the dot product with d tip in [START_REF] Moreira | The miriam robot: A novel robotic system for mr-guided needle insertion in the prostate[END_REF] ensures that only the component of the motion of the handle along the tip direction is taken into account to compute the norm of the applied tip velocity, while lateral components are discarded in order to keep visual consistency. This way V t can also be negative if the user wants to retract the needle.

Finally the third task e 3 is used to automatically ensure that the circular tip trajectory induced by the bevel during the insertion is oriented in the desired direction, reducing the needle bending required to obtain the tip motion. It is defined as the angle γ, represented in Fig. 2, between the desired tip velocity v d tip and the cutting-edge of the bevel:

e 3 = γ = atan2(v d tip,y , v d tip,x ) - π 2 , ( 11
)
where v d tip,x and v d tip,y are the components of v d tip along the x t and y t axis of the tip frame, respectively (see Fig. 2). Assuming a torsionally stiff needle, the rotation of the needle base around the needle shaft is directly transmitted to the needle tip, such that the Jacobian matrix J b ∈ R 1×6 associated to the task e 3 by (1) can be computed as

J b = 0 0 0 0 0 1 . ( 12
)
The desired variation e d 3 of the task is set such that γ is always regulated toward zero according to

ėd 3 = -λ γ γ, ( 13
)
where λ γ is a positive gain that tunes the decrease rate of γ. The final velocity screw vector applied to the needle base is then computed using [START_REF] Misra | Mechanics of flexible needles robotically steered through soft tissue[END_REF]. The control law enables the user to control the state of the insertion in real-time, while leaving the complexity of the needle base manipulation to the system.

The Jacobian matrices associated to the first and third tasks can be computed based on geometry. However, the Jacobian matrix J tip associated to the second task needs to be estimated using a model of the needle interaction with the tissues. We describe in the following section II-B the mechanics-based model that we use to compute J tip .

B. Needle insertion modeling

The 6 DOF of the needle base are controlled in this framework, so that lateral motions can be applied to the base. Due to the needle flexibility and its interaction with the tissues all along its shaft, the needle shape can be modified inside and outside of the tissues, so that the widely used kinematic models of unicycle [START_REF] Webster | Nonholonomic modeling of needle steering[END_REF] can not be used here.

Instead, we model the whole needle body in order to estimate the effect of the needle base motion on the tip motion. We use the same model as was introduced in [START_REF] Chevrie | Online prediction of needle shape deformation in moving soft tissues from visual feedback[END_REF], consisting of a flexible beam representing the needle that interacts with a second beam representing the path cut in the tissues by the beveled tip (see Fig 2). This second beam is updated online to take into account the non-linear properties and the external motions of the tissues using an unscented Kalman filter. In [START_REF] Chevrie | Online prediction of needle shape deformation in moving soft tissues from visual feedback[END_REF] CCD cameras were used to observe the needle through translucent gelatin and it was shown that good estimation performances with mean accuracy under 1 mm could be obtained even when the acquisition rate was lowered to 1 Hz. Here we use instead 3D US imaging to observe the needle in non translucent soft tissues, as this medical imaging modality is commonly used in medical practice to perform needle insertion procedures. The position of several points along the needle shaft, necessary for the model update, are computed using the needle tracking algorithm in 3D US volumes presented in the following section II-C.

A finite difference method is then used with the model to compute the Jacobian matrix J tip linking the motion of the needle base to the tip velocity as defined in [START_REF] Adagolodjo | Inverse real-time finite element simulation for robotic control of flexible needle insertion in deformable tissues[END_REF].

C. Needle and target tracking in 3D ultrasound

As the needle and target tracking is not the main focus of this paper, we only briefly describe the principle of the algorithms that we use. More details can be found in [START_REF] Chevrie | Flexible Needle Steering using Ultrasound Visual Servoing[END_REF]. The needle is tracked in each US volume using an iterative algorithm to locally fit a polynomial curve to the best location of the needle. The curve is first initialized using the current state of the model presented in section II-B. At each iteration, the curve is deformed to maximize an objective function that is computed from the voxel intensity around the curve and that takes into account the typical needle-induced artifacts in the US images. A spherical target is also tracked in the 3D US volumes using a 3D version of a Star algorithm [START_REF] Friedland | Automatic ventricular cavity boundary detection from sequential ultrasound images using simulated annealing[END_REF].

D. Teleoperation and user interface

In order to close the control loop and to ensure the success of the needle insertion, feedback needs to be provided to the user all along the procedure. Visual and haptic feedback are usually used together to provide complementary information on the state of the insertion. In the following we describe the visual and haptic information that we provide to the user in order to perform the teleoperated needle insertion.

1) Visual feedback: In order to monitor the good proceeding of the insertion, a visual feedback is provided to the user in the form of three orthogonal slices of the last acquired US volume, as can be seen in Fig. 3. The slices are automatically selected such that they intersect each other at the current location of the tip of the needle model, such that the tip is always visible during the insertion. Additional information are also displayed on each view in order to help the surgeon interpreting the images and performing the insertion toward the target. The needle model is projected onto each view to give the currently estimated position of the needle body. The results of the needle and target tracking algorithms are also projected on each view, so that it is easy to identify the direction in which the needle should go to reach the target.

Since the controller defined in section II-A applies a tip velocity that is proportional to the position of the handle of the haptic device, it is important to provide a visual feedback of this position. The point corresponding to the default position p h,0 of the haptic device is mapped on each view such that it corresponds to the current position of the needle tip. The position p h of the handle is then displayed on each image using a similar method as was used in [START_REF] Xia | Looking beyond the imaging plane: 3d needle tracking with a linear array ultrasound probe[END_REF], i.e. using a cross with a size equal to the orthogonal distance of p h from the plane and with a color that depends on which side of the image plane p h is lying. A transformation is first applied to p h such that the motion applied by the user on the haptic device is aligned with the motion of the cross in the XY view of the US volume (see Fig. 3). This leads to an intuitive control of the direction of the tip velocity in the XY view, independently of the real orientation of the US probe, haptic device or needle in the world frame. Note that the cross will move at the same speed as the needle tip if the user does not change the position of the haptic device.

2) Haptic feedback: We compare three different kinds of force feedback to provide various levels of guidance to the user during the insertion. Each feedback consists in a force f h ∈ R 3 that drives the handle of the haptic device toward its default position p h,0 defined previously in [START_REF] Kaye | Robotic ultrasound and needle guidance for prostate cancer management: review of the contemporary literature[END_REF]. The applied force is linked to the current position p h of the haptic device by a positive definite stiffness matrix

K h ∈ R 3×3 such that f h = -K h p h -p h,0 . (14)
Using this method, the intensity of the feedback provided to the user is thus proportional to the velocity that is currently applied to the needle tip by the controller defined in [START_REF] Kaye | Robotic ultrasound and needle guidance for prostate cancer management: review of the contemporary literature[END_REF]. This ensures the safety of the framework by natural driving the haptic device toward a position where no tip velocity is applied by the controller. Depending on the exact definition of K h , which can be non-diagonal, different guiding effects can be created by applying a different stiffness in different directions. We describe in the following the three different methods used to define K h and represented in Fig. 4.

a) Method 1: We define a stiffness matrix K iso h such that it provides an isotropic stiffness, i.e. the force provided to the user is directed toward the haptic center and is proportional to the distance from this center. K iso h is thus defined as a diagonal matrix with a constant factor k in the diagonal. This way the feedback is directly proportional to the velocity applied to the tip. Using this definition, the user has a full control over the trajectory of the needle tip and does not feel constrained in any direction. However, no haptic information is provided on the target position or on the degree of feasibility of the tip motion.

b) Method 2: We define a stiffness matrix K tip h such that the stiffness is anisotropic with a greater stiffness applied in directions in which it is difficult to move the tip without deforming the tissues and with a low stiffness in the other directions. This allows the user to be naturally guided in the direction in which the needle tip can be inserted without bending the needle and pushing laterally on the tissues. This method should increase the safety of the insertion procedure by limiting the lateral motions of the needle.

We use the Jacobian matrix J tip,t linking the base velocity to the tip velocity to identify the manipulability ellipsoid at the tip [START_REF] Yoshikawa | Manipulability of robotic mechanisms[END_REF]. J tip,t is directly extracted from the left half of the Jacobian matrix J tip computed in [START_REF] Adagolodjo | Inverse real-time finite element simulation for robotic control of flexible needle insertion in deformable tissues[END_REF]. The manipulability ellipsoid can be characterized by three unitary orthogonal axis u i ∈ R 3 each associated to one of the singular values σ i of J tip,t . A low value for σ i indicates that moving the tip along axis u i requires a large motion of the base and vice versa. Due to the high axial stiffness of the needle, a direction for which the base velocity is transmitted to the tip always exists, so that one of the σ i should always be equal to 1 and the associated u i indicates the natural insertion direction of the tip. Then we associate a stiffness k i for each axis such that

k i = k + σ i , ( 15
)
where k is the default isotropic stiffness defined previously and < 1 is used to bound the value of the stiffness when σ i becomes too small. This way, the natural insertion direction is associated to a low stiffness near k, while the other directions are associated to a stiffness that can grow up to a maximum value of k . We empirically chose = 0.1 in the experiments, which was enough to provide guidance without inducing stability issues in the haptic force feedback. The stiffness matrix K tip h can then be computed as

K tip h = tip R T h U   k 1 0 0 0 k 2 0 0 0 k 3   U T tip R h , (16)
where U ∈ R 3×3 is a matrix containing the column vectors u i and tip R h was defined in [START_REF] Kaye | Robotic ultrasound and needle guidance for prostate cancer management: review of the contemporary literature[END_REF].

Using such a haptic feedback, the user can feel which direction is likely to cause the least bending of the needle and the least tissue deformation. At the beginning of the insertion lateral motions of the tip can be easily achieved, so that the applied stiffness is almost isotropic and the user is not much constrained. On the contrary, as the insertion depth increases, lateral motions of the tip require more lateral base motions, so that the user feels more guided in a specific direction which is the natural insertion direction of the needle. c) Method 3: We define the anisotropic stiffness matrix K tar h in a similar way to Method 2, except that the low stiffness axis is realigned with the target (see Fig. 4):

K tar h = tar R T h U   k 1 0 0 0 k 2 0 0 0 k 3   U T tar R h , (17)
where tar R h ∈ SO(3) is the rotation matrix from the tip frame realigned with the target to the frame of the haptic device and the other terms were defined in [START_REF] Mignon | Beveled-tip needle-steering using 3d ultrasound, mechanical-based kalman filter and curvilinear roi prediction[END_REF]. Using this feedback the user is naturally guided toward the target during the insertion.

III. EXPERIMENTS In this section we describe the experiments that we performed to validate the teleoperation framework and we present the results of the comparison between the different control and feedback methods.

A. Experimental setup

The setup used for the experiments can be seen in Fig. 1. We used a Chiba biopsy needle (Angiotech MCN2208) with internal diameter 0.48 mm and external diameter 0.7 mm attached to the end effector of a 6 DOF manipulator Target Fig. 4. Representation of the different definitions of the stiffness for the haptic force feedback, overlaid on the needle tip. From left to right: isotropic stiffness (Method 1), anisotropic stiffness in the tip direction (Method 2) and anisotropic stiffness in the target direction (Method 3). The ellipses represent a set of velocity vector that would induce a force feedback of constant intensity. Green and red arrows indicate velocity directions in which the intensity of the provided force feedback is low or high, respectively.

Viper S650 (Adept). The length of the needle that can bend is 12.6 cm long. The needle Young's modulus is assumed to be that of steel, i.e. 200 GPa. The insertion is done in a homemade porcine gelatin phantom embedding several spherical targets. Visual feedback is provided to the controller using a SonixTOUCH Research 3D US scanner (Ultrasonix Medical Corporation) with a 4DC7-3/40 motorized 3D US probe maintained fix during the experiments. The acquisition parameters are set such that a volume has a depth of 8 cm, a frame field of view of 78 • and a wobbling angle of 37 • . A new volume is acquired every 408 ms. The US station allows us to get access to the prescan data acquired by the probe. The 3D volumes are reconstructed in cartesian coordinates [START_REF] Lee | Intensity-based visual servoing for non-rigid motion compensation of soft tissue structures due to physiological motion using 4d ultrasound[END_REF] after a 3 × 3 × 3 median filter is applied to remove some noise in the prescan data . A Falcon haptic interface (Novint Technologies) is used to get the user input and to provide the haptic feedback.

B. Experimental scenario

We compare four different combinations of user tasks and haptic feedback. The first three combinations use the 3 DOF user task (see ( 8)) along with one of the methods to compute the haptic feedback defined in section II-D.2 (with the same indexes). The fourth combination (Method 4) uses the user task limited to the insertion velocity (see [START_REF] Piccin | A force feedback teleoperated needle insertion device for percutaneous procedures[END_REF]) along with the isotropic stiffness Note that the anisotropic stiffness is not considered in this method since the insertion direction is already imposed by the controller. In each case the default stiffness is set to k = 500 N.m -1 and the position to velocity scaling factor is set to α = 0.1s -1 . For each insertion, the needle tip is first placed at the surface of the tissues and the interaction model is initialized. The needle is then inserted 8 mm into the tissues so that the user can initialize the needle and target tracking algorithms in the US volume displayed on the screen. The user can then start and stop the controller by clicking on a button on the haptic device.

C. Results

In this section we present and discuss the results of the experiments. Five insertions were performed for each method. Trajectories obtained during one of the insertions performed with Method 1 are shown in Fig. 5. We can see that the user is able to reach the target at the end of the insertion. Note that the trajectory obtained by integrating the 5. Trajectories in the frame of the ultrasound probe during a typical example of insertion using the 3 DOF tip control and without haptic guidance (isotropic stiffness). The black curve represents the trajectory of the target, the red curve is the trajectory of the tip of the needle model, the green curve is the trajectory of the tracked needle tip and the blue curve is the position obtained by integrating the velocity input given by the user through the haptic interface. input of the user diverges from the actual trajectory of the tip, mainly due to the damping introduced in the control law [START_REF] Misra | Mechanics of flexible needles robotically steered through soft tissue[END_REF]. However this does not seem to hinder the targeting performance since the visual loop is closed by the user. We first compare the performance in terms of targeting accuracy, segmented manually in the US volume at the end of the insertion. Table I gives the mean final targeting accuracy obtained across the 5 insertions for each method. Good targeting performance can be obtained with all methods, with a mean targetting error under 2.5 mm and a mean lateral error under 2.0 mm, which is on par with the accuracy obtained in current research work on needle steering under US feedback and is sufficient for most medical applications.

In order to increase the safety of the insertion procedure, it is important to avoid breaking the needle or tearing the tissues. We thus compared the mean and maximum values of the needle bending and of the tissue deformation obtained during the insertion with the different methods. The needle bending is taken as the integral of the absolute value of the needle curvature along the whole needle shaft. The tissue deformation is taken as the mean displacement of the tissue from its initial position along the inserted needle (see Fig. 2). Both are estimated using the current state of the interaction model presented in section II-B. The average and maximum values obtained for each method are summarized in Table II.

As expected, the amount of needle bending and tissue deformation is greatly reduced when the tip motion is fully constrained by the controller (Method 4) compared to the 3 DOF control. This is mainly due to the fact that the controller aligns the needle with the target as soon as the insertion begins, so that little lateral motion of the needle base is then required for the remaining of the insertion. Guiding the user toward the target (Method 3) also seems to decrease the induced amount of needle bending and tissue deformation compared to the non-guided Method 1. This validates the usefulness of providing a haptic force feedback to the user to guide the insertion. However, contrary to what could be expected, guiding the user along the direction of least lateral needle base motion (Method 2) seems to actually increase the amount of induced needle bending and tissue deformation. This can actually be explained by the same reasons that justified the good behavior of Method 4. Using this kind of guidance at the beginning of the insertion will indeed tend to delay the alignment of the needle with the target, which will then require more tissue deformation to be achieved later when the needle is deeper inserted.

Overall the results validate the good performance of the framework in terms of targeting accuracy and stress out the importance of carefully choosing the way haptic guidance is provided to the surgeon in order to ensure a safe procedure.

IV. CONCLUSIONS AND FUTURE WORK

In this paper we presented a proof-of-concept for a teleoperation framework to perform the real-time robotic insertion of a beveled-tip flexible needle under 3D ultrasound guidance and haptic feedback. The framework enables the user to intuitively guide the trajectory of the needle tip in the image while the controller handles the complexity of the 6D motion that needs to be applied to the needle base. A mean targeting accuracy of 2.5 mm could be achieved in gelatin phantoms, validating the whole framework. We compared different ways to provide the haptic feedback as well as different levels of control given to the user on the tip trajectory. Limiting the user input to the insertion speed while automatically controlling the trajectory of the needle tip seems to provide a safer insertion process, however it may be too constraining and can not handle situations where more control over the tip trajectory is required, for example if unpredicted obstacles need to be avoided. On the contrary, giving the full control of the 3D tip velocity to the user and applying a haptic feedback to guide the user toward the target proved to maintain a low level of needle bending and tissue deformation.

Since the results are user-dependent, future work would include a larger scale user study with experts and beginners in order to confirm the obtained performance. As the different components of the framework can already handle the effect of external tissue motion with little to no adaptation, tests should also be conducted to validate its ability to perform automatic motion compensation during the teleoperation.

Fig. 1 .

 1 Fig.1. Control framework and setup used to perform a teleoperated needle insertion. Elements in blue interact together to provide the visual interface: the target and needle tracking algorithms, the ultrasound probe and the needle/tissue interaction model send data to be displayed on the screen. Elements in red interact together to provide the haptic interface: the target tracking algorithm and the needle/tissue interaction model are used to compute the haptic feedback sent to the haptic device, while the haptic device sends position data to the controller.

Fig. 2 .

 2 Fig.2. Illustration of the different elements used in the task functions for the control: angle θ between the needle base axis and the initial insertion point p in , angle γ between the cutting edge of the bevel and desired tip velocity v d tip . Tissue deformation is computed as the mean displacement of the tissue from its initial position along the inserted part of the needle.

Fig. 3 .

 3 Fig. 3. Picture of the screen view providing the visual feedback from the ultrasound probe. Three orthogonal slices of the current US volume passing by the model needle tip are displayed. The needle model is represented by the green lines, the needle tracking is represented by the red lines and the target center and boundaries are represented by the green dots and green circles. The current position of the haptic device is represented by blue crosses when it is in front of the image plane and by yellow crosses when it is behind the image plane, and the size of the crosses corresponds to the orthogonal distance from the image plane.

  Fig.5. Trajectories in the frame of the ultrasound probe during a typical example of insertion using the 3 DOF tip control and without haptic guidance (isotropic stiffness). The black curve represents the trajectory of the target, the red curve is the trajectory of the tip of the needle model, the green curve is the trajectory of the tracked needle tip and the blue curve is the position obtained by integrating the velocity input given by the user through the haptic interface.

TABLE I

 I

		COMPARISON OF THE FINAL TARGETING ERROR
	Method Axial error (mm) Lateral error (mm) Total error (mm)
	1	1.3±0.7	1.7±0.5	2.5±0.6
	2	-0.4±1.5	1.5±0.5	2.0±0.6
	3	0.1±1.8	1.6±0.4	2.3±0.5
	4	-1.3±1.1	2.0±0.6	2.5±1.0


Grégoire David

Frank Wennmohs

Frank Neese

Nicolas Ferré
email: nicolas.ferre@univ-amu.fr

Chemical Tuning of Magnetic Exchange Couplings Using Broken-Symmetry Density Functional Theory

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

Introduction

Magnetic exchange coupling J is one of the most important quantum interactions between two (or more) unpaired electrons. [START_REF] Kahn | Molecular magnetism[END_REF] In the case of a finite-size systems, this interaction is important for designing molecular compounds with controlled magnetic properties. When J is non-zero in a molecule, it favors either the parallel alignment of the magnetic electron spin angular momenta (ferromagnetic interaction) or the anti-parallel one (antiferromagnetic interaction), eventually driving the magnetic properties of the system. Accordingly, the determination of J, either experimentally [START_REF] Gatteschi | Molecular Nanomagnets[END_REF][START_REF] Pedersen | Single-molecule magnet engineering: buildingblock approaches[END_REF][START_REF] Demir | Radical ligand-containing single-molecule magnets[END_REF] or by means of quantum-chemistry calculations, [START_REF] Hay | Orbital Interactions in Metal Dimer Complexes[END_REF][START_REF] Noodleman | Valence bond description of antiferromagnetic coupling in transition metal dimers[END_REF][START_REF] Caballol | Remarks on the Proper Use of the Broken Symmetry Approach to Magnetic Coupling[END_REF][START_REF] Barone | Through-bond and through-space effects in the magnetic properties of nitroxide biradicals by an integrated QM/MM approach including solvent effects[END_REF][START_REF] Ruiz | Further Theoretical Evidence for the Exceptionally Strong Ferromagnetic Coupling in Oxo-Bridged Cu(II) Dinuclear Complexes[END_REF][START_REF] Ciofini | DFT calculations of molecular magnetic properties of coordination compounds[END_REF][START_REF] Illas | Extent and limitations of densityfunctional theory in describing magnetic systems[END_REF][START_REF] Moreira | A unified view of the theoretical description of magnetic coupling in molecular chemistry and solid state physics[END_REF][START_REF] Queralt | On the applicability of multireference second-order perturbation theory to study weak magnetic coupling in molecular complexes[END_REF][START_REF] Vancoillie | Theoretical Description of the Structure and Magnetic Properties of Nitroxide-Cu(II)-Nitroxide Spin Triads by Means of Multiconfigurational Ab Initio Calculations[END_REF][START_REF] Neese | Prediction of molecular properties and molecular spectroscopy with density functional theory: From fundamental theory to exchange-coupling[END_REF][START_REF] Phillips | Magnetic Exchange Couplings from Semilocal Functionals Evaluated Nonself-Consistently on Hybrid Densities: Insights on Relative Importance of Exchange, Correlation, and Delocalization[END_REF][START_REF] Malrieu | Magnetic Interactions in Molecules and Highly Correlated Materials: Physical Content, Analytical Derivation, and Rigorous Extraction of Magnetic Hamiltonians[END_REF][START_REF] De Graaf | Magnetic Interactions in Molecules and Solids (Theoretical Chemistry and Computational Modelling)[END_REF][START_REF] Malrieu | Communication: Proper use of broken-symmetry calculations in antiferromagnetic polyradicals[END_REF][START_REF] Roemelt | Exchange Coupling Interactions from the Density Matrix Renormalization Group and N-Electron Valence Perturbation Theory: Application to a Biomimetic Mixed-Valence Manganese Complex[END_REF] is of tremendous importance for understanding the behavior of magnetic open-shell molecular systems.

Different physical mechanisms may compete and determine the sign and intensity of J. In the case of two spin-1/2 magnetic centers, the most important ones are the direct exchange and the kinetic exchange between the magnetic electrons, and the polarization of the nonmagnetic electrons. While high level quantum-mechanical methods can accurately evaluate the importance of these mechanisms, eg by selectively considering classes of excitations in configuration interaction calculations, [START_REF] Terencio | Physical analysis of the through-ligand long-distance magnetic coupling: spin-polarization versus Anderson mechanism[END_REF] their computational cost limits their applicability to small molecules. On the other hand, the Kohn-Sham density-functional theory combined with the so-called broken-symmetry approach (BS-KSDFT) is a popular and effective method which often provides semi-quantitative evaluations of J. [START_REF] Caballol | Remarks on the Proper Use of the Broken Symmetry Approach to Magnetic Coupling[END_REF][START_REF] Neese | Prediction of molecular properties and molecular spectroscopy with density functional theory: From fundamental theory to exchange-coupling[END_REF][START_REF] Sinnecker | Calculating the Electron Paramagnetic Resonance Parameters of Exchange Coupled Transition Metal Complexes Using Broken Symmetry Density Functional Theory: Application to a MnIII/MnIV Model Compound[END_REF] By means of selective freezing of a subset of orbitals, some of the present authors recently introduced a decomposition scheme [START_REF] Coulaud | Analysis of the physical contributions to magnetic couplings in broken symmetry density functional theory approach[END_REF][START_REF] Coulaud | Additive Decomposition of the Physical Components of the Magnetic Coupling from Broken Symmetry Density Functional Theory Calculations[END_REF][START_REF] Ferré | Spin decontamination of broken-symmetry density functional theory calculations: deeper insight and new formulations[END_REF] which has been essentially applied to centro-symmetric molecules until now. [START_REF] Guennic | Analysis of the Magnetic Exchange Interaction in Halide-Bridged Cu(II) Binuclear Complexes: Deciphering the Paths[END_REF] Thanks to the recent implementation of this method in the Orca package, [START_REF] Neese | The ORCA program system[END_REF] the J decomposition is now available for molecules belonging to various symmetry groups, hence expanding its range of application to systems where the magnetic centers are bridged by chemically different moieties, eg π-conjugated or not, with electron-donating or electron-withdrawing substituents, etc. Hereafter, we report J decompositions for various non-centrosymmetric copper dinuclear compounds. All of them feature oxygen-or carboxylic-bridged cores, with different chemical substitutions. In all, it is demonstrated that the deep understanding of the physical contributions to J can be used to in-silico design molecular compounds with selected magnetic properties.

Molecular models

Among all the possible, available, non-centrosymmetric copper dinuclear complexes that exhibit significant magnetic exchange couplings, we have selected [{Cu(dmen

)} 2 (μ-OMe){μ-[O 2 C(η 5 -C 5 H 5 )
(hereafter identified as PATFIA) and its closely related analogue [{Cu(tmen

)} 2 (μ-OH)(μ-O 2 ClO 2 ){μ-[O 2 C
(hereafter identified as PATFOG). Both PATFIA and PATFOG, whose synthesis have been reported in 2005, [START_REF] López | Magneto-structural correlations in binuclear copper(ii) compounds bridged by a ferrocenecarboxylato(-1) and an hydroxo-or methoxo-ligands[END_REF] feature two different bridges (Figure 1). Moreover, inspection of the corresponding experimental crystallographic structures reveal a close contact with 1 or 2 perchlorate counterions. From experiment, the two copper ions in PATFIA have been shown to be weakly antiferromagnetically coupled (J = -5.5 cm -1 ), while PATFOG shows ferromagnetic coupling (J = 14.7 cm -1 ), both J values being defined through the H = -2JS A •S B spin Hamiltonian. In the same article, 28 quantum chemical calculations have been used to understand the origin of the different couplings in PATFIA and PATFOG. The out-of-plane angle of the hydroxo-based bridge side-chain is found to be the principal structural param-eter controlling the ferro-vs antiferromagnetic character of PATFOG vs PATFIA. Since the publication of this article, PATFIA has been repeatedly used as a benchmark molecule for evaluating the accuracy of various quantum-mechanical methods. The corresponding J values are reported in Table 1, together with the model details. In almost every reported studies, the perchlorate counterions have been ignored, while the ferrocene moiety is simply replaced by a methyl group or a hydrogen atom. Looking at Table 1 entries, it is quite obvious that the correct evaluation of PATFIA J value depends on a large number of parameters: the DFT xc-functional, the geometry, the basis set, the method (BS vs other ones) etc. While some models perform wrong, even qualitatively (ferromagnetic instead of anti-ferromagnetic), some other ones overshoot the coupling by two orders of magnitude. On the other hand, besides the original calculations performed on PATFOG, [START_REF] López | Magneto-structural correlations in binuclear copper(ii) compounds bridged by a ferrocenecarboxylato(-1) and an hydroxo-or methoxo-ligands[END_REF] we are aware of only a very recent theoretical study dealing with this system, 36 for which a large set of DFT functionals have been tested, with J values ranging from -3.6 cm -1 to 59.9 cm -1 . Notice that all but one DFT functional predict PATFOG as ferromagnetic.

In the present work, the starting structures have been retrieved from the Cambridge Structural Database. 37 Comparing PATFIA and PATFOG structures, several differences show up. First, one of the bridges is hydroxo in PATFOG, while it is methoxo in PATFIA.

In the case of PATFOG, the hydroxo moiety is oriented in such a way it can form a hydrogen bond with one perchlorate molecule. This orientation is thought to govern the magnetic coupling character of these systems. Second, external ligands are dimethylethylenediamine in PATFIA, while they are tetramethylethylenediamine in PATFOG, hence feature two more methyl groups. Accordingly, it would be interesting to assess the influence on J of these two extra chemical groups. Third, the perchlorate anions are distributed differently in PATFIA and PATFOG, the latter featuring a hydrogen bond with the hydroxo bridge. While these anions are usually removed from the theoretical models, the original work from López et al [START_REF] López | Magneto-structural correlations in binuclear copper(ii) compounds bridged by a ferrocenecarboxylato(-1) and an hydroxo-or methoxo-ligands[END_REF] insisted on the structural importance of perchlorate counterions, constraining the orientation of the hydroxo bridge, ultimately modifying PATFOG magnetic properties. Together with the evaluation of the ferrocene effect on J (which has been recently confirmed to be innocent [START_REF] Orms | Singlet-triplet energy gaps and the degree of diradical character in binuclear copper molecular magnets characterized by spin-flip density functional theory[END_REF] ), we intend to study how these chemical differences and substitutions influence the different physical contributions to the magnetic exchange couplings in PATFIA and PATFOG, thanks to the BS-KSDFT-based decomposition scheme, as presented in the next section.

Methods

Our J decomposition scheme is based on successive BS-KSDFT calculations, optimizing either restricted open-shell (RO) or unrestricted (U) Kohn-Sham wavefunctions under the constraint that a subset of orbitals (either the magnetic or the non-magnetic ones) are kept frozen. In the case of two magnetic centers A and B, each of them carrying a single unpaired electron, we make use of the following effective spin-Hamiltonian:

H = -2JS A • S B (1)
for which J can be decomposed using the Heisenberg-Dirac-van Vleck model as:

J = J 0 + ∆J ke + ∆J cp (2)
ie, J is the sum of three important contributions: the direct exchange, the kinetic exchange and the so-called core polarization of the non-magnetic orbitals. The first step of the procedure consists in the calculation of the magnetic orbitals. Using a high-spin state (a triplet, in the case of two spin-1/2 sites) RO solutions (T,RO), both singly-occupied molecular orbitals are localized, using one of the available localizations schemes in Orca. [START_REF] Neese | The ORCA program system[END_REF] Such schemes include the popular Pipek-Mezey 38 and Boys 39 localization schemes as well as the more recently developed method of intrinsic atomic orbitals. All schemes are combined with a robust and efficient search algorithm that is designed to avoid local minima of the localization functional to the largest extent.

The remaining steps are hereafter summarized:

1. Flip the spin of one magnetic orbital to produce the broken-symmetry M s = 0 wavefunction (BS,RO), without optimizing any orbital, and evaluate the direct exchange contribution J 0 :

J 0 = E BS,RO -E T,RO (3)
2. Relax the magnetic orbitals in the M s = 0 wavefunction, in the field of the frozen non-magnetic orbitals, allowing the magnetic orbitals to delocalize from one center to the other (BS,UFC):

∆J ke = 1 S 2 T,RO -S 2 BS,U F C (E BS,U F C -E T,RO ) -J 0 (4)
3. Relax the core orbitals in both high-spin (M s = 1) and broken-symmetry M s = 0 wavefunctions (respectively T,UFM and BS,UFM):

∆J cp = (E BS,U F M -E T,U F M ) 2 -S 2 BS,U F C + S 2 BS,U F M /2 + S 2 BS,U F C S 2 BS,U F M -S 2 BS,U F C /2 -J 0 -∆J ke (5)
4. Total magnetic coupling is obtained with the complete relaxation of orbitals in both high-spin (M s = 1) and broken-symmetry M s = 0 wavefunctions (respectively T,U and BS,U):

J T ot = (E BS,U -E T,U ) 2 -S 2 BS,U F C + S 2 BS,U /2 + S 2 BS,U F C S 2 BS,U -S 2 BS,U F C /2 (6)
Freezing orbitals is achieved thanks to the Local Self-Consistent-Field method, [START_REF] Assfeld | Quantum chemical computations on parts of large molecules: the ab initio local self consistent field method[END_REF] in which the relaxed orbitals are variationally optimized under the constraint of being orthogonal to the frozen ones.

Because the various orbital relaxations are performed in the unrestricted formalism, the resulting Kohn-Sham wavefunctions are spin-contaminated, ie characterized by a S 2 value different from 2 (high-spin) or 0 (low-spin). Consequently, we make use of spindecontamination schemes in all the three steps. While the Yamaguchi procedure can be applied safely to the direct and kinetic exchange contributions, it is no longer the case for core polarization, as demonstrated elsewhere. [START_REF] Ferré | Spin decontamination of broken-symmetry density functional theory calculations: deeper insight and new formulations[END_REF] The selection of a xc functional is known to be very delicate in BS-KSDFT J calculations. Even the physical interpretation of the contributions to J depends on the chosen xc functional. [START_REF] David | What Are the Physical Contents of Hubbard and Heisenberg Hamiltonian Interactions Extracted from Broken Symmetry DFT Calculations in Magnetic Compounds?[END_REF] In this work, we have followed the very recent advice by Illas and coworkers, [START_REF] Costa | Post-B3LYP Functionals Do Not Improve the Description of Magnetic Coupling in Cu(II) Dinuclear Complexes[END_REF] recommending to use the very common B3LYP functional for obtaining qualitative (ferromagnetic vs anti-ferromagnetic) to semi-quantitative (a given order of magnitude) results.

Our choice is further motivated by the main goal of the current study, in which we will compare chemically-related compounds which differ only by one or two substitutions. Regarding the basis set, we have also followed recommendations by others, selecting a fairly large basis set for copper atoms (def2-QZVPP) while keeping reasonable ones for the other atoms (def2-SVP). The quality of the basis set and the possible influence of relativistic effects have been considered and are reported in Supporting Information.

Our results are discussed on the basis of the Mulliken atomic spin populations, complemented with spin density maps calculated at different stages of the decomposition scheme.

Results

We have first tested the influence of the ferrocene and solvent molecules on the J decomposition results. Similarly to what has been reported by other (references in Table 1), our results available as Supporting Information indicate we can safely ignore them while retaining a good description of PATFIA and PATFOG magnetic properties.

Effect of chemical substitutions

PATFIA and PATFOG feature two different bridges which chemical nature may modify the coupling between their two copper-centered magnetic electrons. One of them is either hydroxo or methoxo, with only one oxygen atom. Oxygen-bounded moieties will be called R1 hereafter (Figure 2). As already pointed in reference 28, their different orientations are In the following, all results have been obtained with tight convergence thresholds: 10 -8 a.u.for energy and 10 -5 a.u.for one electron energy change. In Table 2, only hydrogen and methyl substitutions have been considered. Test calculations have also been performed with (i) hydrogen in place of the methyl group in R1, re-optimizing the O-H bond length, (ii) methyl in place of hydrogen in R1, re-optimizing the O-C bond length (and the methyl internal coordinates). Original PATFIA corresponds to R 1 --CH 3 and R 2 --H. As previously found, 28 R1 is the main responsible for the magnetic character of these PATFIA-like models. Indeed, when R1 is hydrogen, the J value is consistent with PATFOG one, even if the PATFIA geometry is kept. The direct exchange J 0 and the kinetic exchange ∆J ke are also governed by R1 only, being reduced when hydrogen replaces methyl. However, the kinetic exchange contribution is more intense than the direct exchange when R1 is a methyl group, while this is the contrary when R1 is a hydrogen atom. Polarization of the non-magnetic electrons (∆J cp ) remains small. Regarding the nature of R2, hydrogen or methyl, it looks like all the contributions to J are rather insensitive to it. The last two columns in Table 2 show two interesting features. The magnitude of J 0 depends mainly on the chemical nature of R1, not on its position. Conversely, the R1 position/orientation affects ∆J ke .

In order to maximize the R1 and R2 substitution effects, we now report J decomposition results when either R1 or R2 is a methyl group, while the other bridge is complemented with electro-donating or electro-withdrawing chemical groups. We first consider R2 = CH 3 and some carbon-centered substituents for R1 (Table 3). ∆J cp is predicted to be negligible. Three different kinds of R1 effects can be highlighted.

• Methyl and tertiobutyl: similar direct and kinetic exchange contributions, the latter larger than the former, resulting in an overall anti-ferromagnetic coupling.

• Trimethylformate: reduced direct and kinetic exchange contributions, the latter being still larger than the former, the total magnetic exchange remaining anti-ferromagnetic.

• Trifluoromethyl: even more reduced direct and kinetic exchange contributions, the latter smaller than the former, resulting in an overall ferromagnetic coupling.

Because these R1 groups are all positioned in equivalent orientations, roughly parallel to the CuOCu plane, we can't rely on the usual explanation that was previously used to rationalize the different magnetic couplings in PATFIA and PATFOG. [START_REF] López | Magneto-structural correlations in binuclear copper(ii) compounds bridged by a ferrocenecarboxylato(-1) and an hydroxo-or methoxo-ligands[END_REF] Indeed, the chemical nature of R1 is responsible for the variations of all contributions to J.

Now we consider R1 = CH 3 , as it is the case in PATFIA, and a more extended set of substituents as R2 (Table 4). Whatever the nature of R2, ∆J cp remains negligible with respect to the other contributions, while the direct exchange contribution is almost constant (124 ≤ J 0 ≤ 132 cm -1 ) for the carbon-centered R2 (or H). In the case of the phosphine, this contribution slightly increases to 144 cm -1 . However, the chemical nature of R2 dramatically impacts on the kinetic exchange contribution. Always antiferromagnetic, it decreases (in absolute value) from 173 cm -1 (R2 = CF 3 ) to 50 cm -1 (R2 = P(CH 3 ) 2 ). In two cases (R2 = P(NH 2 ) 3 and P(CH 3 ) 2 ), the PATFIA-like molecule turns out to be ferromagnetic instead of anti-ferromagnetic! This striking result, and, more generally, the chemical tuning of J contributions, deserve further analysis in the next section.

Analysis and discussion

As evidenced in Table 2, the magnetic character of PATFIA-like systems is mainly determined by the nature and the position of R1: anti-ferromagnetic when the bridge is methoxo, while ferromagnetic when it is hydroxo. In the original work by Lopez, [START_REF] López | Magneto-structural correlations in binuclear copper(ii) compounds bridged by a ferrocenecarboxylato(-1) and an hydroxo-or methoxo-ligands[END_REF] but already evidenced earlier, [START_REF] Ruiz | Toward the Prediction of Magnetic Coupling in Molecular Systems: Hydroxo-and Alkoxo-Bridged Cu(II) Binuclear Complexes[END_REF] this observation was explained as resulting from the different orientation of R1, more precisely to its tilt angle with respect to the CuOCu plane. If R1 is close to this plane, the compound is anti-ferromagnetic. If the out-of-plane angle is larger than 35 • , it is ferromagnetic. In our calculations, the same change of magnetic character is due to a larger decrease of the kinetic exchange contribution (about 50%) with respect to the decrease of the direct exchange contribution (about 25%), when going from methoxo to hydroxo. On the other hand, substituting methyl (resp. hydrogen) with hydrogen (resp. methyl) in R1, keeping the same original orientation, indicates that the direct exchange does only depend on the nature of the R1 group. However, the kinetic exchange does depend on both the orientation and nature of R1.

Both contributions can be rationalized by tiny but sizable variations of atomic spin populations (Table 5).

Direct exchange (J 0 ) ferromagnetic contribution is an intrinsic property of the magnetic Table 5: Mulliken atomic spin populations on selected atoms, at various stages of the J decomposition procedure. If R1 is methyl, then the total CH 3 group spin population is reported. Values in italics correspond to non-optimized positions of R1 (out-of-plane CH 3 , in-plane H). Atom labels are given in Figure 2.

Cu 1 Cu 2 O 1 O 21 O 22 R1 R1 = CH 3 T, RO
0.634 0.638 0.132 0.064 0.073 0.012 0.621 0.626 0.168 0.063 0.072 0.028 BS, RO 0.624 -0.627 0.005 0.062 -0.070 -0.000 0.611 -0.616 0.006 0.060 -0.069 0.002 BS, U F C 0.620 -0.624 0.005 0.062 -0.070 -0.000 0.609 -0.614 0.006 0.060 -0.068 0.002 R1 = H T, RO 0.638 0.640 0.138 0.065 0.075 0.006 0.644 0.647 0.118 0.064 0.073 0.010 BS, RO 0.619 -0.621 0.004 0.063 -0.071 0.000 0.628 -0.631 0.004 0.062 -0.070 0.000 BS, U F C 0.617 -0.620 0.004 0.063 -0.071 0.000 0.625 -0.628 0.004 0.062 -0.070 0.000 orbitals as obtained in the T,RO calculation. The more localized they are, the less they interact. Hence the corresponding copper atomic spin populations in the T,RO solution can be used to rationalize the variations of J 0 . When the methyl position is modified from in-plane to out-of-plane, there is a slight redistribution of spin population from the copper atoms to the methoxo ones. Conversely, when the hydrogen position is modified from out-ofplane to in-plane, the spin population redistribution is directed from hydroxo to the copper atoms. When methyl is replaced by hydrogen, keeping the same in-plane orientation, the spin population in copper atoms flows from R1. Conversely, when hydrogen is replaced by methyl, keeping the same out-of-plane orientation, the spin population in copper atoms flows towards R1.

Regarding the kinetic exchange contribution, ie the delocalization of a magnetic orbital towards the other magnetic site, copper spin populations are the most affected. This change is twice larger in the case of R1 = CH 3 than in R1 = H. This effect translates to a larger ∆J ke contribution when R1 is methoxo.

Overall, the methoxo vs hydroxo effect remains small from the electronic point of view.

For this reason, we now turn to the analysis of spin populations in PATFIA-like systems, in which several substitutions in R1 or R2 positions have been considered. Table 3 has shown that both the direct and kinetic exchange contributions are strongly dependent of the nature of the R1 substituent. If methyl and trifluoromethyl are compared, J 0 is divided by two while ∆J ke is reduced by a factor 4. These variations can be understood by inspection of Figure 3, in which iso-contours of the spin density are reported. Regarding the direct exchange contribution, Figure 3a indicates that spin density is more delocalized when R1 = CF 3 , the fluorine centers carrying a noticeable population. In other words, the F electronegative centers attract electron density outside from the bridge Cu -O -Cu, which in turns lowers J 0 . About the kinetic exchange contribution, the comparison of the two pictures in Figure 3b clearly shows that the spin density variation is reduced when R1 = CF 3 , in agreement with a smaller ∆J ke value. In both compounds, this variation is restricted to the copper centers and the atoms directly connected to them. Accordingly, the spin density attracted by the fluorine atoms is not modified by the kinetic exchange contribution.

Table 4 has shown that the kinetic exchange is the only contribution which depends on the R2 substituent. Its intensity is divided by more than 3 when going from R2 = CF 3 to P(CH 3 ) 2 . In Figure ( 4), we have reported the variations of the spin density between the BS,RO and the BS,UFC solutions. In all cases, the spin density redistribution includes Cu, These differences are perfectly in line with the reported ∆J ke changes.

In closing this section we note that the particular effect of trifluoromethyl on the kinetic exchange should be highlighted. When placed at the R1 position, it causes the smallest ∆J ke among all the considered substituents. Conversely, when placed at the R2 position, it causes the largest ∆J ke contribution. In other words, CF 3 electron-rich moiety in R1 drastically reduces the magnetic orbital delocalization, while it enhances it when CF 3 is in R2. In principle, this huge effect could be understood by comparing the Hubbard t (hopping energy) and U (on-site minus different-sites energy difference) values. These parameters are strictly defined in the case of centro-symmetric compounds. [START_REF] David | What Are the Physical Contents of Hubbard and Heisenberg Hamiltonian Interactions Extracted from Broken Symmetry DFT Calculations in Magnetic Compounds?[END_REF] However, the PATFIA-like models being close to symmetric, we have calculated approximate t and U values and checked their relevance by calculating a new value of the kinetic exchange contribution, ∆J CI ke :

∆J CI ke = U - √ U 2 + 16t 2 4 ( 7
)
Perfect agreement between both kinetic exchange contributions fully justifies the relevance occupies both the R1 and R2 positions, it appears that t is mostly affected, and reduced, when CF 3 is in R1 position. In other words, the presence of this electron-rich substituent at position R1 decreases a lot the probability to transfer one magnetic electron from one Cu center to the other one, while it slightly increases it when placed at position R2.

Conclusions

In this article, we have presented how the new implementation in Orca of the magnetic exchange coupling decomposition, based on broken-symmetry Kohn-Sham DFT calculations, opens the way to the computer-assisted design of substituted dinuclear non-centrosymetric compounds with either a ferro-or an anti-ferromagnetic character.

In details, we have confirmed that reduced models of PATFIA and PATFOG related systems are characterized by magnetic contributions qualitatively similar to the ones obtained with the experimental crystallographic structure. In agreement with the observation already made by López and coworkers, [START_REF] López | Magneto-structural correlations in binuclear copper(ii) compounds bridged by a ferrocenecarboxylato(-1) and an hydroxo-or methoxo-ligands[END_REF] the orientation of the hydroxo vs methoxo substituents explains the different magnetic characters of PATFIA and PATFOG. However, the decomposition has allowed us to refine this interpretation by pointing out the importance of the chemical nature of these substituents, especially in the case of the direct exchange contribution.

By means of selective chemical substitutions on each bridge, we have designed PATFIAlike models which, if one can experimentally synthesize them, would exhibit either a ferromagnetic or an anti-ferromagnetic character by tuning one or several contributions to the magnetic exchange coupling. Last, we have rationalized the effect of one particular substituent, trifluoromethyle, which has been shown to drastically reduce the kinetic exchange contribution in one particular position.

Figure 1 :

 1 Figure 1: Molecular structures of PATFIA (left) and PATFOG (right).

Figure 2 :

 2 Figure 2: Molecular structures of modified PATFIA models. R1 and R2 are the sites of substitution. H atoms are omitted for clarity.

  (a) T,RO (isovalue = 0.0015 a.u.) (b) BS,UFC-BS,RO (isovalue = 10 -5 a.u.)

Figure 3 :

 3 Figure 3: Spin densities (a) and spin density variations (b) for two PATFIA-like compounds, with either C(CH 3 ) 3 (left) or CF 3 (right) as R1 substituent.

Figure 4 :

 4 Figure 4: Spin density (isovalue = 10 -5 a.u.) variations (BS,UFC-BS,RO) for three PATFIAlike compounds, with CF 3 , C(CH 3 ) 3 or P(CH 3 ) 2 as R2 substituent.

  the t and U parameters. Whatever the CF 3 position, the U parameters keeps the same value, implying that the energy of the ionic configurations in which both magnetic electrons are on the same Cu center are not affected by the position of the substituent. Conversely, the hopping energy t increases by a factor of two between R1 and R2, explaining the corresponding increase of ∆J ke . If compared to the unsubstituted PATFIA model in which CH 3

Table 1 :

 1 Experimental (in bold) J value complemented with some examples of PATFIA models reported in the literature, with the corresponding level of theory, the type of geometry (either experimental (exp) or optimized (opt), the inclusion of ferrocene (if not, the nature of the substituted fragment) and the corresponding J values in cm -1 . BS: Broken-Symmetry, TD: Time-Dependent, REKS: Restricted-Ensemble Kohn-Sham, SF c : collinear Spin-Flip, SF nc :non-collinear Spin-Flip, CV: Constricted Variational. DFT xc-functional is given in parenthesis, while other model details and basis sets are available in the above references.

	Model	Level of theory	Geometry Ferrocene	J
	Ref 28	H = -2JS A • S B	exp	yes	-5.5
	Ref 28	CASPT2	opt	no, H	-9
	Ref 28	BS(B3LYP)	exp	no, H	-31
	Ref 29	BS(BH&HLYP)	exp	no, H	+5
	Ref 29	REKS(B3LYP)	exp	no, H	+70
	Ref 29 REKS(BH&HLYP)	exp	no, H	+16
	Ref 30	BS(M06)	exp	no, H	-8
	Ref 30	BS(M06-2X)	exp	no, H	-9
	Ref 30	BS(M06-HF)	exp	no, H	-19
	Ref 30	BS(PBE0)	exp	no, H	-18
	Ref 31	BS(HSE)	exp	no, H	-19
	Ref 31	BS(LC-ωPBE)	exp	no, H	-19
	Ref 32	SF c -TD(B3LYP)	exp	no, H	-200
	Ref 32	SF c -TD(B3LYP)	exp	yes	-210
	Ref 32 SF c -TD(B3LYP40)	exp	no, H	-200
	Ref 32 SF c -TD(B3LYP40)	exp	yes	-49
	Ref 32	SF c -TD(BMK)	exp	no, H	-36
	Ref 32	SF c -TD(BMK)	exp	yes	-43
	Ref 32	SF c -TD(M06-2X)	exp	no, H	-11
	Ref 32	SF c -TD(M06-2X)	exp	yes	-12
	Ref 33	SF c -CV(B3LYP)	exp	no, H	-43
	Ref 33 SF c -CV(BH&HLYP)	exp	no, H	-10
	Ref 34	BS(LDA+U)	exp	no, H	+18
	Ref 34	BS(GGA+U)	exp	no, H	-36
	Ref 16	BS(STPSS)	exp	no, H	-307
	Ref 16	BS(PBETPSS)	exp	no, H	-246
	Ref 16	BS(TPSSTPSS)	exp	no, H	-188
	Ref 35 SF c -TD(B5050LYP)	exp	no, H	-56
	Ref 35 SF nc -TD(B5050LYP)	exp	no, H	+2
	Ref 35 SF nc -TD(B5050LYP)	opt	no, H	-25
	Ref 35 SF c -TD(B5050LYP)	exp	yes	-60

Table 2 :

 2 J decomposition (in cm -1 ) for PATFIA-like models, with R1 and R2 as hydrogen or methyl. -H * means that hydrogen is positioned along the original methoxo O-C bond (in-plane), while -CH * 3 means that methyl is positioned along the original hydroxo O-H bond (out-of-plane).

	R1 -CH 3 -CH 3 -H	-H	-H * -CH * 3
	R2 -CH 3	-H	-H -CH 3 -CH 3 -CH 3
	J 0	123	125	93	91	95	127
	∆J ke -154	-146	-71	-75	-121	-111
	∆J cp	-1	0	3	3	-4	9
	J 0 + ∆J ke + ∆J cp	-32	-21	25	19	-30	25

Table 3 :

 3 J decomposition (in cm -1 ) for PATFIA-like models, with R2 = CH 3 and various R1 groups.

	R1 -CF 3 -C(OCH 3 ) 3 -C(CH 3 ) 3 -CH 3
	J 0	58	91	135	123
	∆J ke	-37	-99	-148	-154
	∆J cp	-1	0	4	-1
	J 0 + ∆J ke + ∆J cp	20	-8	-9	-32

Table 4 :

 4 J decomposition (in cm -1 ) for PATFIA-like models, with R1= CH 3 and various R2 groups.

	R2 -CF 3 -CH 3 -H -C(CH 3 ) 3 -C(NH 2 ) 3 -P(CH 3 ) 2
	J 0	125	123	125	124	132	144
	∆J ke -173	-154 -146	-126	-95	-50
	∆J cp	-1	-1	-1	-2	-7	-6
	J 0 + ∆J ke + ∆J cp	-49	-32	-22	-4	30	88

Table 6 :

 6 Kinetic exchange contributions and Hubbard parameters (in cm -1 ) for CF 3 in R1 or R2 position. For comparison, the last line reports the same values for CH 3 in both positions.

		∆J ke -∆J CI ke	U	t
	R1	5.10 -3	31074 -756
	R2	3.10 -3	31077 -1651
	CH

Acknowledgement

GD and NF are grateful to Nathalie Guihéry and Jean-Paul Malrieu (Université Paul Sabatier, Toulouse 3, France) for all the interesting discussions and the method development over the last 10 years. GD thanks Deutscher Akademischer Austauschdienst for generous Marseille-Mulheim travel grant in 2017. Mésocentre of Aix-Marseille Université is acknowledged for allocated HPC resources.

Supporting Information Available

Influence of ferrocene and solvent molecules on J decomposition. Effect of the basis set quality and influence of relativistic effects. Spin densities of some of the compounds. Coordinates of all the compounds.


Crosslinguistic differences in spatial expression have been shown to influence speakers' focus on particular event components in discourse (Slobin 2004; Talmy 2000). For example, Satellite-framed languages express Manner of motion in verb roots and Path in satellites (English 1), whereas Verb-framed languages lexicalize Path in the verb leaving Manner implicit or peripheral (e.g., French 2):

1. He walked into the room

Il est entré dans la pièce [en marchant]).

A debated question is whether such typological differences also influence non-verbal cognition (Gennari et al., 2002; Papafragou & Selimis 2010). The present study addresses this question by comparing verbal and non-verbal responses produced by adult native speakers in two language groups, English and French, differing with respect to motion expression.

RESULTS

The production task (Fig. 2) shows crosslinguistic differences in the structures used by speakers in the two language groups.

-In English: Manner verbs with Path adjuncts (ex. 1).

-In French: Path verbs, less frequent Manner (ex. 2 & 3).

Production

Fig. 2 Information locus in production To measure the relative role of languageindependent and language-specific factors, we compared the verbal and non-verbal performance of adult English and French native speakers.

Subjects (16 per language) performed three tasks involving motion events (Fig. 1).

(

1)Categorization, non-verbal condition:

Participants saw a target cartoon (e.g. a cat walking up a hill), then two variants that differed from it with respect to Manner or Path (walking down vs. jumping up). They then had to choose which variant best matched the target, while simultaneously performing a syllable repetition task that prevented them from internally verbalizing the stimuli.

(2)Categorization, verbal condition: The target was a sentence presented orally (There's a cat walking up a hill), rather than a video (no interference task).

( Interference task « babibo » (continuous syllable repetition)

1.

The horse trotted down the hill.

(a) Un ours qui traverse les rails.'

A bear that crosses the tracks.' (b) Un ours rentre dans la maison en sautillant.

'A bear enters the house whilst hopping'

3.

Un lion court à côté de la voie ferrée.

'A lion runs next to the rail track.'

REFERENCES CONCLUSIONS AND PERSPECTIVES

• Manner types (Fig. 4a)

Both groups rely more on Manner when the Manner in the stimuli is salient:

salient > not salient (jump > walk,run)

• Path types (Fig. 4b)

Both groups rely more on Manner when the Path in the stimuli is INTO/OUT than with other paths.

• Interaction Path x Condition (Fig. 5)

Boundary crossings (INTO, OUT OF, ACROSS) elicit more Manner choices than other Paths for both groups and in both conditions, but more so in French than in English in the verbal condition.

(1) Language properties influence verbal cognition, but do not seem to impact non-verbal cognition • Language differences occur in the production task that explicitly implies language use, but not in the categorization tasks (neither in verbal nor in non-verbal conditions).

(2) Manner and Path components are differentially accessible • Path is the main criterion chosen for categorization in both groups and in both conditions.

• However, relative focus on Manner depends on event type (boundary crossing > vertical; M salience) • Interactions also occur between event type, condition, and language.

(3) Methodological issues to take into account when testing language effects on non-verbal cognition (in progress) -Stimuli: It is necessary to use more ecological motion (humans, videos) and to balance the salience of Path and Manner.

-Measures: It is necessary to test on-line processes of attention allocation (eye-tracking). In addition, other interactions show that:

-INTO/OUT-OF elicit most manner choices in both conditions and in both languages.

-ALONG elicits most Manner choices in the verbal condition.

-UP/DOWN elicit Manner choices in English in both conditions but only in the verbal condition in French.

-In both languages ACROSS elicits Manner choices in the verbal condition but less so in the non-verbal condition in French.

Speakers' choices of Manner vs. Path criteria depended on conditions and stimuli properties (Manner and Path).

• Conditions (Fig. 3) Overall, both groups rely less on Manner in both conditions, notwithstanding two tendencies (not significant):

• English group ≥ French group;

• non-verbal ≥ verbal condition (English group).

INTRODUCTION

  Fig. 1. Stimuli: target and choices. Triads were used in categorization, targets only in production

Fig. 4 .

 4 Fig. 4. Manner choices as a function of Manner (a) and Path (b) in stimuli (collapsing conditions)

  Fig. 5. Manner choices by Path type & Condition

Fig. 3 .

 3 Fig. 3. Manner choices in verbal vs. non-verbal conditions