Alexandre Krupa

A new duty-cycling approach for 3D needle steering allowing the use of the classical visual servoing framework for targeting tasks

We present a method for three dimensional steering of a beveled-tip flexible needle that can be used in medical robotics for percutaneous assistance procedures. The proposed solution is based on the design of a new duty-cycling control strategy that makes possible to control three degrees of freedom (DOFs) of the needle tip. These DOFs correspond to two rotations that are performed around the needle tip lateral axes and one translation along the needle shape axis. This approach has two advantages compared to existing methods. First it does not rely on a trajectory planning and opens therefore numerous closed-loop control scheme possibilities as for example the implementation of the visual servoing framework based on the task function approach. Second, it does not constraint the needle to follow a succession of planar arcs but allows non planar 3D trajectories. Experimental results of a targeting task performed by visual servoing demonstrate the feasibility of this new concept and its robustness to needle kinematic model errors.

I. INTRODUCTION

Physicians perform needle insertion into human tissues for various medical interventions such as biopsy or injecting therapeutic drugs into local tumours. This kind of operations requires a lot of skills from the clinician in order to reach the target successfully. In fact, positioning errors with respect to the target can imply injuries of tissues and compromise the effectiveness of diagnosis or therapy, or bias the results of studies that require accurate needle placement. That is the reason why robotic assistance for needle steering is a fascinating subject of ongoing research as reported in [START_REF] Reed | Robot-assisted needle steering[END_REF]. The needles can be either rigid or very flexible depending on the medical applications considered. In practice, all needles tend to bend when they are inserted due to the forces that the tissues exert on them. This bending effect is strongly increased for flexible needle and becomes an advantage for targeting anatomical elements that were not possible to reach with a standard rigid needle. In fact, bending the needle can be very useful in order to avoid obstacles such as sensitive organs or bones. The material of this type of needle is flexible and its tip is beveled so that the tissues exert an asymmetrical force on the needle tip causing it to bend naturally. In opposite, the tip of a standard rigid needle is symmetric, and the material used is generally stiff. A flexible needle with asymmetric beveled tip, when inserted straight forward follows a circular arc trajectory [START_REF] Webster | Nonholonomic modeling of needle steering[END_REF], [START_REF] Misra | Mechanics of flexible needles robotically steered through soft tissue[END_REF]. The trajectory taken depends on the geometry of the beveled tip and on the relative stiffness between the needle and the tissue. In order to steer the needle, it is needed to use a kinematic model.

Alexandre Krupa is with Inria Rennes -Bretagne Atlantique and IRISA, France, alexandre.krupa at inria.fr Generally the model considered in most of research studies related to needle steering is the one of a bicycle or a unicycle [START_REF] Webster | Nonholonomic modeling of needle steering[END_REF]. An alternative approach was also proposed in [START_REF] Glozman | Image-guided robotic flexible needle steering[END_REF] where the needle is modeled as a linear beam supported by virtual springs. In all existing approaches a path planning initial step is required to predict the needle trajectory before steering it. The planner must take into account the nonholonomic constraints of the needle in order to work. Different solutions were proposed to reach this goal. Most approaches use inverse kinematics to plan paths in 2D or 3D [START_REF] Duindam | Threedimensional motion planning algorithms for steerable needles using inverse kinematics[END_REF]. In [START_REF] Alterovitz | Planning for steerable bevel-tip needle insertion through 2d soft tissue with obstacles[END_REF] a finite element mesh is used to predict soft tissue deformations and the needle path planning is done thanks to numerical optimization. Probabilistic methods were also tested, like the Rapidly-Exploring-Random Tree (RRT) algorithm [START_REF] Lavalle | Randomized kinodynamic planning[END_REF], or the Stochastic roadmap [START_REF] Alterovitz | The stochastic motion roadmap: A sampling framework for planning with markov motion uncertainty[END_REF]. This last one gives very good results for preoperative path planning, however it requires a lot of computing and therefore cannot be used to adapt the path in real-time intraoperatively.

Adjusting the path intraoperatively in a control-loop has been done in [START_REF] Wood | Needle steering system using duty-cycled rotation for percutaneous kidney access[END_REF] for needle steering in a 2D plane thanks to the use of an extented Kalman filter that estimates the current pose of the needle from image feedback. The path planning consists then in computing a planar arc-trajectory of the needle that links the current tip position to the target thanks to the use of the inverse kinematic model of the needle. This path is then used to determine the input of a duty-cycling control technique [START_REF] Minhas | Modeling of needle steering via duty-cycled spinning[END_REF] that is applied to steer the needle along the desired arc-trajectory. Other approaches based on the RRT algorithm were proposed to take into account obstacles during the intraoperative path planning. For example, a RRT algorithm with backchaining was presented in [START_REF] Xu | Motion planning for steerable needles in 3d environments with obstacles using rapidlyexploring random trees and backchaining[END_REF] to steer the needle in 3D with obstacle avoidance. More recently, a combination of the RRT algorithm and duty-cycling control technique was proposed in [START_REF] Bernardes | Robot-assisted automatic insertion of steerable needles with closed-loop imaging feedback and intraoperative trajectory replanning[END_REF] for needle steering in a 2D plane. In this last work the 2D path planning is updated from the current needle tip position. This later is automatically tracked in an image provided by a camera observing the needle evolving in a translucent phantom.

In this paper, we propose a new approach to steer a needle in 3D that has several advantages. First, it does not require a trajectory planning that can be time consuming. Second, it does not constraint the needle to follow a succession of planar arcs but allows non planar 3D trajectories. Third, with the proposed method the direct control of the lateral angular velocities ω x , ω y and the insertion velocity v z of the needle extremity is possible. This last point is crucial since, in opposite to path planning approaches, it opens numerous control scheme possibilities as for example visual control using the classical framework of visual servoing [START_REF] Chaumette | Visual servo control, part i: Basic approaches[END_REF] that is based on the task function approach [START_REF] Samson | Robot control: The task function approach[END_REF]. This work is a preliminary study toward the development of a needle steering robotic system using 3D ultrasound feedback. In our previous work [START_REF] Nadeau | Active stabilization of ultrasound image for robotically-assisted medical procedures[END_REF], in collaboration with the authors of [START_REF] Bernardes | Robot-assisted automatic insertion of steerable needles with closed-loop imaging feedback and intraoperative trajectory replanning[END_REF], we already addressed the compensation of the patient physiological motion by stabilizing the ultrasound image thanks to an ultrasound visual servoing technique [START_REF] Nadeau | Intensity-based ultrasound visual servoing: modeling and validation with 2d and 3d probes[END_REF]. We also recently proposed in [START_REF] Chatelain | Real-time needle detection and tracking using a visually servoed 3d ultrasound probe[END_REF] an efficient image processing method to automatically detect and track in real-time a flexible needle in a sequence of ultrasound volumes.

The paper is organized as follows. Section II presents the kinematic model of the flexible needle considered in our study and section III recalls the principle of the classic dutycycling control technique allowing to bend the needle with a desired curvature in a plane. We present in Section IV our new 3D duty cycling strategy that allows to control the two lateral angular velocities and the insertion velocity of the needle. Section V presents the control scheme we propose for automatic soft tissue targeting using visual servoing. Finally, experimental results on the visual targeting task are reported and discussed in section VI.

II. NEEDLE KINEMATIC MODEL

A. Needle natural bending

The bevel-tip needle can be driven by applying two velocity inputs to the needle base. These velocities are the insertion speed and the angular spin velocity. The rotational velocity causes the needle to turn about its shaft, reorienting the bevel of the tip. If we neglect the torsional compliance of the needle, the tip rotates at the same speed as the needle base. As previously mentioned, when only a translation motion is applied to the needle, the soft tissue exerts asymmetric forces on the beveled-tip that deflect the needle trajectory. As a consequence the needle follows an arc in a 2D plane. Moreover, if the needle is appropriately selected in terms of material compared to the tissue mechanical properties (needle less stiff than the tissue), the needle shaft trajectory will be the same as the needle tip. Thus, tracking the position of the needle can be performed by tracking only the needle tip.

The needle can be modeled as a nonholonomic vehicle owning 2 degrees of freedom DOF (forward translation and rotation around its shaft). In fact, as experimentally demonstrated in [START_REF] Webster | Nonholonomic modeling of needle steering[END_REF], inserting the needle straight forward in soft tissue is equivalent to ride an unicycle along a planar circular arc of constant radius r = 1 K as illustrated in Figure 1. K is the natural curvature of the needle which depends both on needle and tissue mechanical properties. Applying a rotation around the needle shaft allows to reorientate the plane where the needle will further evolve.

B. Nonholonomic constraints

Let F n and F w be respectively a Cartesian frame attached to the needle tip and the world frame as shown in Figure 1. We describe the 3D pose p (with p ∈ SE(3)) of the needle tip r with respect to F w by the homogeneous 4×4 transformation matrix :

-! y n -! z n -! x w -! y w -! z w ω z v z F n F w
w M n = w R n w T n 0 [1×3] 1 (1)
where w R n is a 3 × 3 rotation matrix representing the orientation of F n with respect to F w and w T n is the 3 × 1 translation vector defining the origin of F n in F w . We also denote the instantaneous velocity of the needle frame F n with respect to the world F w by the velocity screw vector:

(n) v n = ( (n) v x , (n) v y , (n) v z , (n) ω x , (n) ω y , (n) ω z ).
Note that the upper script (n) indicates that this velocity is expressed in the frame F n . Here (n) v x , (n) v y , (n) v z are the translation velocities along each axis of F n , and (n) ω x , (n) ω y , (n) ω z are the angular velocities around each axis. When applying a forward insertion motion to the needle, this later can only move in the plane y-z of F n . Thus, the velocity components (n) v x , (n) v y , (n) ω y are constraint to be null and only the insertion velocity (n) v z , and angular spin velocity (n) ω z can be applied by the user to steer the needle. The bending effect introduces however a rotation velocity around the x axis of F n that depends on the insertion velocity and natural curvature of the needle such as

ω x = K (n) v z .
The screw velocity vector is therefore linked to the needle input velocities by the following kinematic relation that expresses the nonholonomic contraints:

(n) v n = 0 0 1 K 0 0 0 0 0 0 0 1 T (n) v z (n) ω z (2)

III. DUTY-CYCLING FOR STEERING INSIDE A PLANE

A. Principle

The duty cycling method [START_REF] Minhas | Modeling of needle steering via duty-cycled spinning[END_REF] was introduced to allow the proportional control of the curvature of a bevel-tip needle during insertion. It works by alternating between time period T trans of pure insertion where the needle will follow an arc of maximum curvature K, with time period T rot of both insertion and rotation, where the needle will follow a helical curve. This later becomes a straight line if the rotation velocity is highly superior than the insertion velocity. Those two phases compose a cycle of period T = T rot + T trans .

The ratio of the rotation period T rot to the cycle period T is defined by:

DC = T rot T (3)
With this technique the effective curvature K ef f that will take the needle is given by the following formula:

K ef f = K(1 -DC) (4)
Therefore, it is possible to set the curvature of the needle between 0 and the natural curvature K. In previous works [START_REF] Minhas | Modeling of needle steering via duty-cycled spinning[END_REF], [START_REF] Bernardes | Robot-assisted automatic insertion of steerable needles with closed-loop imaging feedback and intraoperative trajectory replanning[END_REF] this method has been implemented by moving during each cycle the needle on a constant fixed insertion distance ∆z, and by applying a constant rotational velocity ω z = ω c . As a result, the cycle can be defined by:

T rot = 2π ω c , T trans = T rot (1 -DC) DC , v z = ∆z T (5)

B. Discussion on the limitations

With this technique the needle does a complete rotation (2π rad) during each cycle of period T which allows it to stay in the working plane. However, using this technique restricts the needle to evolve only in a plane. To overcome this limitation some colleagues proposed in [START_REF] Bernardes | 3d robust online motion planning for steerable needles in dynamic workspaces using duty-cycled rotation[END_REF] to steer the needle along a 3D trajectory composed from a succession of 2D planar arcs. The principle is to use a path planner based on a RRT algorithm that pre-computes from the current needle position and the target position a succession of planar arcs with their relative orientation. The needle is then steered by a stop and turn technique which consists at the starting point of each new arc to first stop the needle and rotate it about its z-axis. This allows to re-orientate the needle beveled tip in the new planned working plane before applying the duty cycling technique. However, this approach is computationally expensive since it requires to perform an initial path planning and update it during the process in order to compensate the model errors . Moreover, the trajectory is constrained to follow a succession of planar arcs that may result to important bends at the arc junction points if the planning is not frequently updated.

IV. 3D DUTY-CYCLING NEW STRATEGY

A. Method for controlling the needle velocities ω x , ω y , v z

The principle relies on a modification of the duty-cycling technique in order to allow continuous control of the velocity components ω x , ω y and v z of the needle. Let us consider a new cartesian frame F ñ attached to the needle that corresponds to the current needle frame rotated with an angle θ around the z axis of F n . By expressing the rotational velocity of F ñ in F n and using the relation (2) with the effective curvature K ef f we obtain the following relation that links the angular velocity of F n to the needle input velocities (ñ) v z and (ñ)

ω z expressed in F ñ:   (n) ω x (n) ω y (n) ω z   =   cos(θ) -sin(θ) 0 sin(θ) cos(θ) 0 0 0 1     K ef f (ñ) v z 0 (ñ) ω z   (6)
As F ñ and F n have the same origin and their z-axis are aligned, we have (n) v z = (ñ) v z that we can substitute in [START_REF] Alterovitz | Planning for steerable bevel-tip needle insertion through 2d soft tissue with obstacles[END_REF] in order to formulate the angular velocity of F n in function of the needle translational velocity and effective curvature:

(n) ω x = cos(θ)K ef f (n) v z (7) (n) ω y = sin(θ)K ef f (n) v z
One can notice that it is therefore possible to control the ratio between the velocities (n) ω x and (n) ω y at each cycle period T by applying a given value of angle θ around its z-axis before inserting the needle. Section IV-B will describe how this angle accumulation is performed in practice during the duty-cycling. Now, we detail the computation of the desired curvature value K * ef f , the desired angle θ * that represent with (n) v * z the control inputs of the duty cycling low-level controller. We denote by

( (n) v * z , (n) ω * x , (n) ω * y )
the desired control velocity values that we want to apply to the needle frame. Dividing the second line by the first line of ( 7) leads directly to :

θ * = atan( (n) ω * y / (n) ω * x ) (8)
and from simple calculus on (7) we obtain:

K * ef f = (n) ω * 2 x + (n) ω * 2 y / (n) v * z (9)
Note that if (n) v * z is null we set K * ef f = 0 since the needle is not moving in this case.

B. Duty cycling low-level controller with angle accumulation

The duty-cycling classic technique introduced in section III-A is not so easy to implement in practice since it is not possible to accurately apply a constant rotation velocity (n) ω z = ω c during a given time T rot and stop immediately the rotation when the spin angle reaches 2π rad before beginning the pure translation phase T trans . In fact, performing like this would require a robotic system with a rotating actuator that provides an infinite acceleration to start and stop the needle rotation without generating transient phases. Therefore to cope with this important issue, we propose in our new 3D duty-cycling approach to implement a positionbased controller that regulates the needle spin angle in order to accurately follow a desired angular trajectory. It is based on a motion profile generator that computes in real-time a trapezoidal acceleration profile that is applied to the PID controler of the actuator rotating the needle. The advantage of this position-based controller is to avoid any angular drift on the rotation around the z-axis. In fact, this is crucial for accurately controlling the needle spin angle during the duty cycling process. With this approach, the desired rotation velocity (n) ω * z is not constant during the acceleration and deceleration phases of the rotating actuator but an accurate final angle positioning is achieved.

Here we detail the successive steps that are performed during one period of our modified duty-cycling approach that allows the application of the velocity control inputs

( (n) v * z , (n) ω * x , (n) ω * y )
to the needle frame F n . 1) measure the current absolute needle spin angle: θ a 2) compute the desired absolute spin angle to achieve by accumulating 2π rad and the control angle input θ * given by eq. ( 8) : θ * a = θ a + 2π + θ * 3) compute the desired period T * rot and desired spin angular velocity profile (n) In practice, we also add a saturation mechanism on the desired curvature value K * ef f in order to guarantee that the pure translation phase T trans will not take infinite time. In fact this could happen when the desired curvature is not reachable due to the limitation of needle natural curvature K. We propose therefore to saturate K * ef f to 80 % of the K value.

V. AUTOMATIC TARGETING BY VISUAL SERVOING

In this section we present the visual control law that we propose for automatically guiding the needle in order to reach a target point O defined by the user at location w T o = (t ox , t oy , t oz ) in the world frame.

As the target coordinates are expressed in the world frame F w and not directly in the needle frame F n , it is more suitable to express the needle lateral angular control velocities in the fixed world frame F w . The lateral angular velocities applied to the needle frame will then be obtained from the lateral control velocities ( (n) ω *

x , (n) ω * y ) expressed in F w thanks to the following relation:

(n) ω * x (n) ω * y = w R n [2×2] T (w) ω * x (w) ω * y (10) where w R n [2×2]
T is the transpose of the [2 × 2] upper-left part of w R n . It follows that in order to control the needle in the world frame reference, we need to know the rotation matrix that describes the orientation of the needle frame with respect to the world frame. We will detail now the method we propose to estimate it and also how to reconstruct the full needle 3D pose from measurements.

A. Needle 3D pose estimation

We consider that a current 3D volume observing the needle is available thanks to the use of a 3D ultrasound probe or a TDM or MRI imaging device. It is therefore possible to extract the voxels belonging to the shape of the needle by performing some image processing on the volume. We will not detail the method to extract the needle but the reader can refer to our previous work [START_REF] Chatelain | Real-time needle detection and tracking using a visually servoed 3d ultrasound probe[END_REF] where we presented a RANSAC-based algorithm to extract in real-time a flexible needle from a sequence of 3D ultrasound volumes. This extraction gives us the current location of the needle tip point denoted N of coordinates w T n = (t nx , t ny , t nz ) and its unitary orientation vector w u = (u x , u y , u z ) at its extremity both expressed in the volume cartesian frame that corresponds to the world frame F w . From these visual measures we can therefore only partially determine the current 3D pose of the needle frame with respect to the world frame that is described by the homogeneous matrix w M n (see [START_REF] Reed | Robot-assisted needle steering[END_REF]). In fact, the missing information corresponds to the orientation of the needle around its z-axis since it is not possible to identify it from vision due to the tubular needle shape. However, this angle information is directly accessible from the duty-cycling low-level controller and corresponds to the absolute angle θ a that is continuously measured by the actuator encoder. By describing the orientation of the needle w R n with the Euler angles using the x-y-z convention where (φ, ψ, γ) are respectively the rotation angles around the x, y and z axis we obtained the following expression:

w R n =   c ψ c γ -c ψ s γ s ψ s φ s ψ c γ + c φ s γ -s φ s ψ s γ + c φ c γ -s φ c ψ -c φ s ψ c γ + s φ s γ c φ s ψ s γ + s φ c γ -c φ c ψ   (11)
with the notation c x = cos(x) and s x = sin(x). As the third column corresponds directly to the needle unitary orientation vector w u extracted from the volume we can easily identify the φ and ψ angles by the relations:

φ = atan(-u y /u z ) ψ = atan(u x /(c φ u z -s φ u y )) (12)
If the needle frame axes are initially aligned with the world frame axes before starting the 3D duty-cycling control process, then the angle γ will correspond directly to the absolute angle θ a measured by the actuator encoder during the needle trajectory evolution. Therefore by substituting γ = θ a and ( 12) in [START_REF] Xu | Motion planning for steerable needles in 3d environments with obstacles using rapidlyexploring random trees and backchaining[END_REF] we obtain the estimated rotation matrix w R n . Consequently we also obtain the full homogeneous matrix w M n that describes the estimated current pose of the needle frame with respect to the world frame.

B. Visual servoing

To perform the task we propose to apply a constant needle insertion velocity (n) v z = v c and online adapt the angular control velocities ( (w) ω *

x , (w) ω * y ) by visual servoing in such a way to orientate the needle z-axis toward the target. The set of visual features that we chose is s = (l xz , l yz ) where l xz and l yz are respectively the distances between the target point O and their orthogonal projections onto the needle line direction projected in the x-z plane and y-z plane of F w as shown in Figure 2a. Let us denote d xz and d yz the distances of the projections of vector --→ N O into the xz plane and y-z plane of F w . We also define θ x and θ y as being respectively the angles around the world frame x and y axes that describe the relative orientation of the vector --→ N O with respect to needle unitary orientation vector w u (see Figure 2a). Note that all these geometrical features can be easily computed from the measured 3D positions of the needle and target. The current visual features are then obtained by: l xz = d xz sin(θ y ) and l yz = -d yz sin(θ x ). The visual servoing task consists then in aligning the needle line direction on the target by regulating s to the desired value s * = (0, 0). As we previously mentioned, one important advantage of our new 3D duty-cycling strategy is that we can apply the classical framework of visual servoing [START_REF] Chaumette | Visual servo control, part i: Basic approaches[END_REF]. Therefore we first determine the image jacobian J that relates the variation ṡ of the selected visual features to the needle velocity inputs and obtain the following expression from basic kinematic relationships:

ṡ = lxz lyz = J ωxy J vz =J   (w) ω x (w) ω y (n) v z   ( 13
)
with

J ωxy = 0 d xz cos(θ y ) -d yz cos(θ x ) 0 , J vz = 0 0 (14)
We can note that the variation of the features are totally decoupled from the needle translational velocity applied in the needle frame and that they are independently coupled to the needle lateral angular velocities expressed in the world frame. This shows that the selected features are strongly suited for the required target alignment task since the constant translation velocity will note disturb the needle alignment task. Finally we applied the following classical control law that gives the desired control angular velocity inputs in such a way to obtain an exponential decrease of the visual error e = ss * :

(w) ω *

x

(w) ω * y = -λJ -1 ωxy e (15)
where λ > 0 is a proportional gain. The angular velocities are then expressed in the needle frame F n thanks to the relation [START_REF] Minhas | Modeling of needle steering via duty-cycled spinning[END_REF]. The desired needle spin angle θ * to accumulate and the desired effective curvature K * ef f are then computed using ( 8)-( 9) and applied with the constant insertion velocity to our 3D duty-cycling low-level controller presented in Section IV-B.

VI. EXPERIMENTAL RESULTS

In this section we present experimental results of the automatic visual targeting task using our 3D duty-cycling new technique. We consider the images provided by two cameras in order to demonstrate the feasibility of our approach but other kinds of image modality as for example intraoperative 3D ultrasound images could be employed. The experiment setup is presented in Figure 2b. It is composed of a 6-DOF anthropomorphic robot holding a 22 gauge bevel-tip flexible needle whose length is 20 cm. A home-made translucent gelatin phantom contained in a translucent box is employed to simulate a piece of organic soft tissues and is observed by two orthogonal video cameras. These cameras are connected on the firewire interface of a PC equiped with a Dual-core 2.4 Ghz Intel Pentium that grabs the images at a frame rate of 15 images/second. A flat screen monitor is positioned behind the phantom in order to provide a scattered light inside it. Another firewire connexion of the PC is used to send the joint velocity control inputs to the robot and receive the joint measures. A software that we designed in C++ (based on ViSP open source library [START_REF] Marchand | Visp for visual servoing: a generic software platform with a wide class of robot control skills[END_REF]) is running on the PC to perform the needle motion control. This software realizes in a a fist thread the image processing and computes the visual servoing control law output [START_REF] Nadeau | Active stabilization of ultrasound image for robotically-assisted medical procedures[END_REF]. This later is continuously sent to a second thread that performs the 3D duty-cycling low-level process.

F w F w -→ x w -→ y w -→ z w -→ z w θ y θ x O O N N needle needle target target Front view X-Z Side view Y-Z d xz d yz l yz l xz (a)
The intrinsic parameters of the two cameras were calibrated in order to allow the conversion from pixels to meters for a given distance to the phantom. In order to observe the whole trajectory of the needle, these distances were respectively fixed to 14 cm and 19 cm for the left and right cameras since their objective focal length were not identical. To perform a 3D localization of the needle during its insertion, we extract the needle tip meter coordinates from the two orthogonal images thanks to an image processing algorithm that automatically detects the shape of the needle. This algorithm consists in thresholding the images and applying few morphological operators (erosion and dilatation). For a lack of place we do not detail it more in depth since it is out of the topic of this paper. For all experiments the needle frame F n was initially aligned with the frame F w before the insertion.

In a first experiment, we do not apply our visual servoing scheme that relies on our 3D duty-cycling technique but we apply only a constant insertion velocity (n) v * z = 1 mm/s to the needle. Figures 3a and3b show respectively the resulted trajectory observed in the left image (corresponding to the y-z plane of F w ) and the right image (corresponding to the x-z plane of F w ). One can observe that the needle follows as expected its natural curvature. Note that the needle current location is surrounded in red by the image processing algorithm and that the other visible lines correspond in fact to few holes performed by the needle in the phantom during previous insertions. Moreover, the final vertical locations in pixels are not the same in the two images since the distances from the cameras to the phantom are different but the final vertical coordinates in meter that are obtained after pixel-meter conversion (thanks to the calibrated intrinsic parameters of the cameras) are identical. From this first experiment, we coarsely estimate from the location of three 3D points on the needle arc trajectory the natural curvature of the needle and obtain K = 2 meter -1 (corresponding to a radius of 0.5 meter).

In a second experiment, we test our targeting task by applying the visual control law (15) and using our 3D dutycycling low-level controller with a constant insertion velocity of 1 mm/s. The nominal needle spin velocity that we use to compute the profile of the spin velocity (n) ω * z that is applied by our 3D duty-cycling low-level controller (see Section IV-B) was fixed to ω c = 2π rad/second. For this experiment the control gain was set empirically to the low value of λ = 0.02 in order to obtain a smooth trajectory. Since we did not put a physical target in the phantom, we first inserted the needle manually in order to reach a position that we defined as being the target to reach. This allows us to set the target coordinates by clicking on the needle tip in the two images. Then we remove the needle and slightly displace the phantom in order to insert the needle in a new entry point of the phantom. Finally we apply the automatic targeting by visual servoing. Figures 3c and3d x and (w) ω * y of the needle expressed in the world frame are reported in Figures 4c-4d. One can see that the velocities decrease during the first 70 seconds and then increase when the needle reaches the target. This is due to the fact that the control law actively reacts to compensate the final increase of angles θ x and θ y when the needle tip get more and more close to the target in order to refine the alignment, but this is not an issue since the needle insertion motion is automatically stopped when the z component of the euclidean distance between the tip and target is null. Figures 4e-4f show the projection of the needle trajectory onto the y-z and x-z planes of F w and the 3D trajectory is depicted in Figure 5a. Again, we can see that the target depicted by the red circles is accurately reached as it is also demonstrated in the Figure 5b that presents the evolution of the needle-target distance. The evolution of the duty-cycling control ratio DC * applied by our 3D duty-cycling low-level controller is reported in Figure 5c. We can observe that at the start of the targeting task DC * is low in order to apply an important curvature to the needle, then in the period from 40 to 70 seconds, it progressively increases to realize a more straight trajectory of the needle since this one becomes well aligned on the target. Then the ratio DC * drastically decreases when the needle reaches the target due to the final increase of θ x and θ y that we explained before. Figure 5d shows the evolution of the desired spin angular velocity profile (n) ω * z applied by the duty-cycling low-level controller. We can see that the duration T trans of the pure insertion motion (phases when (n) ω * z = 0) progressively decreases while the needle direction becomes more and more aligned on the target. These results also experimentally demonstrate that the targeting task by visual servoing is robust to modelling errors since we do not consider the torsional effect of the needle in its kinematic model and we employ a coarse estimated value of its natural curvature that is not on-line adapted during the control.

VII. CONCLUSION

In this paper we presented a new duty-cycling strategy that allows to control the insertion velocity and the two lateral angular velocities of the needle tip. This control is performed continuously at the rate of the duty-cycling period and allows therefore to steer the needle in 3D in opposite to the classical duty-cycling approach that is limited to planar bending. Experimental results demonstrate the efficiency of our approach to steer the needle to reach a 3D target without relying like other methods on a path planning that generates a 3D trajectory composed from successive oriented planar arcs. In fact, as our strategy allows to implement the classical task function approach, it will be possible in the future to design a multi-task control law to perform several hierarchical tasks at the same time like for example reaching a soft-tissue target while insuring obstacle avoidance and needle mechanical constraints. We are currently implementing the approach presented in this paper to automatically steer a flexible needle by ultrasound visual servoing thanks to the use of a 3D ultrasound probe.

Fig. 1 .

 1 Fig. 1. Needle beveled-tip modeled by an unicycle following an arc of radius r.

  ω * z needed to achieve θ * a based on the trapezoidal acceleration profile with a nominal needle spin velocity of ω c 4) compute the desired insertion distance ∆z * = (n) v * z T * rot 5) perform the rotation and insertion thanks to the positionbased robot controller to reach the desired angle θ * a and translation ∆ * z 6) measure the real period T rot that was taken by the system 7) compute the desired cycle ratio DC * = 1 -K * ef f /K with K * ef f given by (9), the desired pure translation period T * trans = T rot (1 -DC * )/DC * and the desired insertion distance ∆z * = (n) v * z T * trans 8) perform the insertion thanks to the position-based robot controller to reach the desired translation ∆ * z 9) go back to step 1) to start a new cycle.

Fig. 2 .

 2 Fig. 2. (a) Geometry of selected visual features. (b) Experimental setup.

Fig. 3 .

 3 Fig. 3. (a) left camera (y-z plane) and (b) right camera (x-z plane) views of the needle trajectory obtained with a pure needle insertion motion. (c) left camera (y-z plane) and (d) right camera (x-z plane) views of the needle trajectory obtained with the targeting task by visual servoing. The target is represented by the red cross and needle tip by the green one.

  present the

Fig. 4 .

 4 Fig. 4. Results of the targeting task by visual servoing: (a)-(b) visual features errors, (c)-(d) control angular velocities expressed in the world frame, (e)-(f) obtained needle trajectory in the y-x and x-z world frame planes. The red circle represents the target to reach.obtained trajectory that is observed in the left and right images. One can see that the target depicted by the red cross is perfectly reached by the needle tip in the two images. The evolutions of the visual errors are presented in Figures4a-4band correspond as expected to exponential decreases toward zero. Note that the noise is essentially due to the image processing algorithm and could be reduced with the use of a more efficient algorithm to extract the needle from the images. At the convergence the final positioning errors are 0.4 mm for l xz and 0.5 mm for l yz . The lateral angular control velocities (w) ω *x and (w) ω * y of the needle expressed in the world frame are reported in Figures4c-4d. One can see that the velocities decrease during the first 70 seconds and then increase when the needle reaches the target. This is due to the fact that the control law actively reacts to compensate the final increase of angles θ x and θ y when the needle tip get more and more close to the target in order to refine the alignment, but this is not an issue since the needle insertion motion is automatically stopped when the z component of the euclidean distance between the tip and target is null. Figures4e-4fshow the projection of the needle trajectory onto the y-z and x-z planes of F w and the 3D trajectory is depicted in Figure5a. Again, we can see that the target depicted by the red circles is accurately reached as it is also demonstrated in the Figure5bthat presents the evolution of the needle-target distance. The evolution of the duty-cycling control ratio

  Duty-cycling velocity ω * z (deg/s) Fig. 5. Results of the targeting task by visual servoing: (a) obtained needle 3D trajectory, (b) evolution of the distance between the needle and target, (c) duty-cycling ratio and (d) spin angular velocity (n) ω * z applied by the low-level 3D duty-cycling controller.


M Romanelli

L Wang

M Brunel

M Vallet

Measuring the universal synchronization properties of coupled oscillators across the Hopf instability

Keywords:

When a driven oscillator loses phase-locking to a master oscillator via a Hopf bifurcation, it enters a bounded-phase regime in which its average frequency is still equal to the master frequency, but its phase displays temporal oscillations. Here we characterize these two synchronization regimes in a laser experiment, by measuring the spectrum of the phase fluctuations across the bifurcation. We find experimentally, and confirm numerically, that the low frequency phase noise of the driven oscillator is strongly suppressed in both regimes in the same way. Thus the long-term phase stability of the master oscillator is transferred to the driven one, even in the absence of phase-locking. The numerical study of a generic, minimal model suggests that such behavior is universal for any periodically driven oscillator near a Hopf bifurcation point.

I. INTRODUCTION

Synchronization, i.e. the tendency of coupled oscillators to oscillate at the same frequency, is recognized as "one of the most pervasive drives in universe" [START_REF] Strogatz | Sync : how order emerges from chaos in the universe, nature and daily life[END_REF][START_REF] Pikovsky | Synchronization : A Universal Concept in Nonlinear Sciences[END_REF]. Countless examples of this fact have been reported, from the synchronous blinking of fireflies to the coherent oscillation of coupled Josephson junctions [START_REF] Wiesenfeld | Synchronization transition in a disordered Josephson series array[END_REF] or of chemical micro-oscillators [START_REF] Toiya | Synchronization of chemical microoscillators[END_REF]. Synchronization plays also a role in memory processes [START_REF] Fell | The role of phase synchronization in memory processes[END_REF] and in the complex networks arising, for instance, in neurosciences or in economics [START_REF] Arenas | Synchronization in complex networks[END_REF]. Given that synchronization occurs in so many different contexts, the study of a specific system can acquire a broader significance by helping to identify features that are universal, i.e. common to a large class of different systems displaying synchronization.

In laser physics, issues related to synchronization are typically encountered in optically injected lasers. Optical injection has been studied for a long time, but it is still a subject of interest today, both for the great complexity and variety of the possible dynamics [START_REF] Wieczorek | The dynamical complexity of optically injected semiconductor lasers[END_REF][START_REF] Erneux | Laser Dynamics[END_REF], and for its importance for applications such as high-speed modulation [START_REF] Naderi | Modeling the injection-locked behavior of a quantum dash semiconductor laser[END_REF] or microwave photonics [START_REF] Hung | Optical doublesideband modulation to single-sideband modulation conversion using period-one nonlinear dynamics of semiconductor lasers for radio-over-fiber links[END_REF]. In particular, recently a synchronization regime of frequency locking without phase locking has been experimentally isolated in injected laser oscillators [START_REF] Kelleher | Phasor plots in optical injection experiments[END_REF][START_REF] Thévenin | Resonance assisted synchronization of coupled oscillators : frequency locking without phase locking[END_REF]. In this regime, the relative phase of two interacting oscillators is not constant in time, yet its variations are bounded, so that the systems keep oscillating at the same average frequency. It is quite obvious that this regime differs from the more common situation in which the oscillators are phase-locked, but how should one characterize and quantify this difference ? In the present paper, we propose to answer this question by studying experimentally and numerically a driven, optically-based radiofrequency (RF) oscillator. Furthermore, the bounded-phase behavior, predicted long ago in driven laser systems [START_REF] Braza | Constant phase, phase drift, and phase entrainment in lasers with an injected signal[END_REF][START_REF] Solari | Laser with injected signal : perturbation of an invariant circle[END_REF][START_REF] Yeung | Nonlinear dynamics of a solid-state laser with injection[END_REF] and more recently also in quantum dot lasers [START_REF] Pausch | Optically injected quantum dot lasers : impact of nonlinear carrier lifetimes on frequency-locking dynamics[END_REF][START_REF] Lingnau | Feedback and injection locking instabilities in quantum-dot lasers : a microscopically based bifurcation analysis[END_REF], is by no means specific to lasers, since it occurs also in models of circadian [START_REF] Kronauer | Mathematical model of the human circadian system with two interacting oscillators[END_REF] or coupled van der Pol [START_REF] Chakraborty | The transition from phase locking to drift in a system of two weakly coupled van der Pol oscillators[END_REF] oscillators. Actually, bounded-phase phenomena are extremely generic, since they are associated to a Hopf bifurcation [START_REF] Kelleher | Bounded phase phenomena in the optically injected laser[END_REF], and are to be expected when amplitude effects cannot be ignored in the dynamics of coupled oscillators [START_REF] Aronson | Amplitude response of coupled oscillators[END_REF]. Thus, it can be expected that our study of a driven opto-RF oscillator will be relevant for coupled oscillators in general. Indeed, we will show in the following that our system provides an example of a generic, universal behavior of driven oscillators close to a Hopf bifurcation point.

II. EXPERIMENTAL SETUP AND RESULTS

A. Experimental setup

The experimental setup (see Fig. 1) is similar to the one described in [START_REF] Thévenin | Phase and intensity dynamics of a twofrequency laser submitted to resonant frequency-shifted feedback[END_REF]. We build an optically-based RF oscillator, and drive it with a stable RF synthesizer, that is coupled to the slave oscillator by optical frequency-shifted feedback. The core of the experiment is a diode-pumped, solid-state dual-frequency laser (DFL) [START_REF] Brunel | Generation of tunable high-purity microwave and terahertz signals by two-frequency solid state lasers[END_REF], represented in Fig. 1(a). The laser cavity, of length L = 75 mm, is closed on one side by a high-reflection plane mirror, coated on the 5-mm long Nd :YAG active medium, and on the other side by a concave mirror (radius of curvature of 100 mm, intensity transmission of 1% at the lasing wavelength λ = 1064 nm). The active medium is pumped by a laser diode emitting at 808 nm. A 1 mm-thick silica étalon ensures single longitudinal mode oscillation. Two eigenmodes E x and E y , polarized along x and ŷ, with eigenfrequencies ν x and ν y respectively, oscillate simultaneously. An intracavity birefringent element (here two quarter-wave plates QWPs) induces a frequency difference, finely tunable from 0 to c/4L = 1 GHz by rotating one QWP with respect to the other [START_REF] Brunel | Tunable two-frequency lasers for lifetime measurements[END_REF]. The typical output power of the two-frequency laser is 10 mW when pumped with 500 mW. When the laser output is detected by a photodiode after a polarizer at 45 • , an electrical signal oscillating at the frequency difference ∆ν 0 = ν y -ν x is obtained. The DFL can thus be seen as an opto-RF oscillator.

In order to lock this oscillator to an external reference signal, we use optical frequency-shifted feedback [START_REF] Kervevan | Beatnote jitter suppression in a dual-frequency laser using optical feedback[END_REF] (Fig. 1(b)). The feedback cavity contains an acoustooptic modulator (AOM), driven by a stable RF synthesizer, which provides an external phase reference. Next, a quarter-wave plate at 45 • followed by a mirror flips the x and ŷ polarizations, and finally the laser beam is reinjected in the laser cavity after crossing again the AOM. As a result, a x-polarized field oscillating at the frequency ν y + 2f AO and a ŷ-polarized field oscillating at the frequency ν x + 2f AO are reinjected in the laser. We choose the value of 2f AO so that ν x + 2f AO is close to ν y . Under suitable feedback conditions, ν y locks to the injected beam frequency ν x + 2f AO , i.e. the frequency difference ν y -ν x locks to 2f AO . We note that the optical reinjection has no direct effect on E x , because the frequency difference between ν x and ν y + 2f AO is too large. For the same reason, multiple round trips in the feedback cavity have no effect on the dynamics. The laser output is detected with a photodiode (3 GHz analogic bandwidth) after a crossed polarizer, thus providing an electrical signal proportional to

I = |E x + E y | 2 .
The signal is then analyzed with an electrical spectrum analyzer, a digital P-Q signal analyzer and an oscilloscope.

B. Experimental results

An electrical signal, oscillating in the absence of feedback at the frequency difference ∆ν 0 = ν y -ν x around 200 MHz, is obtained by detecting the interference signal after a polarizer. The whole system can thus be seen as an opto-RF oscillator, driven by the RF synthesizer (master oscillator) when the optical feedback is present. The beatnote signal I can be visualized at each instant in time as a vector having two components P = I 0 (t) cos φ(t) and Q = I 0 (t) sin φ(t) in a reference frame rotating at the master signal frequency 2f AO . The phase φ(t) is defined as φ(t) = φ y (t) -φ x (t) -2π2f AO t, where φ x,y are the optical phases of the two electric fields E x and E y . For small enough detuning ∆ν = ∆ν 0 -2f AO between the slave and the master oscillator, phase-locking occurs, and the corresponding experimental phasor plot consists of just a fixed point (Fig. 2(a), regime I). If ∆ν is increased above the Adler frequency f A , the fixed point solution loses stability through a Hopf bifurcation. The point representing I then rotates at the frequency ∆ν on a limit cycle (Fig. 2(a), regime II). The relative phase exhibits bounded, periodic temporal oscillations, whose amplitude increases with ∆ν. Finally, when ∆ν is bigger than a second characteristic frequency f B , the origin of the plane lies inside the limit cycle and the phase variations are not bounded anymore (Fig. 2(a), regime III). In the time and frequency domains (Fig. 2(b-c)), phaselocking (I) appears as pure sinusoidal oscillation, and the corresponding power spectrum consists of a single, sharp Fourier peak. On the contrary, in the bounded (II) and in the unbounded-phase regimes (III), the power spectrum shows two peaks at the natural frequencies 2f AO and ∆ν 0 of the uncoupled master and slave oscillator respectively, and the time series displays a beating, i.e. a slow modulation of the amplitude of the oscillations. From Fig. 2, the question arises about the exact status of the bounded-phase regime, occurring when f A < ∆ν < f B . On the one hand, since the relative phase is bounded, the average frequency of the oscillators is the same, and they can be considered as synchronized. On the other hand, the bounded-phase regime appears very similar to the unbounded-phase regime. In particular, the power spectrum displays two separated peaks at the natural frequencies of the oscillators, as it happens when synchronization is absent. Furthermore, while a qualitative change of the dynamics occurs when passing from phase locking to the bounded phase through a Hopf bifurcation, in the {P, Q} plane in Fig. 2(a) the bounded-phase and unbounded-phase regimes appear topologically equivalent and no bifurcation is seen between them [START_REF] Wieczorek | The dynamical complexity of optically injected semiconductor lasers[END_REF].

I = |E x + E y | 2 between

the laser fields on a photodiode

To investigate quantitatively such questions, we propose to use the power spectral density (PSD), also called phase spectrum, of the relative phase between the oscillators as a meaningful measure of synchronization. The PSD measures the strength of the phase fluctuations of an oscillator at all frequencies, and thus characterizes its quality in a complete and rigorous way. Precisely, the PSD is defined as follows [START_REF] Rubiola | Phase noise and frequency stability in oscillators[END_REF][START_REF]IEEE Standard Definitions of Physical Quantities for Fundamental Frequency and Time Metrology[END_REF]. For a stochastic, ergodic process φ(t), the autocorrelation function is :

C(τ ) = φ(t)φ(t + τ ) = lim T →∞ 1 T T /2 -T /2 φ(t)φ(t + τ )dt.
(1) Brackets indicate the ensemble average. The singlesideband power spectral density of the process is defined as

P SD(f ) = 2θ(f ) +∞ -∞ C(τ ) exp (-i2πf τ )dτ, (2)
where θ(f ) is the Heaviside function. Indeed, the negative frequency components can be discarded because C(τ ) is a real, even function. P SD(f ) is frequently expressed in decibels as 10 log [P SD(f )], and its units are dBrad 2 /Hz. This is the quantity we have measured, and indicated as S φ (f ) = 10 log [P SD(f )]. In practice, the above definitions are not directly used for the computation of S φ (f ). One rather uses the equality

P SD(f ) df = 2θ(f ) |dΦ(f )| 2
, and estimates |dΦ(f )| 2 by taking the square modulus of the Fast-Fourier Transform (FFT) of a sampled time series of φ(t). In our setup, the relative phase φ(t) can be extracted from the time series of the two quadratures P (t) and Q(t). The spectrum S φ (f ) of φ(t) is then computed. Note that, when the phase is unbounded, the average slope of the drifting phase has to be substracted in order to get a stationary, ergodic process. The experimental spectra are shown in Fig. 3. One sees that, in the phase-locking regime, the phase fluctuations of the free running laser are strongly suppressed at low frequency. For instance, Frequency (Hz) S (f ) (dBrad /Hz) at 1 kHz from the carrier one has a noise reduction of more than 30 dB with respect to the free-running case. The phase of the slave oscillator is locked to the stable master, which ensures the long-term phase stability. The absolute value of the phase noise that we obtain, -100 dBrad 2 /Hz at 10 kHz from the carrier, is comparable to what is obtained with active stabilization techniques using electronic servo-loops [START_REF] Brunel | High-spectral purity RF beat note generated by a two-frequency solid-state laser in a dual thermooptic and electrooptic phase-locked loop[END_REF]. In the bounded-phase regime, it can be seen that the long-term phase stability is exactly the same as in the phase-locking regime. The low frequency part of the spectrum can hardly be distinguished in the two cases. In this case too, despite the occurrence of fast temporal oscillations, the average value of the slave oscillator phase is locked to the stable master phase, and low frequency drifts and fluctuations are removed with respect to the free running laser. The only differences in the spectra arise at high frequency. Indeed, in the bounded-phase regime one has a narrow peak at the detuning frequency ∆ν, and its harmonics. These peaks are the coherent signature of the Hopf bifurcation, and correspond to the onset of the deterministic oscillation of the phase beyond the bifurcation point f A . At the transition between the bounded and unbounded-phase regime, for ∆ν ≃ f B , there is an enormous increase of the phase noise. This can be understood by observing in Fig. 2(a) that the experimental limit cycle has a finite thickness, because of uncontrolled fluctuations of the parameters and of the setup. Therefore, when ∆ν approaches f B , erratic phase slips occur in a random fashion, driven by the noise present in the system. These stochastic, abrupt 2π phase jumps close to the point f B explain the huge amount of phase noise. Finally, the spectrum of the unbounded-phase regime corresponds essentially to the one of the free-running opto-RF oscillator, apart from the high-frequency peaks. Figure 4(a) reports the PSD value at 1 kHz from the carrier as a function of ∆ν. It is striking to compare this plot to a numerically calculated bifurcation diagram (Fig. 4(b)). Indeed, when looking at the long-term phase stability of the oscillator, the occurrence of a Hopf bifurcation at the point f A is completely undetectable. Conversely, at the point f B a sharp transition occurs as a signature of the unlocking of the slave oscillator. The conclusion that can be drawn is that the relevant boundary for the synchronization range is not the Adler frequency f A , but the bounded-phase frequency f B . From the viewpoint of dynamical systems, phase-locking and bounded-phase are two qualitatively different regimes, separated by a bifurcation ; yet this distinction turns out to be of little relevance if one is concerned with synchronization.

At this point, one may ask whether the bounded/unbounded-phase transition is a bifurcation or not. Mathematically, the answer actually depends on the choice of variables [START_REF] Wieczorek | The dynamical complexity of optically injected semiconductor lasers[END_REF]. If the beat-note signal is described using the amplitude I 0 (t) and the phase φ(t), then its phase space is the surface of a half-cylinder. Bounded-phase and unbounded-phase regimes correspond to limit cycles that do not wrap the half-cylinder (librations) or that wrap it (rotations) respectively, and are not topologically equivalent. On the contrary, if the signal is described using the quadratures P (t) and Q(t), the two regimes are topologically equivalent, as we have seen. This difference stems from the fact that the polar coordinates I 0 (t) and φ(t) are not homeomorphic with the planar coordinates P (t) and Q(t). So, mathematically the nature of the bounded/unbounded-phase transition is, to some extent, a matter of choice, and it is worth noticing that, in a theoretical bifurcation analysis, it looks more convenient to work with P (t) and Q(t), so that no bifurcation appears at f B [START_REF] Wieczorek | The dynamical complexity of optically injected semiconductor lasers[END_REF]. However, our results show clearly that the physics of the system is not ambiguous : a sharp transition occurs at the point f B , namely the transition from synchronization to desynchronization.

III. THEORY : LASER AND GENERIC OSCILLATOR MODELS

We compare the experimental observation with the following laser model [START_REF] Thévenin | Phase and intensity dynamics of a twofrequency laser submitted to resonant frequency-shifted feedback[END_REF] : where e x,y are the normalized amplitudes of the optical modes, m x,y describe the population inversions, η is the pump parameter, β accounts for the cross saturation in the active medium, γ is the feedback strength, ǫ is the inversion lifetime, ∆ = ∆ν/f R = (∆ν 0 -2f AO )/f R is the detuning, and f R is the relaxation oscillation frequency. The scaled time s is related to the physical time t by s = 2πf R t. The normalization has been chosen so to remove stiffness from the class-B laser equations [START_REF] Erneux | Laser Dynamics[END_REF].

de x ds = (m x + βm y ) 1 + β e x 2 , (3)
de y ds = (m y + βm x ) 1 + β e y 2 + i∆ e y + γe x , (4)
For more details on the model, and for the correspondence between the normalized and the physical quantities, see [START_REF] Thévenin | Phase and intensity dynamics of a twofrequency laser submitted to resonant frequency-shifted feedback[END_REF]. We have injected some noise in the system by replacing the pump parameter η with η(1 + 0.05ξ(s)), where ξ(s) is a δ-correlated, normally distributed stochastic process. We notice that pump fluctuations induce fluctuations of f R , and thus of the normalized quantities ∆ and γ. This fact must be taken into account when integrating the laser equations. The PSDs calculated from numerical integration of the laser model are reported in Fig. 5(a). As in the experiment, it is observed that optical feedback is effective in reducing the low frequency phase fluctuations, both when the phase is locked and when it is bounded. In both cases, the opto-RF oscillator is thus synchronized to the master, and the long-term phase stability is not reduced by the onset of self-sustained oscillations. The huge increase of phase noise around the point f B is also reproduced, as well as the fact that, in the unbounded-phase regime, the phase noise is essentially equal to that of the free-running oscillator. One can ask to which extent these findings are specific to our system.

In order to answer this question, we have numerically investigated the following minimal model [START_REF] Pikovsky | Synchronization : A Universal Concept in Nonlinear Sciences[END_REF] :

dy ds = (1 -i∆)y -(1 + iα)|y| 2 y + e. (6)
y represents the complex amplitude of a periodically driven oscillator, ∆ is the detuning between the master and the slave, α accounts for the possible dependence of the oscillator's phase on the amplitude of the oscillations, and e is the strength of the drive. Equation 6contains just a linear term accounting for the growth of oscillations, a saturating nonlinearity, and an external forcing, and is a universal equation describing any periodically driven oscillator close to a supercritical Hopf bifurcation point [START_REF] Pikovsky | Synchronization : A Universal Concept in Nonlinear Sciences[END_REF]. We have introduced some noise by replacing ∆ with ∆ [1 + 0.033 ξ(s)]. In the following, we have taken α = 0 because this seems more appropriate for comparison with a solid-state laser. However, some numerical simulations, not shown here, indicate that this parameter does not influence phase noise, when α 2 < 1/3. For higher values of α, the bifurcation structure of the generic model changes qualitatively [START_REF] Pikovsky | Synchronization : A Universal Concept in Nonlinear Sciences[END_REF], and we have not explored this case. The PSDs of the phase of y are plotted in Fig. 5(b), for different values of ∆. The agreement of the spectra in Fig. 5(b) with the laser model is striking. All the observed regimes are reproduced and appear thus to be generic for any driven oscillator undergoing a Hopf bifurcation.

IV. CONCLUSION

In conclusion, we have characterized the synchronization of an opto-RF oscillator to a reference oscillator by measuring the noise spectrum of the relative phase, and shown that the long-term stability of the master oscillator is transferred to the driven oscillator even when the latter loses phase-locking and enters the bounded-phase regime. The synchronization range thus extends up to the edge of the bounded-phase regime, and can be substantially larger than the phase-locking range (of a factor √ 2 typically [START_REF] Yeung | Nonlinear dynamics of a solid-state laser with injection[END_REF]). A generic, minimal model gives strong evidence that the robustness against a Hopf instability is a universal feature of synchronizing systems.

Given the ubiquitous occurrence of synchronization and the genericity of the mechanism discussed here, the present results can be relevant for many areas. In a broader context, recently synchronization phenomena have been highlighted in optomechanical [START_REF] Heinrich | Collective Dynamics in Optomechanical Arrays[END_REF] or micromechanical [START_REF] Agrawal | Observation of Locked Phase Dynamics and Enhanced Frequency Stability in Synchronized Micromechanical Oscillators[END_REF] oscillators. Most of the recently observed results could be explained in the framework of pure phase equations such as the Adler model, which is the case when the coupling is too weak to affect the amplitude. However, as soon as this condition is not fullfilled, Hopf dynamics is to be expected, therefore we believe that the results presented here will soon prove relevant also for these research fields. One can also expect the present results to have an impact in the context of delay-coupled systems at large [START_REF] Soriano | Complex photonics : Dynamics and applications of delay-coupled semiconductors lasers[END_REF]. In general, it is known that a delayed coupling can be a source of oscillations, and in some cases it prevents phase locking [START_REF] Erneux | Applied delay differential equations[END_REF]. However, frequency locking is more general than phase locking, and the oscillatory behavior of coupled, delayed systems is not incompatible with bounded-phase synchronization. Under our experimental conditions, the effect of the delay could not be observed, because the feedback is seen as instantaneous by the laser. Indeed, the typical time scale of the laser dynamics is 1/f R ≃ 10 µs, which is much longer than the external-cavity round-trip time τ ≃ 5 ns. On the contrary, the effect of delay could be conveniently studied using semiconductor lasers, which evolve on much faster time scales [START_REF] Sciamanna | Polarization synchronization in unidirectionally coupled vertical-cavity surface-emitting lasers with orthogonal optical injection[END_REF][START_REF] Ozaki | Leader-laggard relationship of chaos synchronization in mutually coupled vertical-cavity surface-emitting lasers with time delay[END_REF].

Figure 1 :

 1 Figure 1: (a) Dual-Frequency Laser. M1,2 : cavity mirrors. QWP : quarter-wave plate. (b) Experimental setup that allows synchronizing the DFL beat-note frequency to a RF synthesizer. DFL : Dual-frequency laser (slave oscillator). M : feedback mirror. QWP : quarter-wave plate. AOM : acoustooptic modulator. P : polarizer. D : detector. P-Q analyzer : digital vector signal analyzer, permitting to measure the quadratures of I with respect to the reference signal delivered by the RF synthesizer.

Figure 2 :

 2 Figure 2: (a) Experimental phasor plots of the output signal I = |Ex + Ey| 2 , in the reference frame of the local oscillator at 2fAO. Each plot contains 2000 points, recorded over 100 µs. (I) Phase-locking regime. (II) Bounded-phase regime. (III) Unbounded-phase regime. (b) Corresponding experimental time series. (c) Power spectra.

φ 2 Figure 3 :

 23 Figure 3: Measured phase spectra S φ (f ) for different values of the detuning ∆ν. Black : ∆ν = 60 kHz, phase locking regime. Blue : ∆ν = 220 kHz, bounded-phase regime. Violet : ∆ν = 240 kHz ≃ fB. Red : ∆ν = 260 kHz, phase-drifting regime. Green : freerunning laser (no feedback from the external cavity).

Figure 4 :

 4 Figure 4: (a) Measured phase spectra S φ (f ) at 1 kHz from the carrier frequency, as a function of the detuning ∆ν. The arrows correspond to the spectra of Fig. 3 (except for the freerunning laser situation, not represented here). The red line is an eye-guide. (b) A phase bifurcation diagram, calculated using the laser model and the measured values of fA and fB, is shown for comparison.

dm x,y ds = 1 -

 1 |e x,y | 2 + β|e y,x | 2 + (5) -ǫ m x,y 1 + (η -1)(|e x,y | 2 + β|e y,x | 2 ) ,

Figure 5 :

 5 Figure 5: (a) Simulated phase spectra S φ (f ) using the laser model. Black : ∆ν = 100 kHz, phase locking regime. Blue : ∆ν = 180 kHz, bounded-phase regime. Violet : ∆ν = 230 kHz ≃ fB. Red : ∆ν = 400 kHz, phase-drifting regime. Green : free-running laser (no feedback from the external cavity). The values of the parameters used in the simulations are : β = 0.6, fA = 160 kHz, fR = 70 kHz, γ = fA/fR, ǫ = 0.97 10 -2 , η = 1.2. (b) Simulated phase spectra S φ (f ) for the generic model. Black : e = 1.7, ∆ = 1.4, phase locking regime. Blue : e = 1.7, ∆ = 2.6, bounded-phase regime. Violet : e = 1.7, ∆ = 2.9 ≃ fB. Red : e = 1.7, ∆ = 5.7, phase-drifting regime. Green : e = 0, no forcing.

Acknowledgments.

We thank F. Bondu for his help on the phase noise calculations. This work has been partly supported by the "Contrat de Projet Etat Région Bretagne" PONANT.


Bess Fang

Giuseppe Carleo

Isabelle Bouchoule

Quench-induced breathing mode of one-dimensional Bose gases

Keywords: numbers: 03, 75, Kk, 67, 85, -d

We measure the position-and momentum-space breathing dynamics of trapped one-dimensional Bose gases at finite temperature. The profile in real space reveals sinusoidal width oscillations whose frequency varies continuously through the quasicondensate to ideal Bose gas crossover. A comparison with theoretical models taking temperature into account is provided. In momentum space, we report the first observation of a frequency doubling in the quasicondensate regime, corresponding to a selfreflection mechanism due to the repulsive interactions. Its disappearance through the crossover is mapped out experimentally, giving insights to the dynamics of the breathing evolution.

The field of ultracold atomic and molecular gases has gained increasing significance in the study of nonequilibrium phenomena in quantum many-body systems. The accurate time-dependent control of microscopic parameters enables the realization of prototypical nonequilibrium processes, complementary to those studied in condensed matter physics. In particular, isolated quantum many-body systems can be driven out of equilibrium and accurately monitored. The realization of such experimental simulators [1] therefore constitutes a unique tool to investigate fundamental questions. Remarkable examples are the observation of a light-cone effect in the spreading of correlations [2,3], as well as the studies of quantum ergodicity [1,[4][5][6], where prethermalisation was observed in (nearly) integrable systems.

One-dimensional (1D) systems are ideal test benches for the study of out-of-equilibrium phenomena because of their intrinsic strong correlations and the possibility to realize integrable models. The experimental investigation of their quantum dynamics is of particular relevance in clarifying questions for which a theoretical treatment is challenging. The understanding of the time-dependent behavior of interacting quantum systems is indeed restricted to a few limiting cases [7], and plagued by the lack of systematic ab initio approaches. This stands in contrast to the well consolidated understanding of thermodynamic properties of such systems, for which exact analytical and numerical results are available [8].

Among the simplest out-of-equilibrium situation is the dynamics of a gas after a sudden change (quench) of the external harmonic potential in which the atoms are confined. The resulting coherent, breathing-like oscillations of the atomic cloud are governed by the collective excitations of the quantum gas, and serve to identify the quasicondensate (qBEC) regime [9] and strongly correlated phases in 1D [10]. The sum-rule approach has had a remarkable success in the prediction of the breathing frequencies at zero temperature [11,12]. Nonetheless, numerous fundamental questions remain open. These questions concern the general mechanism of mass transport in interacting quantum systems as well as its counterpart in momentum space. In particular, the effect of temperature on both the lifetime and the frequency of the collective mode cannot be accessed in terms of the standard sum-rule approach. Moreover, breathing oscillations of the momentum distribution provide complementary insight on the dynamics of interacting quantum gases, information which is not yet investigated experimentally nor accessible theoretically to date.

In this Letter, we study experimentally the breathing mode of harmonically confined 1D Bose gases at finite temperature, comparing the evolution in real and momentum space. We report the first observation of a self-reflection mechanism due to repulsive interaction, expected of a strongly interacting system [13]. We measure how the dynamics evolves through the qBEC to ideal Bose gas (IBG) crossover, and develop two complementary theoretical treatments for the breathing modes, taking finite-temperature effects into account. The combination of an exact short-time expansion of the dynamics and a long-wavelength hydrodynamic analysis reveals that the finite-temperature mass transport happens in an isentropic way, and provides a qualitatively good agreement with the experimental observations.

System and quench protocol. The Lieb-Liniger (LL) model of 1D bosons in continuum, with a Hamiltonian

H LL = j p 2 j 2m + j<k g 1D δ(z j -z k )
, where m is the particle mass, and g 1D (> 0) is the coupling constant, is realized experimentally for gases confined in elongated traps, provided the transverse degrees of freedom are frozen out [14]. Various experimental probes have been developed for such gases, permitting not only thermometry methods [15][16][17][18][19][20], but also the study of nonequilibrium behavior. Of special importance is the measurement of the momentum distribution, now feasible via Bragg spectroscopy [START_REF] Fabbri | Inguscio[END_REF] or focusing technique [20].

Our experiment prepares a single tube of 87 Rb gas using an atom-chip setup [17]. The final samples typically constitute 800 to 8000 atoms for this experiment. The in situ density profile [15] indicates a temperature around By varying τ , we map out the evolution in real and momentum space. Quasicondensate regime. We start by considering the qBEC regime: Figures 1(b) and(c) show ρ(z, τ ) and n(p, τ ) for a data set that lies in this regime (data A). The time evolution of these two quantities reveals a remarkable frequency mismatch. In order to elucidate its origin, we use the hydrodynamic equations (HDE) [36]

∂ τ ρ + ∂ z (ρv) = 0, (1)
∂ τ v + v∂ z v = -∂ z ( 1 2 ω 2 D z 2 ) -1 mρ ∂ z P, where ρ = ρ(z, τ ) and v = v(z, τ
) are the density and velocity fields, and P is the pressure. For a qBEC, P = g 1D ρ 2 /2, so that the scaling solution ρ(z, τ ) = [23] is valid provided ρ 0 is the steady state inverted parabola, and the scaling factor b obeys b +

1 b(τ ) ρ 0 z b(τ )
ω z (τ ) 2 b = ω 2 0 b 2 . Such ρ(z, τ ) is periodic in time with a frequency ω Bz ≃ √ 3ω D .
Neglecting thermal fluctuations, it follows that the momentum distribution satisfies [23] n(p, τ ) ∝ ρ 0 ( bp ḃm ), which is periodic in time with a frequency ω Bp = 2ω Bz . Minimal momentum width occurs when the real-space distribution is both the largest and the thinnest. The latter scenario, which does not occur for a non-interacting gas, corresponds to a selfreflection of the could due to repulsive interactions. The frequency doubling (ω Bp = 2ω Bz ) is thus a direct consequence of atomic interactions. In higher dimensions, the scaling solutions also predicts a frequency doubling of the breathing modes [37]. What may have prevented its observation is perhaps the fact that early experiments using time-of-flight (TOF) techniques do not measure the true momentum distribution due to the contribution of interaction energy.

It is worth mentioning that such a momentum-space frequency doubling in the oscillatory behavior of w p , the half width at half maximum (HWHM) in momentum space, is expected to occur for a Tonks gas (g 1D → +∞) [25] : the momentum distribution for large α oscillates between a Fermi-like and a Bose-like distribution, a thin Bose-like distribution occuring both when the cloud is the thinnest and the largest [13], while in the expansion and compression phases the momentum distribution is dominated by the large Fermi-like hydrodynamic component. This marks the difference between the breathing behaviour of a strongly interacting 1D Bose gas and that of a noninteracting gas, even though the breathing frequency in real space is ω Bz = 2ω D for both cases. Similar collective oscillations of strongly interacting 1D Bose gases have been studied in two experiments [4,10] to our knowledge, neither reporting any frequency doubling. In the former [4], a different excitation scheme is employed, potentially rendering the frequency doubling inaccessible. In the latter [10], the finite TOF may have prevented the direct measurement of the momentum distribution.

In our experiment the momentum-space frequency doubling is clearly seen in Fig. 1(b-c). Moreover, the data are in good agreement with the zero-temperature theory above, whose predicted real-and momentum-space density evolutions for a system initially at equilibrium in the trap of frequency ω 0 are shown in Fig. 1(d-e), with n(p, τ ) broadened to account for finite resolution [38]. For a more quantitative comparison, we show two frames of the instantaneous densities. At τ = 1 ms (minimal in situ width, Fig. 1(f)), the scaling solution predicts a vanishing momentum width, but thermal fluctuations dominate such that the finite-temperature momentum distribution computed using quantum Monte Carlo (QMC) methods [20] (dashed lines) accounts for the data better than the sole effect of resolution [START_REF]The division is heuristically justified by the increasingly significant wings for greater values of γ0[END_REF]. However, thermal fluctuations are small corrections at τ = 21 ms (maximal momentum width, Fig. 1 in excellent agreement with the inverted parabolas of the scaling solution. Self-similarity of the shapes is due to the fact the momentum distribution is dominated by the hydrodynamic velocity field, which is linear in position. The same phenomena is at the origine of the dynamic fermionization for a strongly interacting system [13].

Ideal Bose gas regime. For an IBG, a single-particle description suffices and breathing amounts to a rotation of the phase-space density, so that the widths in real and momentum space oscillate out of phase at the same frequency ω B = 2ω D . Fig. 2 shows a data set close to this regime (data C), where the in situ meansquare width z 2 (momentum HWHM w p ) is obtained from Gaussian (Lorentzian) fit. The antiphase is apparent from the plots. Fitting both time evolution with damped sinusoids, we measure identical breathing frequency, ω B /ω D = 1.84 ± 0.04, where ω D is determined by monitoring the center-of-mass (dipole) oscillations.

Crossover in real space. A central point to be understood is how the breathing frequency varies through the qBEC to IBG crossover. To address this question, we vary the total atom number in order to explore different regimes. The inset of Fig. 3 shows the region spanned by the data in the (γ, t) phase diagram of the LL model [28], where γ = mg1D 2 ρ is the (local) interaction parameter, and

t = 2 2 kB T mg 2 1D
is the reduced temperature. Each sample is characterized by γ 0 , evaluated at the peak density. We extract z 2 from fitting ρ(z) with either an inverted parabola for γ 0 < 0.004 [40] or a Gaussian otherwise. We obtain the real-space breathing frequency ω Bz by fitting z 2 (τ ) with a damped sinusoid. The measured ω Bz /ω D as a function of γ 0 , shown in Fig. 3(a), displays a smooth crossover between the asymptotic theoretical limits √ 3 and 2 [11]. Such a frequency shift was first observed in [9], where the authors probed the breathing behaviour of an ensemble of 1D gases in real space. In addition to the issue of thermalization mentioned in [9], the intrinsic en-semble averaging renders it problematic to characterize the system with a single value of temperature T or interaction parameter γ. In contrast, our current experiment with a single 1D system allows for a more quantitative study of the problem at hand.

In order to provide a theoretical treatment of the crossover, we devise two complementary approaches. On one hand, we model the crossover using the HDE Eq. ( 1). Since long-wavelength density waves in a fluid are adiabatic [START_REF] Eq | for a quench of finite amplitude (α = 1) a priori[END_REF], we use the isentropic pressure curves derived numerically from the YY EoS [START_REF] Yang | [END_REF]. The breathing mode frequency measured experimentally does not depend on the oscillation amplitude for the explored parameter range. We thus linearize the HDE for small displacement and ω Bz is obtained by solving an eigenvalue problem. Results evaluated at t = 1100 are shown as a dashed line in Fig. 3(a). On the other hand, we provide a microscopic treatment of the breathing frequency that accounts for the effect of temperature. Assuming the system at τ = 0 is at thermal equilibrium with a Hamiltonian H = H LL + H pot , H pot = mω 2 0 j z 2 j /2, the expansion of the Heisenberg equation of motion after

a quench ∆H = m(ω 2 D -ω 2 0 ) j z 2 j /2 gives ∆H (τ ) = ∆H T -τ 2 2 C 2 T + τ 4 4! C 4 T + .
. ., where C 2 and C 4 are the second and fourth order nested commutators with H f = H + ∆H, and the thermal average T is taken over the thermal state of H. Suppose the time evolution ∆H (τ ) ∝ z 2 (τ ) is purely sinusoidal at the frequency ω Bz , we have

ω 2 Bz ω 2 D = C 4 T C 2 T = 4 - 1 2
H int T H pot T , (2)
where H int = g 1D j<k δ(z j -z k ) is the interaction part of the Hamiltonian. The final equality comes from explicitly working out C 2 T and C 4 T , and taking the limit of an infinitesimal quench amplitude α → 1, where the sinusoidal approximation is valid [42]. The shaded region shows the crossover criteria γco = t -2/3 (the dashed line in the inset), at t ≃ 1100 given by the in situ profile with 20% uncertainly. Inset: the phase diagram of the LL model [28], where all lines represent smooth crossovers. The horizontal line segment shows the region explored by the data [20].

It is remarkable that both approaches (HDE and the exact short-time expansion) give the same results for the breathing frequencies through the crossover. This indicates that the finite-temperature mass transport is indeed isentropic. Moreover, Eq. ( 2) extends the regime of validity of the existing sum-rule approaches, which erroneously predicts ω Bz = √ 2ω D in IBG, whereas the true asymptotic limit is ω Bz = 2ω D .

Our theories provide a qualitatively good agreement with the experimental data, despite an overall overestimation of the breathing frequency. This quantitative departure most likely indicates that the atomic ensemble is not at thermal equilibrium at the end of the preparation (evaporation). Although the in situ profiles calibrated with the YY EoS generally yields about 100 nK, the atom-number fluctuation indicates approximately 35 to 45 nK at the center of the cloud. We thus include in Fig. 3(a) our prediction at t = 400 (dash-dotted line). We believe that such a lack of equilibrium is a direct consequence of evaporating into the integrable (1D) regime. Further investigation in this direction is underway.

Disappearance of the self-reflection mechanism. To address the behavior of the breathing mode in momentum space, we investigate the time evolution of w p , ex-tracted from a Lorentzian fit of n(p). Since w p shows a periodic behavior at the frequency ω Bz , it can be expanded in a discrete Fourier spectrum. The relative weight of the Fourier components varies through the crossover, and we obtain quantitative information by fitting w p with a function

y = Ae -τ τ 1 +Be -τ τ 2 √ K cos(ω Bz τ )- √ 1 -K cos(2ω Bz τ ) , (3)
shown in Fig. 2 (bottom, solid line), with A fixed at the average width during the initial cycle and the phase corresponds to a minimal in situ width at the start of the oscillations. The relative power of the first harmonic K (squares) and second harmonic 1 -K (diamonds) are shown in Fig. 3(b) as a function of γ 0 . In qBEC, the second harmonic dominates, as predicted by the scaling solution, and signals the self-reflection mechanism. In IBG, the first harmonic dominates, as expected for a noninteracting gas where the self-reflection is absent. Both weights vary gradually through the crossover, indicating a smooth disappearance of the self-reflection mechanism. This can be seen as the effect the breathing mode has on the thermally excited Bogoliubov modes of high energy, e.g. their frequency and wavefunction would be modulated in time, such that w p is larger at minimal z 2 than that at maximal z 2 , see Fig. 2 (data B). The periodicity at 2ω Bz is then broken and the first harmonic at ω Bz emerges. Figure 3 shows that the first harmonic starts to gain weight at a value of γ 0 significantly smaller than that where the frequency shift takes place in real space.

The breathing mode observed has lifetime estimated to be on the order of seconds. This is in stark contrast with situations in 3D where damping of the collective modes occurs, mainly via Landau damping mechanism [START_REF] Jin | [END_REF]34]. The long lifetime of the breathing mode in 1D may be related to the integrability of the underlying LL model.

Conclusions. In this Letter we have probed the breathing mode in real and momentum space through the qBEC to IBG crossover. The shift of the real-space frequency between the asymptotic values √ 3ω D and 2ω D is demonstrated. Our theory models that assume thermal equilibrium before the quench do not agree with the measurement quantitatively, indicating the possibility of a non-Gibbs initial state produced by evaporation. We report the first observation of a momentumspace frequency doubling in qBEC, corresponding to an interaction-induced self-reflection mechanism that is expected of a Tonks gas. We experimentally map out the disappearance of the self-reflection through the crossover, for which no theoretical predictions exist to date. This illustrates the richness of out-of-equilibrium dynamics and the importance of the momentum-space observation to unveil the underlying physics.

FIG. 1 .

 1 FIG. 1. (Color online) Quench sequence and subsequent density evolution in position and momentum space for a qBEC (data A). (a) Longitudinal trapping frequency ωz as a function of time τ . (b) ρ(z, τ ), (c) n(p, τ ), experimental data (in atoms/pixel). (d), (e), corresponding plots from ab initio scaling calculation. (f) and (g) show instantaneous ρ(z), n(p) at τ = 1 and 21 ms. Experimental data (dots) are compared with the scaling solutions (solid lines). QMC calculation of n(p, τ ) (dashed line) is shown in (f). z (p) is in units of pixels ∆z (∆p), with ∆z = 2.7 µm, and ∆p = 0.14 /µm.

FIG. 2 .

 2 FIG. 2. Time evolution of the in situ mean-square width z 2 (top) and the momentum HWHM wp (bottom) for 3 data sets. The solid lines show the fit: a damped sinusoid in real space (top), and a two-harmonic model according to Eq. (3) in momentum space (bottom). Statistical error of the widths are shown for the first 100 ms. Data A-C differ in thermodynamic regime. See Fig. 3 and text for details.

FIG. 3 .

 3 FIG.3. Breathing mode through the crossover. (a) breathing frequency in units of dipole frequency, as a function of γ0.Our theory evaluated at t = 1100 (dashed line) and t = 400 (dash-dotted line) are shown. (b) relative power of the first (square) and second (diamond) harmonic of wp. Data labels (A-C) correspond to those in Fig.2. The errorbars account for fitting error only. The dotted lines are the asymptotic limits. The shaded region shows the crossover criteria γco = t -2/3 (the dashed line in the inset), at t ≃ 1100 given by the in situ profile with 20% uncertainly. Inset: the phase diagram of the LL model[28], where all lines represent smooth crossovers. The horizontal line segment shows the region explored by the data[20].

The authors would like to thank A. Minguzzi, T. Roscilde, S. Stringari, P. Vignolo, M. Zvonarev for stimulating discussions. This work is financially supported by Cnano IdF, the Austro-French FWR-ANR Project I607, and the FP7-Marie Curie IEF grant 327143.