Matthieu Fruchard

Controllability and Control Synthesis of Underactuated Magnetic Microrobots

Keywords: controllability, underactuated multi-agent system, nonlinear control synthesis, medical microrobotics

The control of multiple magnetic microrobots is of particular interest for therapeutic applications. Yet the controllability of such a system is not straightforward since, on most of such systems designs, there is only a single control input per axis. The paper addresses the controllability and control synthesis for several magnetic microrobots navigating in blood vessels. First, controllability requires the microrobots to magnetically interact in order to achieve trajectory tracking along an admissible reference trajectory, whether swimming at low or high Reynolds. Then the resulting nonlinear system is shown to be diffeomorphic to different nonlinear canonical forms depending on the choice of the microrobots, so that the stabilization of the system along any admissible reference trajectory can be achieved using a backstepping controller synthesis, yet sometimes at the price of a zero dynamics stability analysis. Simulation results illustrate the efficiency of the proposed approach.

Introduction

The control of untethered microrobots using the vascular network is receiving a growing interest since they can perform minimally invasive and targeted surgery or diagnosis, whilst accessing to remote places with lessened side effects. Deported magnetic actuation is appealing since it avoids an embedded energy source, resulting in an improved miniaturization and payload ratio. Different propulsion designs have been studied: elastic flagellated [START_REF] Lagomarsino | A simulation study of the dynamics of a driven filament in an aristotelian fluid[END_REF]; [START_REF] Evans | Propulsion by passive filaments and active flagella near boundaries[END_REF], helical tailed [START_REF] Dreyfus | Microscopic artificial swimmers[END_REF]; [START_REF] Zhang | Artificial bacterial flagella for micromanipulation[END_REF], or bead pulled robots [START_REF] Abbott | How should microrobots swim[END_REF]; [START_REF] Mathieu | Method of propulsion of a ferromagnetic core in the cardiovascular system through magnetic gradients generated by an MRI system[END_REF]. This emerging field has at first mainly focused on feasibility [START_REF] Nelson | Microrobots for minimally invasive medecine[END_REF], then on controllability [START_REF] Giraldi | Local controllability of the two-link magneto-elastic micro-swimmer[END_REF] and control and observation [START_REF] Fruchard | Estimation of the blood velocity for nanorobotics[END_REF] issues for a single robot. However some applications, such as drug targeting or brachytherapy, require that enough payload is released in the vicinity of a tumor. On the other hand, avoiding arterial embolization hazards induces limitations on the microrobot size, and so on the available payload, hence the interest in controlling multiple robots.

Literature abounds with micro multi-agent strategies, Email address: matthieu.fruchard@univ-orleans.fr (Matthieu Fruchard). 1 Corresponding author. yet there is an obstruction in the present case: the microrobots are often controlled using a single control input, either using electrostatic setups [START_REF] Donald | Planar microassembly by parallel actuation of MEMS microrobots[END_REF], electromagnetic setups [START_REF] Floyd | Two-dimensional contact and noncontact micromanipulation in liquid using an untethered mobile magnetic microrobot[END_REF]; [START_REF] Frutiger | Small, fast, and under control: wireless resonant magnetic micro-agents[END_REF]; [START_REF] Servant | Controlled in vivo swimming of a swarm of bacteria-like microrobotic flagella[END_REF], or Magnetic Resonance Imaging (MRI) devices [START_REF] Martel | Using a swarm of selfpropelled natural microrobots in the form of flagellated bacteria to perform complex micro-assembly tasks[END_REF] actuators. Consequently, multiple magnetic robots modeling results in a nonlinear underactuated system. Vartholomeos et al. (2012) parameterized the input as a sequence of pulses with different widths that exploited the induced motion difference to control independently different robots. However, this work aimed at a position open-loop control and therefore suffered from a lack of robustness. [START_REF] Floyd | Control methodologies for a heterogeneous group of untethered magnetic micro-robots[END_REF]; [START_REF] Diller | Control of multiple heterogeneous magnetic microrobots in two dimensions on nonspecialized surfaces[END_REF] proposed to exploit the natural frequency differences between some robots to actuate one or part of the considered magnetic microrobots with a point to point control objective. In [START_REF] Salehizadeh | Two-agent formation control of magnetic microrobots in two dimensions[END_REF] the authors took advantage of a rotating magnetic field to control the distance and orientation of two agents operating in close vicinity through the modulation of their attractive-repulsive interaction. [START_REF] Eqtami | MRI-powered closed-loop control for multiple magnetic capsules[END_REF] addressed more formally the control of two microrobots using model predictive control to perform a fixed point stabilization. To exhibit a solution, they assumed the existence of a solution to the initial optimal control objective, i.e. that the system is locally controllable from the initial point, but did not focused on controllability issues. The controllability of two non inertial micro-robots has already been studied in Salehizadeh andDiller (2017, 2021), yet only considering their relative position controllability; a full state controllability analysis has been provided in [START_REF] Fruchard | Local controllability, trajectory planning, and stabilization of a two-agent underactuated microrobotic system[END_REF] using the linear test, yet only a sufficient condition is given.

The present paper addresses the control of magnetic interacting microrobots facing the pulsatile blood flow in order to stabilize their full state along a reference trajectory, using the same single control input. The paper contribution is twofold: i) conditions under which these systems are controllable are investigated, which is -to our knowledge-the first full state commandability result for such an underactuated multi-agent microrobotic system, and ii) Lyapunov stabilizing control laws are synthesized exploiting transformations to nonlinear canonical forms, depending on the system design and parameters. Finally, simulations results illustrate the stability of the proposed controller, for an inertial 2-microrobot system and for a non inertial 5-microrobot system.

B M v f (t) Fm k+1 Fm k F k+1 p k-1 p k r k-1,k F k Microrobot Rk Microrobot Rk+1 F dk Fm k-1 F dk-1 Microrobot Rk-1 Pulsatile Blood Flow velocity vf(t) p k+1 ı 0 F k-1 u = ∇Bx ∈ IR F dk+1 r k,k+1
Fig. 1. Forces exerting on a magnetic microrobot (R k ) and its closest neighbours (R k-1 ), (R k+1 ): the drag force F d k , the magnetic motive force Fm k , and F k the net magnetic force interaction acting on robot (R k ). A single control input u = ∇Bx ∈ R is used to control all the microrobots in a frame F (0, i).

Modeling

Spherical microrobots (R k ) navigating in a blood vessel are localized by their positions p k in the frame F(0, i), see Figure 1. Each microrobot is affected by the drag force, the motive magnetic force and the magnetic interaction between the robots, detailed hereinafter. The magnetic field is supposed stationary, like in MRI devices.

Forces

In the blood, a moving spherical microrobot (R k ) of radius r k experiences a drag force that opposes its motion. Even if the proposed approach can be extended mutatis mutandis to nonlinear drag force or to non-Newtonian blood modeling, as e.g. in [START_REF] Arcese | Endovascular magnetically-guided robots: navigation modeling and optimization[END_REF][START_REF] Arcese | Adaptive controller and observer for a magnetic microrobot[END_REF], it is here assumed -for sake of simplicity-that the drag force is linear:

F d k = -m k d k ( ṗk -v f (t)) ı, d k = 9η 2β k ρ k r 2 k (1)
where β k is a dimensionless ratio related to the partial vessel occlusion by the robot [START_REF] Haberman | Motion of Rigid and Fluid Spheres in Stationary and Moving Liquids Inside Cylindrical Tubes[END_REF], and ṗk -v f (t) denotes the relative velocity between the robot (R k ) and the fluid. η, ρ k and m k denote respectively the blood viscosity, the microrobot (R k ) density and mass.

A robot of magnetization M k and ferromagnetic ratio τ m k exposed to a magnetic field gradient input u = ∇B x is pulled by a magnetic force [START_REF] Arcese | Endovascular magnetically-guided robots: navigation modeling and optimization[END_REF]:

F m k = m k a k u ı, a k = τ m k M k /ρ k .
(2)

Let µ 0 and M k = 4 3 πτ m k r 3 k M k denote the vacuum permeability and the (R k ) microrobot magnetic moment, and r k,i = (p i -p k ) ı. The magnetic interaction force F k exerted by other microrobots on (R k ) is given by [START_REF] Vartholomeos | In silico studies of magnetic microparticle aggregations in fluid environments for MRI-guided drug delivery[END_REF]. Since B = B 0 ı, and since this force is short ranged, we consider only a twoclosest neighbours interaction, so:

F k = m k j k-1 r k-1,k 4 - i k r k,k+1 4 ı i k = 3µ0M k M k+1 2πm k r k+1,k • ı r k,k+1 , j k = i k m k m k+1 , i k j k ≥ 0.
(3)

Two Inertial Microrobots

Let x T = p 1 p 2 ṗ1 ṗ2 ∈ X denote the state vector where X is any connected subset of either

X + = {x ∈ R 4 : x 1 -x 2 > 0} or X -= {x ∈ R 4 : x 1 -x 2 < 0}
depending on the relative initial positions x 1 and x 2 .

Applying the Newton's second law using (1), ( 2) and (3) yields:

(S x )              ẋ1 = x 3 ẋ2 = x 4 ẋ3 = -d 1 (x 3 -v f ) - i1 (x1-x2) 4 + a 1 u ẋ4 = -d 2 (x 4 -v f ) + j1 (x1-x2) 4 + a 2 u (4)
Along any C 1 admissible reference trajectory (x r , u r ) with x T r (t) = x 1r x 2r ẋ1r ẋ2r (t) satisfying (4), there exists a

C 2 distance function d 1,2 (t) ∈ D ⊂ R \ {0}: d 1,2 (t) = x 1r (t) -x 2r (t). Let x = x -x r
denote an error state, the resulting error dynamics is an affine nonlinear control system with drift:

(S i 2 )                    ẋ = f (x) + g ũ f (x) =        x3 x4 -d 1 x3 -i 1 δ 1 (x) -d 2 x4 + j 1 δ 1 (x)        , g =        0 0 a 1 a 2        (5a) (5b) with ũ = u-u r (t), X a connected subset of either X + = {x ∈ R 4 : x1 -x2 + d 1,2 (t) > 0} or X -= {x ∈ R 4 : x1 -x2 + d 1,2 (t) < 0}
, and the nonlinear function:

δ 1 : x → 1 (x 1 -x2 + d 1,2 (t)) 4 - 1 d 4 1,2 (t)
.

(6)

A Finite Number m of Non Inertial Microrobots

In a Stokes flow, R e ≪ 1 so inertial effects can be neglected [START_REF] Purcell | Life at low reynolds number[END_REF]; [START_REF] Tikhonov | Systems of differential equations containing small parameters in the derivatives[END_REF]. Rewriting (5) for m non inertial microrobots leads to the following error system:

(S ni m )                            ẋ = f (x) + g ũ f (x) =            -ǫ 1 δ 1 . . . ε k-1 δ k-1 -ǫ k δ k . . . ε m-1 δ m-1            , g =            α 1 . . . α k . . . α m            (7a) (7b)
where x = p 1 -x 1r . . . p k -x kr . . . p m -x mr T ∈ Xm ⊂ R m is the state vector composed of the position errors of the m robots (R k ), and:

δ k : x → 1 (x k -xk+1 + d k,k+1 (t)) 4 - 1 d 4 k,k+1 (t) (8) with d k,k+1 (t) the (R k ) to (R k+1 ) reference distance, ǫ k = i k d k , ε k = j k d k+1 so that ǫ k ε k ≥ 0 and α k = a k d k . Sys- tem S ni m is defined on any connected subset Xm of X - m = {x ∈ R m : x1 < x2 -d 1,2 < . . . < xm - m-1 i=1 d i,i+1 } or X + m = {x ∈ R m : x1 > x2 -d 1,2 > . . . > xm - m-1 i=1 d i,i+1 }.
Since the magnetic interaction is a short range force, one considers that microrobot (R k ) interacts only with its two closest neighbours (R k-1 ) and (R k+1 ).

Problem statement

In system S x , the positions p = (p 1 p 2 ) and velocities v = ( ṗ1 ṗ2 ) belong to a 2D space, whereas the control input u is scalar; underactuation degree is even higher for system S ni m . Thence we face an underactuated microrobotic system with drift whose full state controllability has not yet been addressed in the literature, even though partial controllability has been addressed in Salehizadeh andDiller (2017, 2021). Another consequence of the underactuation is that all reference trajectories are not feasible by the microrobotic system. Admissible reference trajectories thus have to be defined to solve the trajectory planning issue, e. g. using differential flatness whenever the system is flat or the procedure detailed in [START_REF] Fruchard | Local controllability, trajectory planning, and stabilization of a two-agent underactuated microrobotic system[END_REF]. Finally, how to synthesize a control law that stabilizes the state x along any admissible reference trajectory? One can rely either on linear or on nonlinear approaches. Yet, even with a single microrobot, works based on linear controllers, e.g. [START_REF] Tamaz | Real-time MRI-based control of a ferromagnetic core for endovascular navigation[END_REF]; [START_REF] Mehrtash | Design and implementation of LQG/LTR controller for a magnetic telemanipulation system-performance evaluation and energy saving[END_REF]; [START_REF] Choi | Position stabilization of microrobot using pressure signal in pulsating flow of blood vessel[END_REF], report instabilities and important oscillations for PID or LQR controllers can not deal properly with nonlinear perturbated systems. Besides, in the case of a multi-agent system, the nonlinear terms δ k (x) in ( 5)-( 7) depend on the inverse fourth power of the (R k ) to (R k+1 ) distances. Linearizing these terms ( 6)-( 8) around a given trajectory will lead to severe errors even for small position errors, so this additional nonlinearity will further penalize robustness and stability weaknesses of control laws based on the linearized system. Amongst nonlinear approaches, differential flatness is widely used, yet S x is not generically linearizable by a dynamic feedback, so this approach can not apply to the generic case, as it is proven in section 4.4. Another idea is to apply a state feedback with a time-varying gain K(t) obtained by solving a differential Riccati equation; however solving the latter analytically is hardly possible. As in [START_REF] Olfati-Saber | Nonlinear Control of Underactuated Mechanical Systems with Application to Robotics and Aerospace Vehicles[END_REF], a promising approach is to find how to map the system to normal feedback forms in order to propose backstepping controllers.

Main Results

The twofold contribution addresses first a controllability analysis for both the inertial microrobotic system S i 2 and the non inertial system S ni m , and then a control synthesis to Lyapunov stabilize these systems along any admissible reference trajectory.

System controllability

Proposition 1 System S i 2 is strongly accessible and small-time controllable around x = 0 if and only if both following conditions are satisfied: C1 Magnetic interaction exists: (i 1 , j 1 ) = (0, 0).

C2 The microrobots parameters

(a 1 , a 2 , d 1 , d 2 , i 1 , j 1 )
and the planned distance d 1,2 (t) satisfy:

4(a 2 -a 1 ) 2 (a 1 j 1 +a 2 i 1 ) a 1 a 2 d 5 1,2 (t) = (d 2 -d 1 )(d 1 a 2 -d 2 a 1 ). (9)
Remark 1 Condition C1 clearly requires -contrary to most previous studies-that the microrobots operate in close vicinity. Otherwise, there is no magnetic interaction, i.e. i 1 = j 1 = 0, leading to a noncontrollable linear system. In such a case, the system is simply stabilizable on a line depending on the initial condition x(0). More precisely, both velocities are controllable as well as a linear combination of the microrobots positions: the system position is stabilized on a line λx 1 + µx 2 = 0 instead of the origin x = 0, with λ, µ depending on x(0), i.e. even point stabilization can not be achieved.

Remark 2

The microrobots must have different designs as a consequence of condition C2 since (9) is not satisfied if the microrobots are identical (d 1 = d 2 and a 1 = a 2 ), whatever d 1,2 (t). The resulting system is then not strongly accessible and consequently not controllable. If the two robots share the same magnetic sensitivity a, yet have different drag coefficients d i , condition C2 is fulfilled and the system is controllable for all

d 1,2 (t) ∈ D ⊂ R \ {0}.
When swimming at very low Reynolds Re ≪ 1, the controllability result for an m non inertial microrobots underactuated system is given in the following proposition.

Proposition 2 Consider (m -1) identical microrobots (R k ), so ǫ k = ε k = ǫ and α k = α, ∀k ≤ m -1. Then system S ni m is small-time controllable around x = 0 if and only if ǫ i , ε i = 0 and the last m-th robot (R m ) satisfies α m = α.
Remark 3 In the case of m = 3 non inertial microrobots with no common parameter, the strong accessibility Lie rank condition dim ∆ C0 (x) = 3 is fulfilled if and only if

̟δ (1) 1 δ (1) 2 δ (1) 1 γ 1 (γ 1 ε 1 + γ 2 (ε 1 + ǫ 1 )) -δ (1) 2 γ 2 (γ 2 ǫ 2 + γ 1 (ε 2 + ǫ 2 )) = 0 (10) where γ k = α k -α k+1 , ̟ = ǫ 1 ε 2 α 2 + ε 1 ǫ 2 α 1 + ǫ 1 ǫ 2 α 3 and δ (1) k (x) = ∂δ k ∂ xk (x).
Considering one non magnetic microrobot (ǫ i = ε i = 0 for a given i) hinders the system controllability so for an m = 3 non inertial microrobotic system, C1 is unchanged, while C2 is replaced by (10). Of course, it is possible to satisfy (10) thanks to an appropriate choice of the distances d k ; nonetheless, in the generic case considering m different microrobots, condition (10) last factor is a degree (m-1) polynomial in δ

(1) k (x) whose solving is not straightforward. That is the reason why Proposition 2 focuses on a simpler design consistent with standardized microrobots contruction processes.

Remark 4 If the drag force is nonlinear e. g. because of non-Newtonian effects, the proposed approach can still be used considering a first order Taylor expansion approximation of the drag force, i. e. linearization along the reference trajectory. Since the partially linearized system is controllable, then the nonlinear system is locally controllable.

Proof of Proposition 1

If C1 is not satisfied, nonlinear terms are removed from (5) so ẋ = Ax + B ũ, and controllability is checked using the Kalman rank condition. Let Γ = (B AB . . . A 3 B). Then rank(Γ) = 3 < 4, so system S i 2 is not fully controllable: only three modes are controllable and one is not. In closed-loop, the system poles can be placed at 0, λ 1 , λ 2 , λ 3 , so the system is simply stable and ẋ = A cl x, with e.g. distinct eigenvalues. A cl is thus diagonalizable: there exists T = (t ij ) such that the diagonal matrix

D = T -1 A cl T first entry is null. Since A cl first two lines are (0 2 I 2 ) , then t 31 = t 41 = 0. Besides, since T is in- vertible, it follows that t 2 11 +t 2 21 = 0. Let z = T -1 x, then ż = Dz. As a result, it is clear that ∀t ≥ 0 z 1 (t) = z 1 (0) and lim t→∞ (z 2 z 3 z 4 )(t) = (0 0 0). Since x = T z, it is straightforward that lim t→∞ x(t) = (t 11 t 21 0 0) z 1 (0).
The strong accessibility algebra C 0 is the smallest algebra that contains the control vector fields and such that [f, C 0 ] ⊂ C 0 , with [• , • ] denoting the Lie bracket operator. Let ∆ C0 (x) = span{X(x) : X ∈ C 0 }, the strong accessibility Lie rank condition dim ∆ C0 (x) = n = 4 is a controllability necessary condition Nijmeijer and van der [START_REF] Nijmeijer | Nonlinear Dynamical Control Systems[END_REF]. Let ∆ a = a 1 -a 2 and ∆ d = d 1 -d 2 , computing the successive brackets using the adjoint notation yields

ad 1 f g =        a 1 a 2 -d 1 a 1 -d 2 a 2        , ad 2 f g =        -d 1 a 1 -d 2 a 2 d 2 1 a 1 + i 1 δ (1) 1 ∆ a d 2 2 a 2 -j 1 δ (1) 1 ∆ a        (11) ad 3 f g =        d 2 1 a 1 + i 1 δ (1) 1 ∆ a d 2 2 a 2 -j 1 δ (1) 1 ∆ a i 1 δ (1) 1 (a 2 (d 1 + d 2 ) -2d 1 a 1 ) -d 3 1 a 1 j 1 δ (1) 1 (a 1 (d 1 + d 2 ) -2d 2 a 2 ) -d 3 2 a 2        ad 1 f g, ad 2 f g = 0 0 i 1 -j 1 T δ (2) 1 ∆ 2 a ad 1 f g, ad 3 f g =        -i 1 ∆ a j 1 ∆ a i 1 (2d 1 a 1 -(d 1 + d 2 )a 2 ) j 1 (2d 2 a 2 -(d 1 + d 2 )a 1 )        δ (2) 1 ∆ a (12) with δ (i) 1 = ∂ i δ1 ∂ xi 1
. First, using (5b) and successive brackets ( 11)-( 12) and higher length brakets, span{g, ad 1 f g, ad 2 f g , ad 1 f g, ad 1 f g, ad 3 f g , . . .} is not full rank only for ∆ a = 0. Besides, span{ad i f g} i∈N has a rank loss only if condition C2 is violated. So dim L(f, g) < 4 when both distributions are singular, i. e. when ∆ a = ∆ d = 0, and the system is thus not accessible at 0, and consequently not controllable when the microrobots are identical. Now, when C2 is fulfilled, since dim span{ad i f g} i∈N = 4, then the system is strongly accessible at 0. However, strong accessibility does not imply controllability for a drift system since bad brackets can be an obstruction to controllability if they are not neutralized (spanned) by good brackets of smaller length. Bad brackets are iterated Lie brackets containing an odd number n f of f and an even number n g of g, like f, ad 2 g f, [ad 1 f g, ad 2

f g] and so on. The length of an iterated Lie bracket X is computed as L X = θn f + n g for some θ ∈ [0, 1]. In the present case, f (0) = 0 and ad 2 g f = 0, so these bad brackets are neutralized. The bad bracket [ad 1 f g, ad 2 f g] of length L b = 3θ + 2 is neutralized by lower good brackets (the ad i f g of length L g = iθ+1 < L b ) as long as condition C2 is fulfilled since in such a case dim span{ad i f g} i∈N = 4. Longer bad brackets are also neutralized since these lower length good brackets ad i f g already span the full space.

The only case still to be clarified is controllability when C2 is violated. In [START_REF] Kawski | A necessary condition for local controllability[END_REF] Kawski gives a new necessary condition for controllability: if [ad 1 f g, ad 2 f g] is not spanned by g, ad 1 f g, . . . , [ad k f , ad 3 g f ] for k ∈ N, then the system is not small-time local controllable. It is not difficult to show using (5b) that ad 2 g f = 0 so [ad k f , ad 3 g f ] = 0 for all k. Yet dim span{ad i f g} i∈N = 3 when condition C2 is violated, so [ad 1 f g, ad 2 f g] / ∈ span{ad i f g} i∈N . Consequently the system is small-time local controllable if and only if both conditions C1 and C2 are fulfilled.

Proof of Proposition 2

First note that when α m = α, ad i f g = 0 for i ≥ 1: the only non null brackets in the Lie algebra L(f, g) are f and g so its dimension is only 2 < m and hence the system is not controllable from any point. Lengthy yet simple induction shows that computing successively the Lie brackets h i (x) = ad i f g(x) for i ≤ m -1 using the vector fields (7), the following properties hold for i ≥ 1:

P1: h i j (x) = 0, ∀j ≤ m -i -1; P2: h i m (x) = -εm-1 ǫm-1 m-1 j=1 h i j (x)
for any x provided that m ≤ 3 and at x = 0 for any m;

P3: h i m-i (x) = (-1) i+1 ǫ m-1 ǫ i-1 (α -α m ) i k=1 δ (1) m-k (x).
The dimension of the resulting strong accessibility algebra is thus m as long as det C = det ad 0 f g, . . . , ad m-1 f g = 0. Using properties P1-P2 and adding to the last row of C the sum of all other rows weighted by εm-1 ǫm-1 yields:

det C = α 0 . . . h m-1 1 . . . . . . ⋆ α 0 h 2 m-2 ⋆ α h 1 m-1 ⋆ ⋆ s 0 . . . 0 (13)
with s = α m + (m -1)α εm-1 ǫm-1 . Using the Laplace expansion along the last row, it follows that det C = (-1)

m+1+ (m-1)(m+2) 2 s m-1 i=1 h i m-i so property P3 yields det C = ǫ m-1 α m + ε m-1 (m -1)α (α m -α) m-1 ǫ (m-1)(m+2) 2 ǫ m-2 m-1 m-1 i=1 δ (1) i i . (14)
Since ǫ m-1 ε m-1 ≥ 0 and the δ

(1) i (x), i ≤ m -1, do not cancel on Xm , controllability thus requires a non null magnetic interaction, and also requires that α m = α.

Lyapunov stabilizing controller

The proposed approach consists in two parts: first diffeomorphims mapping systems S i 2 and S ni m to normal forms S i are exhibited depending on the parametric values, so the control synthesis for these normal forms then provides a Lyapunov stabilizing controller for the original systems S i 2 and S ni m full state in any case. In the sequel, A denotes a matrix in prime form, F and G are vector fields whose only last component is not null, and L i f h denotes the i-th Lie derivative of h along f . Lemma 1 Let cases A, B and C denote the parametric sets (a

1 = a 2 , d 1 = d 2 ), (a 1 = a 2 , d 1 = d 2 ) and (a 1 = a 2 and d 1 = d 2 ) of S i
2 , respectively. Let j ∈ {A, B, C}, i = 4 for cases A and B, i = 3 for case C and

           h A = h C : x → a 2 x1 -a 1 x2 h B : x → d 1 x1 -d 2 x2 + x3 -x4 z k = φ j k (x) = L k-1 f h j (x), k ≤ i ξ = φ C 4 (x) = x1 -x2 (15a) (15b) (15c) (15d)
On X , under controllability conditions of Proposition 1, the diffeomorphims φ A and φ B map (5) to a pure feedback form S 1 : ż = Az + F (z) + G(z)ũ in cases A and B, with i = 4. In the generic case C, the diffeomorphism φ C defined for i = 3 maps (5) to

(S 2 ) ż = Az + F (z, ξ) + Gũ ξ = f ξ (z, ξ) (16)
Besides, S 2 is a minimum phase system with respect to the origin (z, ξ) = 0 provided that

a 1 -a 2 d 2 -d 1 (a 2 i 1 + a 1 j 1 ) > 0. ( 17
)
The Lemma proof is given in Appendix A.

Remark 5 In the particular cases A and B of Lemma 1, the system is differentially flat with a flat output y = h(x) in (15a) and (15b). However, the system is not generically flat since it requires that the distributions ∆ 0 = span{g} and ∆ i+1 (x) = span{∆ i , [f, ∆ i ]}(x) are involutive, i.e. stable by Lie bracketing. Using the vector fields given by (5b) and their brackets (11)-( 12), it is not difficult to show that ∆ 2 (x) is involutive only in cases A or B. It is possible to propose a control law synthesis -and even a trajectory planning-using differential flatness in either case, yet we would rather present a control synthesis encompassing all possible cases, i.e. A, B or C.

Lemma 2 Provided that the controllability condition of Proposition 2 is satisfied, let δ

(i) k (x) = ∂ i δ k ∂ xi k (x) and        h(x) = -[(m -2)ε m-1 α + ǫ m-1 α m ]x 1 +αε m-1 m-1 i=2 xi + αǫ m-1 xm z i+1 = L i f h(x), i = 0, . . . , m -1, (18)
then system S ni m is mapped to a pure feedback form S 1 : ż = Az + F (z) + G(z)ũ with A in prime form and F, G whose only non null entry is the last one with

G m (z) = ǫ m-2 ǫ m-1 [(m -1)ε m-1 α + ǫ m-1 α m ] (α -α m ) m-1 i=1 δ (1) i = 0, ∀x ∈ Xm . ( 19
)
The proof of this Lemma is deferred to Appendix B.

Proposition 3 Let x r (t), u r (t) denote any admissible C 1 trajectory and the associated reference control input for S i 2 or S ni m . Provided that assumptions of Propositions 1 and 2 are fulfilled for systems S i 2 and S ni m , respectively, the control law u = κ i (z, ξ) given by

κ i (z, ξ) = u r (t) -G -1 i (z) F i (z, ξ) + i j=1 c j i+1 z j (20)
ensures the semiglobal exponential stability of the origin for any initial bounded state x(0) ∈ X0 where c j i are constants recursively defined for gains k j > 0 by:

c i i = 1, c 1 i+1 = k i c 1 i + c 1 i-1 , c i i+1 = i j=1 k j c j i+1 = k i c j i + c j-1 i + c j i-1 , j ∈ [2, i -1] (21)
with a connected set X0 ⊂ X and i = 4 for system S i 2 in cases A and B, i = 3 for system S i 2 in case C, and X0 ⊂ Xm and i = m for system S ni m , respectively.

Proof of Proposition 3

For system S i 2 , the backstepping procedure is iterated from step 1 to step i = 4 on system S 1 in cases A and B, and up to i = 3 and using the minimum phase property proven in Lemma 1 for system S 2 in case C. For system S ni m , the procedure is iterated up to step i = m.

Differentiating a first candidate Lyapunov function

V 1 = 1 2 z 2 1 leads to V1 = z 1 z 2 . Set z1 = z 1 and z2 = z 2 + k 1 z 1 = 2 r=1 c r 2 z r with k 1 > 0 results in V1 = -k 1 z2 1 + z1 z2 , with c 1 2 = k 1 , c 2 2 = 1 from (21). The last term is counterbalanced at next step. 2. Let V 2 = V 1 + 1 2 z2
2 denote a second candidate Lyapunov function. From z2 definition and using the normal form of systems ż = Az + F (z) + G(z)ũ, it follows that ż2 = z 3 + k 1 z 2 so V2 = -k 1 z2 1 + z2 (z 1 + k 1 z 2 + z 3 ). Let k 2 > 0 and set z3 = 3 r=1 c r 3 z r with c r 3 given by ( 21), then we obtain

V2 = -k 1 z2 1 -k 2 z2 2 + z2 z3 (22) j. At step j ≤ i -1, define V j = V j-1 + 1 2 z2 j , where zj = j r=1 c r j z r (23)
with c r j given by ( 21). From steps 1 to j -1 we have

Vj-1 = - j-1 r=1 k r z2 r + zj-1 zj . (24)
Since, using (21), ς = zj-1 + żj + k j zj simplifies as

ς = (c 1 j-1 + k j c 1 j )z 1 + j-1 r=2 (c r j-1 + k j c r j + c r-1 j )z r +(k j c j j + c j-1 j )z j + c j j z j+1 = zj+1 , (25)
then, ( 23), ( 24) and ( 25) yields Vj = -

j r=1 k r z2 r + zj (z j-1 + żj + k j zj ) = - j r=1 k r z2 r + zj zj+1 .
i. At final step i, we have to discriminate between systems S ni m and S i 2 in special cases A and B on the one hand, and S i 2 in the generic case C on the other hand.

In the former case, żi = F i (z) + G i (z)ũ with i ∈ {4, m} for system S i 2 in special cases A and B, and for system S ni m , respectively. Since żi =

i-1 r=1 c r i z r+1 + F i (z) + G i (z)ũ, differentiating V i = V i-1 + 1 2 z2 i leads to: Vi = - i-1 r=1 k r z2 r + zi (z i-1 + i-1 r=1 c r i z r+1 +F i (z) + G i (z)ũ). (26)
Using (23) for both j = i-1 and j = i in (26

) yields Vi = - i-1 r=1 k r z2 r + zi c 1 i-1 z 1 + i-1 r=1 (c r i-1 -c r-1 i )z r + c i-1 i z i + F i (z) + G i (z)ũ . (27)
Then, since G i given by either (A.3), (A.5) or ( 19) is a never vanishing function on X0 , using (20) in ( 27) and ( 21) again thus results in: Vi = i r=1 k r z2 r . Hence we obtain the exponential stability of z from Vi being negative definite. Since z and z are linearly related, and since (15)-( 18) define diffeomorphisms, it follows that z = 0 and hence x = 0 are exponentially stable.

In the latter case, i = 3 and the error system is diffeomorphic to S 3 . Thus, we have ż3 = F 3 (z, ξ) + G 3 ũ, so that z3 =

3 r=1 c r 3 z r dynamics are ż3 = c 1 3 z 2 +c 2 3 z 3 +F 3 (z, ξ)+G 3 ũ and differentiating V 3 = V 2 + 1 2 z2
3 with V2 given by ( 22) leads to

V3 = z3 (z 2 + 2 r=1 c r 3 z r+1 + F 3 + G 3 ũ) - 2 r=1 k r z2 r .
Note that G 3 given by (A.6) is a non null constant, so choosing the control law ũ = -G -1 3 3 r=1 c r 4 z r + F 3 (z, ξ) with some positive gain k 3 using (21) thus ensures:

V3 (z) = -k 1 z 2 1 -k 2 z2 2 -k 3 z2 3 . ( 28
)
In this generic case C for system S i 2 , i.e. when d 1 = d 2 and a 1 = a 2 , the system state is (z, ξ).

We have proved that the function V3 (z) is negative semi-definite, hence the manifold Z = {(z, ξ) ∈ R 3 × R : z = 0} is exponentially stable. Along this set, the system evolves on the zero dynamics ξ = f ξ (0, ξ). We have shown in Lemma 1 that the only zero dynamics stable equilibrium point was ξ * a = 0. Using (15d), it results in the asymptotic stability of S x around zero. The effectiveness of the proposed approach is illustrated first on an inertial two-agent microrobotic system S i 2 , considering the three different classes of designs A, B and C. The two inertial microrobots (R 1 ) and (R 2 ) parameters are given in Table 1. Then a non inertial microrobotic system control is considered, for 5 non inertial robots (R k ), whose parameters are given in Table 3 for Re ≪ 1. Magnetic material is the same for all robots with M = 1.72 10 6 A.m -1 . Parameters are such that the inertial and non inertial systems S i 2 and S ni 5 are controllable according to Propositions 1 and 2, respectively.

Simulation results

The pulsatile blood velocity is modeled using the Womersley model [START_REF] Womersley | Method for the calculation of velocity, rate of flow and viscous drag in arteries when the pressure gradient is known[END_REF] which results in a truncated Fourier series approximation. To be consistent with the vessel radius and simplify the study but without loss of generality, a first-order truncated Fourier series is considered with a negative mean value so that the microrobots have to swim against the flow. The physiological datas are given in Tables 2 and4 for a small artery and an arteriole.

v f (t) = a 0 + a 1 cos (ωt + φ) (29)

Inertial Microrobots S i 2

In the following simulations, the Reynolds number lies in the range [0; 2.5] so inertial effects can not be neglected.

The microrobot (R 1 ) is the same in the three simulations, whilst the parameters of the second microrobot are modified in each simulation (see Table 1) so as to illustrate the two special cases A and B, and lastly the generic case C. Initial condition is given by

x 0 = 0.005 -0.0025 0 0 T ∈ X + . (30)
The admissible reference trajectories are computed using [START_REF] Fruchard | Local controllability, trajectory planning, and stabilization of a two-agent underactuated microrobotic system[END_REF] for the same given desired distance d 1,2 . With parameters from Table 1, initial conditions on x r (t) are given in Table 5.

- Reference distance [m] d1,2 2 10 -3 2 10 -3 2 10 -3 Reference state xr0     0 -0.002 0 0         0 -0.002 0.014 0.014         0 -0.002 0 0    

First simulation Case A

In this configuration, the microrobots have the same radius and density, and are consequently equally affected by the drag force since d 1 = d 2 . However, the first one has a high ferromagnetic ratio contrary to the second one which is mainly composed of a therapeutic load. Figure 2a illustrates the convergence of both positions and velocities of the two microrobots to their references after a 2 s long transient phase. During this phase, the control effort is high to bring the system along the admissible reference. The (R 1 ) to (R 2 ) distance, initially set to 7.5 mm, quickly reduces to 1.4 mm before stabilizing at the desired distance d 1,2 . The nonlinear magnetic interaction term δ 1 (x) depends on that distance and severely impacts the system dynamics. To overcome the magnetic interaction which increases as the robots get closer down to 1.4 mm, the control input depicted on Figure 2c has an overshot during the first 0.1 s though not reaching saturation. The system error depicted on Figure 2b is then stabilized to zero. The control input is depicted in Figure 2c: saturation is not reached and the input tends to 0.114 T.m -1 . To follow the same trajectory, had the first microrobot been alone, it would have required only a fifth of the control effort (around 0.025 T.m -1 ); but if the second microrobot had been alone, the control input would have reached 0.835 T.m -1 , i.e. far beyond the actuator limitations, because of its very poor ferromagnetic ratio. The present simulation thus illustrates the interest in using the first robot to tow the second one by taking advantage of the interaction magnetic force to relax the control effort.

Second simulation Case B

In the special case B, the microrobots share the same magnetic sensitivity a to the control input but have not the same drag coefficient since the second microrobot is far smaller than the first one.

The control input reaches both the upper and lower saturations for t ≤ 0.04 s, as can be noticed on Figure 2f. Saturation events do not impact severely the tracking and the control effort is far smaller than in case A over this period. However, contrary to the previous simulation, the robots distance does not go below 1.6 mm so the interaction force remains small, what explains why the control effort in case B is smaller than in case A considering the first 0.2 s, despite the initial saturations. After a 1 s long transient phase, Figure 2e shows that the error is stabilized to zero.

Third simulation Case C

This simulation illustrates the generic situation where the microrobots have no common parameter d nor a.

Table 1 shows that microrobot (R 2 ) is chosen so as to satisfy condition (17) of Lemma 1, e.g. having less magnetic material and a higher drag parameter than the first robot, so that the error system is a minimum phase system with respect to the origin. The simulation is quite similar to the previous one, except that the saturation occurs only at t = 0 s (see Figure 2i). The control input stabilizes the set Z = {(z, ξ) ∈ R 3 × R : z = 0} within 3 seconds and the internal state ξ = x1 -x2 is stabilized due to the minimum phase property. The distance between the two microrobots is given by ξ + d 1,2 . Since ξ is not directly controlled, this distance grows quickly to 20 mm thus inducing a far reduced magnetic interaction between the robots, before it converges to the desired d 1,2 = 2 mm value, as depicted in Figure 2g. In arterioles, the radius and the relative velocity of the microrobots in the blood are smaller (see Table 4, resulting in a Reynolds number Re ≈ 2 10 -3 ≪ 1 for the biggest microrobot (R 5 ) and twice smaller for (R k ), k ≤ 4. The microrobots inertia can be neglected and the 5-agent system is well modeled by a non inertial system S ni 5 . We consider a 5-agent system where the leader (R 5 ) is highly magnetized as can be noticed in Table 3. Thus, thanks to the attractive magnetic interaction, (R 5 ) tows the four smaller and mainly payload filled identical microrobots (R k ), k ≤ 4, along a preplanned trajectory computed using reference distances d k,k+1 and initial state x r0 given in Table 6. The control input is computed using Lemma 2 and Proposition 3 with gains of Table 6. The system state x 0 is initialized at (x r0 -25.10 -6 ) m ∈ X - m . Since the initial error is important and because of the underactuation constraints on the system, the microrobots have to manoeuver to reach the planned trajectory as depicted in Figure 2j. During the transient phase, Figure 2k shows that the (R 5 ) position error x5 peaks up to nearly 4 µm despite a stricly decreasing Lyapunov function. Indeed, the choice of a small gain k 5 makes it possible to use mainly the leader as a degree of freedom to manoeuver the followers: (R 5 ) is highly attractive for (R 4 ) so this attraction dominates the system dynamics. Besides, since (R 5 ) contains more magnetic material, it is more sensitive to the magnetic actuation than other robots and the strategy is here to stabilize x5 after all other states have converged. During the first 0.2 s, (R 4 ) is highly attracted by (R 5 ), so error x4 decreases, which in turn attracts (R 3 ) along the reference trajectory. Then, due to the increasing distance between (R 5 ) and (R 4 ), the resulting magnetic interaction on (R 4 ) becomes negative til 0.4 s, so (R 4 ) is attracted by (R 3 ) whilst the control input changes sign to bring back (R 5 ) closer to the other robots. Figure 2l illustrates the control input: on the very beginning, u < 0 which may seem counter-intuitive since xk < 0, yet since (R 5 ) is the more reactive robot it reduces the distance between (R 5 ) and (R 4 ) while not affecting too much other robots, thus inducing a strong attraction on (R 4 ) to recover the planned trajectory. The next manoeuver is clearly visible during the interval [0.2; 0.4] s when u sign changes to bring back (R 5 ) after peaking. After 1 s, there is no more manoeuver and the tracking convergence looks like the tracking of a fully actuated system. Finally, u converges to u r ≈ 0.1T m -1 where the control input mainly compensates for the pulsatile blood flow. It is worth mentioning that controlling a single microrobot of type (R k ), k ≤ 4, along the same trajectory would have required a 300 bigger control effort (more than 29T m -1 ) than used here to control five microrobots simultaneously.

Non Inertial

Discussion and conclusion

The present paper has presented a controllability analysis and a novel approach to synthesize control laws aiming at Lyapunov stabilizing the trajectory of magnetic microrobots navigating in a blood vessel using the same single control input. To this end, a sufficient and necessary controllability condition is given for this nonlinear underactuated system with drift, which was lacking in the existing literature. In particular, microrobots have to both operate in close vicinity and be designed carefully. A unified Lyapunov stabilizing control law has then been proposed, with concerns about the possible configurations of the microrobotic system. This control synthesis relies on a two steps method: the system is first proven to be diffeomorphic to different nonlinear canonical forms, and a backstepping synthesis coupled with a zero dynamics stability analysis is proposed. Simulations have highlighted the interest in using the magnetic interaction so that one of the microrobots can pull or push the other ones, thus reducing the control efforts required to propel the latter.

Further perspectives, which are currently under investigation, are threefold. First, since therapeutic goals may require that the system has to follow a non admissible trajectory, it would be interesting to analyze under what conditions -at least practical-stabilization can be ensured. Another key point for the control of multiple magnetic microrobots, which has not yet been addressed in the nonlinear case except in [START_REF] Sun | Output feedback synthesis for a two-agent nonlinear microrobotic system[END_REF], is the synthesis of observers in order to propose stabilizing output feedbacks. Lastly, biophysical systems are unlikely to be free from parametric and modeling errors. For instance, in the present paper, the blood velocity is supposed known whilst it is likely that this information is hardly accessible to measurement, thus implying errors in its amplitude, number of harmonics, or frequency. Ongoing work is thus led to ensure stability despite these uncertainties, an issue that we have already addressed for a single microrobot using either adaptive approach [START_REF] Arcese | Adaptive controller and observer for a magnetic microrobot[END_REF] or state estimation in [START_REF] Sadelli | 2D observerbased control of a vascular microrobot[END_REF]. Our recently developed low-power high gain observer for fluid mechanics application [START_REF] Ahmed | Lowpower high gain observers for wake flow rebuild[END_REF] is also promising for estimation of the blood flow despite an unknown blood pulse and inherent possible observability defects. z = 0, it follows from (A.6)-(A.7) that:

f ξ (0, ξ) = (a 1 -a 2 )(a 2 i 1 + a 1 j 1 ) a 1 a 2 (d 2 -d 1 ) 1 (ξ+d1,2) 4 -1 d 4 1,2
.

The equilibrium points of the zero dynamics are thus given by ξ * a = 0 and ξ * b = -2d 1,2 . Linearizing the last equation of S 2 along z = 0, ξ = ξ * gives: ξ = -4(a 1 -a 2 )(a 2 i 1 + a 1 j 1 ) a 1 a 2 (ξ * + d 1,2 (t)) 5 (d 2 -d 1 ) ξ, (A.8) so that the equilibrium points ξ * a and ξ * b are respectively asymptotically stable and unstable provided that (17) is satisfied. So S 2 is a minimum phase system with respect to the equilibrium point (z, ξ) = 0 under the aforementioned condition.

B Proof of Lemma 2

If the assumptions of Proposition 2 are fulfilled, system S ni m is small-time controllable. As in the previous proof, to simplify the control synthesis, a diffeomorphism mapping system S ni m to a pure feedback form is sought. Searching for a function h whose differential is a basis of ∆ = span i=0..m-2 {ad i f g} annihilator leads to (m -1) relationships using the first (m -1) columns of (13). Let h i denote the partial differential of h with respect to xi . Due to the triangular form of the ad i f g(x), one proceeds recursively as follows:

• Using the second relationship dh ad f g(x) = 0 that involves only h m-1 and h m , since (αα m )δ

(1) m-1 (x) = 0 on Xm , one gets a first integral w 1 (x) = xm-1 + ǫm-1 εm-1 xm . So h is invariant on the foliation w 1 being constant.

• Let h w1 = ∂h ∂w1 . Since dh ad 2 f g(x) = 0 involves only h m-2 and h w1 , and δ

(1) m-2 (x) = 0 on Xm , another first integral is given by w 2 (x) = xm-2 + w 1 so h(x 1 , . . . , xm-3 , w 2 ).

• At step i < m, i-th equation dh ad i f g(x) = 0 depends only on h m+1-i and h wi-2 where w j (x) = xm-j + w j-1 , ∀j ≤ i -2. Hence one deduces from δ (1) m+1-i (x) = 0 that h(x 1 , . . . , xm-i , w i-2 ).

• Lastly, using the first relationship dh ad 0 f g(x) = 0 and the fact that δ Consequently, defining h and z as in (18) results in a pure feedback system S 1 with F m = L m f h(x). The mapping ψ : x → z is characterized by a jacobian matrix whose determinant is computed adding all columns in the second one, then factorizing on line 2 ≤ j ≤ m-1 by ǫ j-1 ((m-1)ε m-1 α+ǫ m-1 α m ) and performing a Laplace expansion along the second column to get a lower triangular matrix. The Jacobian is then easily computed as: By controllability assumption in Proposition 2, ̺ m = 0, and since the partial derivatives δ

(1)

i do not cancel on Xm , ψ is a diffeomorphism. For the same reasons, it follows that G m , obtained by induction as ( 19), is non null on Xm , and hence can be inverted.

Fig. 2 .

 2 Fig. 2. From top to bottom: simulations for two inertial microrobots S i 2 in cases A, B and C, and simulation for 5 non inertial microrobots S ni 5 . (a), (d), (g): Microrobots (R1) and (R2) positions (x1, x2) and velocities (x3, x4) are depicted in black, blue, green and red solid lines, their reference trajectories are depicted in dashed lines; (j): Microrobots positions (solid lines) and reference trajectories (dotted lines) for robots (R1), (R2), (R3), (R4) and (R5) are depicted in black, blue, green, cyan and red lines, respectively. (b), (e), (h): State of the error system S i 2 ; (k): State of the error system S ni 5 . (c), (f), (i), (l): the single control input of the systems with saturations ±0.4 T.m -1 depicted in red dash-dotted lines.

  never cancels on Xm , it follows that h depends only on w m-1(x) = -(m -2)αε m-1 + α m ǫ m-1 x1 + αε m-1 w m-2 .

  Microrobots S ni

			5
	Table 6
	Non Inertial Microrobots: Planning & Control Parameters
	Parameter	Symbol	Value
	Control gains	k	5 6 15 19 4	T
	Reference dis-	di,i+1	-50 42 38 101.6	T 10 -6
	tances [m]
	Reference	xr0 -231.6 130 92 50 0	T 10 -6
	state [m]

Acknowledgements

The author is grateful to the Editor and anonymous reviewers for their insightful and helpful suggestions to improve the quality of this work.

A Proof of Lemma 1

Under conditions C1 and C2, it has been shown from (5b), ( 11)-( 12) that dim span i=0..3 {ad i f g}(0) = 4. Now let us seek for a function h whose differential is a basis of the annihilator of ∆ = span i=0..2 {ad i f g}(x):

If such a function can be exhibited, then choosing z i+1 = L i f h(x) maps S i 2 to the pure feedback form S 1 . Let h i = ∂h ∂ xi (x), (A.1) first equation leads to:

). Solving the last two equations of (A.1) will depend on the system parameters.

1 If d 1 = d 2 , i. e. in case A, using (A.2), the second equation of (A.1) gives a 1 h 1 + a 2 h 2 = 0, and using (A.2) in the third equation of (A.1) leads to h 3 = 0. One then deduces that h depends only on (a 2 x1 -a 1 x2 ). Besides, the choice (15a) yields

is a non-vanishing function on any compact. (15a)-(15c) thus defines a diffeomorphism mapping S i 2 on S 1 . Note that, in this particular case, the system controllability is ensured only if a 1 = a 2 so as to satisfy condition C1. A straightforward differentiation yields:

2) in (A.1) last equation gives:

-

. Now let solve the partial differential equations inherited from (A.4):

It is not difficult to show, e.g. by a proof by contradiction, that dh is an exact differential form if and only if a 1 = a 2 = a since it leads to γ = 0. So, for case B, one infers that h is given by (15b). Let ι =

0, the diffeomorphism (15b)-( 15c) maps S i 2 to the feedback normal form S 1 with:

(A.5) 3 In the generic case C, d 1 = d 2 and a 1 = a 2 . In the light of the foregoing analysis, there exists no transformation mapping S i 2 into the pure feedback form. However, the choice of h as in (15a) yields by differentiation to S 2 with:

Now let us study the zero dynamics of the system S 2 , i.e. ξ = f ξ (0, ξ). From (15a)-(15d), it is clear that:


L Florentin
email: florentinlea@gmail.com

Etienne Deloule

F Faure

D Mangin

Chemical 3D-imaging of glass inclusions from allende (CV3

Introduction

Chondrules are silicate blebs embedded, with Ca-Al rich inclusions (CAIs) and metal drops, in a fine-grained matrix of variable porosity. Together, they represent the main components of the undifferentiated meteorites that are chondrites [START_REF] Sears | Overview and classification of meteorites[END_REF]Zanda, 2004, Connolly and[START_REF] Connolly | Chondrules: the canonical and noncanonical views[END_REF]references therein). Chondrules measure between 100 µm and 1 cm in diameter and display a wide range of textures, from glassy to porphyritic. They are thought to have formed by brief and repeated phases of heating [START_REF] Hewins | Temperature conditions for chondrule formation[END_REF]. Porphyritic chondrules are often studied because they correspond to approximately 85% of all chondrule textural types in ordinary chondrites [START_REF] Gooding | Relative abundances of chondrule primary textural types in ordinary chondrites and their bearing on conditions of chondrule formation[END_REF]. They are usually composed of olivines and pyroxenes. They can be classified according to either their redox state or their mineralogy. The focus of this study is the "type I" chondrules, i.e. chondrules whose silicates have low Fe-Ni contents. In this case, the olivines are forsteritic (i.e. Mg-rich). Both PO (Porphyritic Olivine) and POP (Porphyritic Olivine and Pyroxene) chondrules were studied. The origin of type I chondrules is still under debate [START_REF] Grossman | Early chemical history of the solar system[END_REF][START_REF] Wasson | Constraints on chondrule origins[END_REF][START_REF] Yu | Transient heating and chondrule formation: Evidence from sodium loss in flash heating simulation experiments[END_REF][START_REF] Connolly | On the origin of the "kleine Kügelchen" called Chondrules[END_REF][START_REF] Alexander | The formation conditions of chondrules and chondrites[END_REF][START_REF] Chaussidon | Oxygen isotopic constraints on the origin of magnesian chondrules and on the gaseous reservoirs in the early Solar System[END_REF][START_REF] Asphaug | Chondrule formation during planetesimal accretion[END_REF][START_REF] Sanders | The origin of chondrules and chondrites: Debris from lowvelocity impacts between molten planetesimals? Meteorit[END_REF][START_REF] Fedkin | Vapor saturation of sodium: Key to unlocking the origin of chondrules[END_REF][START_REF] Johnson | Impact jetting as the origin of chondrules[END_REF], Florentin et al., 2017). In order to understand how chondrules formed, it is necessary to understand the processes that form the minerals they host, in particular olivines.

Any approach aimed at understanding the formation and behavior of olivine within the chondrules requires study of the glass inclusions they host. These glass inclusions are liquid blebs trapped inside the host minerals during mineral growth [START_REF] Roedder | Origin and significance of magmatic inclusions[END_REF]. Such glass inclusions were observed inside olivines from chondrules in early studies [START_REF] Sorby | On the microscopical structure of meteorites[END_REF][START_REF] Dodd | Metamorphism of the ordinary chondrites: A review[END_REF]. The inclusions are isolated and ovoid and may host a shrinkage bubble and sometimes daughter minerals [START_REF] Fuchs | Mineralogy, Mineral-Chemistry, and Composition of the Murchison (C2) Meteorite[END_REF]. These petrological features are also typical of primary terrestrial melt inclusions [START_REF] Roedder | Significance of Ca-Al-rich silicate melt inclusions in olivine crystals from the Murchison type II carbonaceous chondrite[END_REF], suggesting that inclusions from terrestrial and extraterrestrial samples could behave similarly as a recorder of the growth environment of their olivine hosts.

Experimental work on olivine formation [START_REF] Donaldson | An experimental investigation of olivine morphology[END_REF][START_REF] Jambon | Olivine growth rates in a tholeiitic basalt: An experimental study of melt inclusions in plagioclase[END_REF][START_REF] Pack | Partitioning of Ca and Al between forsterite and silicate melt in dynamic systems with implications for the origin of Ca, Al-rich forsterites in primitive meteorites[END_REF][START_REF] Schiano | Primitive mantle magmas recorded as silicate melt inclusions in igneous minerals[END_REF][START_REF] Faure | Experimental investigation of equilibration conditions during forsterite growth and melt inclusion formation[END_REF] established that melt inclusions (i.e. glass inclusions with a magmatic origin) may not necessarily be representative of the bulk magma from which they originate, depending on the process that controls the olivine growth: diffusion or interface attachment (screw dislocation, surface nucleation or continuous growth mechanisms). When cooling rates are slow, olivine grows by interface attachment processes: diffusion of the components in the melt to and from the melt / olivine interface is fast relative to the uptake and rejection of components at the interface [START_REF] Faure | Experimental investigation of equilibration conditions during forsterite growth and melt inclusion formation[END_REF].

Consequently, the liquid trapped within an olivine to constitute the melt inclusion is representative of the original melt. In contrast, when a crystal grows rapidly, olivine growth is controlled by diffusion, and melt inclusions should trap a compositional boundary layer that develops at the crystal-liquid interface. If so, skeletal or dendritic olivine will be produced. The chemical composition of the inclusions will not in this case be representative of the bulk initial melt composition, even if the inclusion is experimentally re-homogenized. The entrapped melt inclusions will in fact be enriched in slowly diffusing incompatible elements rejected from the olivine structure, such as Al and, to a lesser extent, Ca.

Inclusions hosted in Mg-rich olivines from type I chondrules have been extensively studied over the last few decades [START_REF] Fuchs | Mineralogy, Mineral-Chemistry, and Composition of the Murchison (C2) Meteorite[END_REF][START_REF] Grossman | High-temperature condensates in chondrites and the environment in which they formed[END_REF][START_REF] Kurat | Glass inclusions in Renazzo olivine: condensates from the solar nebula ? 60th Ann. meeting of the Met[END_REF][START_REF] Varela | Chemistry of glass inclusions in olivines of the CR chondrites Renazzo, Acfer 182 and El Djouf 001[END_REF]2005;2006;Varela andKurat, 2009, Faure et al., 2012). Glass inclusions can be separated into four groups based on their chemical composition: Al-rich, Al-poor, Na-rich, and Si-rich [START_REF] Varela | Chemistry of glass inclusions in olivines of the CR chondrites Renazzo, Acfer 182 and El Djouf 001[END_REF], Faure et al., 2017). [START_REF] Varela | Chemistry of glass inclusions in olivines of the CR chondrites Renazzo, Acfer 182 and El Djouf 001[END_REF] considered these chemical compositions to be "exotic", meaning that they are not in equilibrium with their host olivine.

More recently, it was shown that the exotic compositions recorded in Al-rich glass inclusions hosted in Mg-rich olivines could be reproduced experimentally by slow growth of olivine from a CMAS liquid [START_REF] Faure | Origin of glass inclusions hosted in magnesian porphyritic olivines chondrules: Deciphering planetesimal compositions[END_REF]. According to [START_REF] Faure | Contrasted liquid lines of descent revealed by olivine-hosted melt inclusions and the external magma[END_REF], the exotic compositions observed inside synthetic melt inclusions result from metastable crystallization of olivine in near-equilibrium conditions on the walls of the inclusions. This suggests that, like for terrestrial melt inclusions, the original magma compositions of exotic glass inclusions hosted in Mg-rich olivines of chondrules can be experimentally determined or recalculated. However, this requires that, again as with terrestrial inclusions, the natural extraterrestrial inclusions did not trap a compositional boundary layer during their formation.

A study of inclusions from Semarkona type II chondrules by nanoSIMS showed that a twomicron chemical zoning was present at the inclusion-olivine interface [START_REF] Tronche | Melt inclusions in olivine crystals from chondrules: New insights from nanoscale concentration profiles in the Semarkona chondrite[END_REF].

However, no such data exist for type I chondrules, whose origin is presently under debate. To answer this important question, we perform a comparative study of synthetic melt inclusions trapped in olivine grown slowly under controlled conditions [START_REF] Faure | Contrasted liquid lines of descent revealed by olivine-hosted melt inclusions and the external magma[END_REF] and natural glass inclusions from the Allende (CV3) meteorite. Both inclusion types were studied in 3D via SIMS in order to construct 3D chemical images of the inclusions. From these, chemical depth profiles of the inclusions can be drawn, to allow investigation of the presence of a boundary layer at the glass inclusion / host olivine interface.

Though 2D high-resolution imaging via SIMS was developed many years ago [START_REF] Levi-Setti | Aspects of high resolution imaging with a scanning ion microprobe[END_REF], no study has so far applied it to the chemical 3D-imaging of melt inclusions. Today, many techniques allow 2D or 3D chemical imaging, but only the ion probe, using depth profiles, provides, a resolution high enough to distinguish an expected few tens of nanometer scale heterogeneity for non-conductive materials [START_REF] Adams | History and present status of imaging analysis[END_REF]. Furthermore, in contrast with the lateral profiling proposed by Tronche et al., 2007, 3D chemical imaging allows us to ensure that no daughter crystal is present in the inclusion. Consequently, this study is also an opportunity to develop an analytical protocol for mapping inclusions in 3D via SIMS.

Analytical methods

Sample preparation

Three synthetic glass inclusions were selected from olivines experimentally produced by slowly cooling a CMAS liquid in a vertical furnace (charge number 16 in [START_REF] Faure | Contrasted liquid lines of descent revealed by olivine-hosted melt inclusions and the external magma[END_REF].

The starting glass material (47.32 wt. % SiO 2 , 11.96 wt. % Al 2 O 3 , 12.21 wt. % CaO, 28.48 wt. % MgO) was chosen to represent the compositions of the primordial planetesimals that could produce Al-rich glass inclusions hosted in type I porphyritic olivine chondrules [START_REF] Faure | Origin of glass inclusions hosted in magnesian porphyritic olivines chondrules: Deciphering planetesimal compositions[END_REF]. The starting material was first heated for 1 h at 1540°C (4°C above the liquidus temperature) and then cooled slowly at 2°C h -1 to 683°C (below the glass transition).

Three natural glass inclusions were selected from thick sections of Allende (CV3) meteorite.

Two of the glass inclusions were sampled from two distinct olivines hosted in two type I chondrules and the other inclusion came from an isolated olivine. The olivines selected exhibited no clear petrological indications of alteration (no sign of Fe penetration within the chondrule or olivine, Fo# above Fo90). The chemical compositions of the glass inclusions and their host olivine crystals are available in Appendix 1. Careful observation allowed us to select only inclusions that met previously established pristine criteria [START_REF] Roedder | Origin and significance of magmatic inclusions[END_REF]: isolated, with a retraction bubble and an ovoid shape, no visible leak, and no proximity of a fracture. They were also selected only when glassy (i.e. no devitrification or visible daughter crystal) (Figure 1).

Each inclusion was brought to the surface by polishing with corundum discs. For the purpose of 3D mapping via Secondary Ion Mass Spectrometer (SIMS), only inclusions with sizes of between 3 and 20 µm were considered. Each inclusion was then carbon coated and observed using SEM in order to confirm its glassy nature and to roughly estimate its chemical composition. It is worth noting that gold coating is not usually performed on meteorite samples because Au cannot be fully removed. Therefore, both natural and synthetic samples were carbon coated for this study.

SIMS analyses

As the four CMAS oxides (Al 2 O 3 , SiO 2 , MgO and CaO) account for 88 to 98 wt. % of the whole inclusion composition for Allende inclusions (as calculated from major element data for glass inclusions from Allende; [START_REF] Florentin | Origin of Na in glass inclusions hosted in olivine from Allende CV3 and Jbilet Winselwan CM2: implications for chondrule formation[END_REF], only these four elements were analyzed.

Analyses were performed using the CAMECA IMS 1280 at CRPG, Nancy (France) and the CAMECA IMS 7F at IJL, Nancy, (France).

The elements were measured at the masses of 24 each sample in order to select the best signal and minimize the risk of C interferences. For both ion probes, a 10 kV Cs + primary ion beam was used to obtain a small beam size of ≈ 300 nm with a current intensity varying from 130 to 200 pA during measurements. Although a smaller beam would have provided a better lateral resolution, it would not have allowed abrasion of enough matter to measure the whole depth of the inclusion within a reasonable length of time.

The secondary ions images were acquired by ion counting on EM in monocollection mode, with each individual image consisting of 256×256 pixels.

For analysis with the CAMECA IMS 1280, the secondary ions were accelerated at 10kV (impact energy at 20 keV) and a mass resolution at 5000. In contrast, with the CAMECA IMS 7F, the secondary ions were accelerated at 5kV (impact energy at 15 keV). The lower impact energy provides a better depth resolution, but the mass resolution of 4000 for the secondary ions is not enough to separate major elements or oxides from carbon interferences.

Each mass (including the background) was measured for 10 seconds with a waiting time of 1 second. This resulted in a cycle lasting around 2 minutes and sputtering of between 15 and 30 nm depth.

As stated in the previous section, each inclusion was brought to the surface of the sample and then analyzed with SEM. Afterwards, the whole remaining inclusion (i.e. the part below the surface) was sputtered in order to collect 3D chemical images. Each cycle provided a 2D chemical image with an area of between 20x20 and 30x30 µm, depending of the inclusion size, corresponding to sputtering of 30 to 15 nm depth. A chemical 3D image of the sample was then obtained by stacking the 2D images. Depending on the inclusion size and depth, the total number of cycles for one inclusion varied between 200 and 2000. This resulted in acquisitions lasting from two to more than twenty-four hours, depending on the settings of the probe and the size of the inclusion. The acquisition time was estimated from the inclusion size and the settings, but as inclusions were sputtered from thick sections (which did not allow us to observe the 3D shape of the inclusion), each analysis was carefully monitored.

The 3D images and depth profiles were collected using the WinImage® and Wincurve® software provided by CAMECA. For each acquired cycle, (i.e. a 2D-image that represents a few tens of nanometers) each pixel on the image contained counting information for the concerned element. Chemical depth profiles for an inclusion were extracted by compiling, for each cycle, data from an area drawn manually on the image of the inclusion (Figure 2a and2b). After analysis, the depth pits were measured using a Nikon Eclipse LV100 optical microscope, calibrated using a Moticam 3000 Camera.

Each depth profile is presented as an oxide/silica ratio in order to account for eventual internal shift from the probe. All data were also normalized using the oxide/SiO 2 ratio of the host olivine, measured using an SEM (JEOL JSM-6510) equipped with a Brucker Quantax EDS, for consistency. Analyses on EDS were calibrated following the protocol detailed in [START_REF] Florentin | Origin of Na in glass inclusions hosted in olivine from Allende CV3 and Jbilet Winselwan CM2: implications for chondrule formation[END_REF].

Results

Results are presented in the form of 3D images of half glass inclusions (Figure 2b). The color scale of the 3D images corresponds to the number of atoms counted per pixel of the image.

Because of the rounded shape of the inclusion and the beam size of about 300 nm, the periphery of the inclusion is depleted in Al relative to its center due to the contribution of the surrounding olivine, which contains no Al. By selecting an area of 1x1 to 2x2 µm size at the center of the inclusion, it is possible to combine the counting for all the pixels in the area, for each element and for each level between the surface and the bottom of the inclusion. The data are reported as AlO/SiO MgO/SiO and CaO/SiO ratios and are normalized to the oxide/SiO 2 mass ratios measured via SEM for the olivine. In this way, it was possible to draw depth profiles for the oxide/SiO 2 ratios. It is worth noting that for all inclusions, the duration and the number of acquired cycles of the depth profile depends on the size of the inclusion.

Profiles from all inclusions typically display three distinct segments. Based on the 3D-imaging, the first can be attributed to the inclusion and the last to olivine The intermediate segment corresponds to the transition zone between glass and crystal.

Study of chemical profiles at glass inclusion-host olivine interface from synthetic samples

Details of 3D images and depth profiles for each synthetic inclusion are available in Appendix 2.

(For SEM and optical microscopic images of these inclusions, see [START_REF] Faure | Contrasted liquid lines of descent revealed by olivine-hosted melt inclusions and the external magma[END_REF].

Typical depth profiles for two synthetic inclusions are presented in Figure 3, one measured with the IMS 7F ion probe, the other with the IMS 1280 ion probe. Both inclusions are hosted in the same synthetic olivine. The two profiles show three distinct segments.

Inclusion 1, analyzed with the IMS 7F ion probe, shows a nearly flat CaO profile with a CaO/SiO 2 ratio of around 0.01. This flat pattern is most likely due to isobaric interferences between CaO and C 2 O 2, caused by carbon coating, whose signals were not separated due to the low mass resolution of the IMS 7F settings. The Al 2 O 3 /SiO 2 ratio decreases from 0.1 to 0.002, and is higher in the melt inclusion than in the olivine. MgO measurements were not used because of the isobaric interference of C 2 O on MgO due to the carbon coating of the section. Because of the low mass resolution of the IMS 7F setting, it was not possible to separate the two signals. For this inclusion (inclusion 1), the transition zone between glass and olivine is not clearly separated from the melt inclusion signal. This is due to the small size of the inclusion: the sputtering removed 2.5 µm in depth but the bottom of the inclusion was at around 1.5 µm depth below the surface. The signal, especially for Al 2 O 3 /SiO 2 , does not correspond to a straight line that drops when the end of the inclusion is reached, but rather to a smooth slope towards the end of the inclusion. As the inclusion is small and rounded (i.e. the bottom of the inclusion is not flat), the signal recorded from the inclusion decreases gradually. The depth profile for the Al 2 O 3 /SiO 2 ratio thus shows a slight decrease towards the glass/olivine transition. The profile is flat in the olivine and shows a steeper slope at the interface. Inclusion 2 was analyzed with the IMS 1280 ion probe. The CaO/SiO 2 ratio remains around 0.03 and slowly decreases towards the bottom of the inclusion. This value is higher than the value measured in inclusion 1, even though they contain the same amounts of CaO and SiO 2 . This slight difference is most likely due to the fact that the measurements were made using two different ion probes. The MgO/SiO 2 ratio for inclusion 2 is clearly an artifact due to a matrix effect: the emission yield of MgO relative to SiO is different and is higher for a glass than for an olivine.

The Al 2 O 3 /SiO 2 ratio depth profile shows the three distinct segments: a gently decreasing slope inside the inclusion, a flat profile in the olivine, and a steep slope in the transition zone.

If we take into account the 1 µm size difference between the two inclusions, these results show that the depth profiles are similar for both inclusions in spite of the use of two different SIMS.

Study of chemical profiles in inclusions hosted in olivine from Allende chondrite

Details of 3D images and depth profiles for natural glass inclusions from Allende are available in Appendix 2. The depth profiles measured on two inclusions are presented in Figure 4. These are representative of the analyses of three Allende glass inclusions. Inclusions were analyzed from two distinct olivines with the 1280 ion probe. The chemical depth profiles show the typical three distinct segments: largely flat in the glass inclusion, flat in the olivine, and sloped in the steep transition zone. Both inclusions show the same ratios of MgO/SiO 2 and CaO/SiO 2 and the inclusion "Allende 1" shows a lower Al 2 O 3 /SiO 2 ratio than the inclusion "Allende 2". This reflects the large compositional range of Allende inclusions (Florentin et al., 2016), and is consistent with Al-rich and Al-poor glass inclusions groups defined previously [START_REF] Varela | Chemistry of glass inclusions in olivines of the CR chondrites Renazzo, Acfer 182 and El Djouf 001[END_REF]. It can also be observed that these profiles have a much better definition than the profiles of synthetic inclusions. This is due to the diameter of the inclusions (15 to 20 µm for natural inclusions, 3 to 5 µm for synthetic inclusions).

Likewise, the slopes at the transition zone of the natural inclusions seem to be much steeper than those of the synthetic inclusions, but this is an effect of the total depth of analysis.

Data processing

Expected chemical compositions at the glass/olivine interface

The elemental profiles recorded in the glass inclusions depend on the processes involved in the olivine growth. Three types of chemical profile can be expected at the periphery of the inclusion (Figure 5). (i) If olivine growth is controlled by interface processes, a flat profile will be expected in the glass (Figure 5a). In this case, the olivine and the liquid are at equilibrium or near equilibrium conditions at all times. (ii) If the olivine grows rapidly, the growth is therefore controlled by diffusion processes and an increase in incompatible elements in the glass should be recorded towards the glass-olivine transition zone (Figure 5b). In this case, Al and Ca, both incompatible elements with regard to olivine, are not incorporated into the olivine and a boundary layer can be observed at the periphery of the inclusion. (iii) Finally, when melting of the olivine walls occurs around the inclusion, a decrease in the incompatible elements should be observed due to their enrichment in the olivine component (Figure 5c). In this case, Mg and Si would be added to the liquid composition, which would increase the MgO content and decrease the Al 2 O 3 and CaO contents in the glass towards the transition zone.

It is now possible to compare these theoretical profiles with our measured profiles from the synthetic and natural inclusions. However, it is worth noting that theoretical profiles, unlike measured ones, are free of any matrix effect. The measured profiles for MgO, CaO, AlO and SiO, depend on the matrix (olivine or glass) they are measured in. Consequently, the ratios cannot be quantitatively representative of the real chemical ratios within the inclusion or the olivines, which are shown for the theoretical profiles.

Our results show no increase in incompatible elements, which rules out the possibility of an incompatible-rich boundary layer due to rapid crystallization trapped in the inclusion. They instead suggest a decrease in incompatible elements.

This could indicate reheating of the inclusion after trapping and melting of olivine on the walls of the inclusion. However, several arguments suggest this is not a reasonable hypothesis. (1) Synthetic and natural inclusions both contain only a few wt. % MgO. In the case of olivine melting, the MgO content of the inclusion would increase dramatically, which is not observed here. (2) Both synthetic and natural glass inclusions show very similar profiles. However, the thermal history of synthetic inclusions is well known, and they were not reheated after the trapping of the inclusions [START_REF] Faure | Contrasted liquid lines of descent revealed by olivine-hosted melt inclusions and the external magma[END_REF]. (3) Calcium diffuses more than two times faster than Al in silicate melts [START_REF] Liang | Diffusion in silicate melts: I. Self diffusion in CaO-Al 2 O 3 -SiO 2 at 1500°C and 1 GPa[END_REF]. Consequently, the Al 2 O 3 /CaO ratio is expected to show the existence of an incompatible-rich boundary layer, as the amount of Al 2 O 3 should rise relative to the amount of CaO near the transition zone. Thus, in the case of such a boundary layer, the Al 2 O 3 /CaO ratio should increase towards the wall of the inclusion due to diffusion. When compared, Al 2 O 3 /CaO ratios for synthetic inclusions and Allende natural glass inclusions (Figure 6) show a similar evolution at the transition zone, with only a slight dependence on the inclusion size.

Analytical artifact or real chemical transition?

For synthetic inclusions, Al 2 O 3 /SiO 2 ratio patterns display two distinct segments from the center of the inclusion towards the border of the olivine: a nearly flat profile inside the inclusion and a steep slope at the transition zone (Figure 3). However, the corresponding theoretical profile (Figure 5a) shows no transition zone. The transition zone can thus be attributed to analytical artifacts. (1) The main part of this artifact is the shape of the glass inclusion. As the bottom of the inclusion is not flat but rounded, the abundances (recorded over an area) would be expected to decrease slowly towards the end of the inclusion, as a growing part of the area becomes olivine (for details, see Appendix 3). (2) Another cause of the artifact is due to a technical limitation, in that depth profiles are calculated from an area rather than a discrete point. As seen in Figure 7, the size of the mixing zone between the olivine and the inclusion can be related to the size of the area (Figure 7a) used to draw depth profiles from the stacking of 2D images: the smaller the area, the smaller the mixing zone (Figure 7c). However, a small area corresponds to a mean calculated from fewer pixels, which increases the noise in the depth profile (see Areas 1 and 2 in Figure 7b). Consequently, in order to obtain a well-defined profile, it is necessary to use an area that is wide enough to prevent noise but small enough to achieve the smallest possible mixing zone. (3) Furthermore, it was established that when sputtering, the beam implants some atoms further into the sample instead of removing them [START_REF] Wilson | Secondary Ion Mass Spectrometry[END_REF]. The tail effect associated with this ion mixing is usually considered to be a hundred nanometers in size, which is much smaller than the mixing area, even if in our case the very high contrast of Al contents in the melt inclusion and in the olivine might enhance the tail depth.

Consequently, the progressive decrease in incompatible elements observed at the olivine-glass inclusion interface is the result of an analytical artifact and does not correspond to an incompatible-rich boundary layer. As such, it will be considered an artifact in the following discussion.

Discussion

Profiles comparison between synthetic glass inclusions and natural glass inclusions from

Allende

At the transition zone, the Al 2 O 3 /CaO ratio depth profiles from synthetic and natural melt inclusions (Figure 6) show similar patterns: a slight decrease when reaching the inclusion walls, a steep decrease in the transition zone, and then a flat profile in the olivine. The synt1 and synt 2 profiles resemble each other closely but are not identical. This may be attributed to the fact that the inclusions are not of the same size, meaning that the size of the transition zones are different.

It could also be due to the use of different probes with different settings. The two profiles obtained from the Allende inclusions (All1 and All2) also do not superimpose. In this case, the differences can be attributed to the wide compositional range of the glass inclusions in natural samples. Nevertheless, all four profiles display the same overall trend, with a transition zone smaller than 2 µm. The nature of this transition zone will be discussed in the following sections.

Implications for the absence of an incompatible-rich boundary layer

Al 2 O 3 /CaO ratios (Figure 6) suggest that both natural and synthetic glass inclusions experienced a similar process of formation. If true, then the inclusions were not affected by the complex history of the meteorite. Several reasons can be invoked to explain this. First, chondrules do not record the occurrence of any high temperature event after their formation. The formation of the parent body and potential alteration events took place at much lower temperatures, and thus only affected the parent body on a local scale. As the olivines were selected from Allende thick sections that showed no petrological evidence of alteration, we can assume that they did not suffer significantly from low temperature alteration. Consequently, the glass inclusions, sealed within the olivine, should not have suffered any low temperature alteration. Collisions between chondrite parent bodies should also leave their marks, especially in terms of thermal changes.

For this reason, olivines were also only selected when they looked pristine and when the finegrained matrix surrounding the chondrules was non-metamorphosed.

Analogy between synthetic and natural glass inclusions suggests that glass inclusions hosted in olivines from the Allende chondrite could have formed by slow crystallization of a molten planetesimal, as proposed by [START_REF] Faure | Origin of glass inclusions hosted in magnesian porphyritic olivines chondrules: Deciphering planetesimal compositions[END_REF]. In this case, the olivine crystals would have grown slowly enough to trap, in the inclusions, the chemical composition of their growth environment. Thus, any post-trapping processes such as daughter mineral crystallization, olivine overgrowth, or formation of a shrinkage bubble, should be reversible. The problem of iron loss sometimes observed for terrestrial melt inclusions [START_REF] Danyushevsky | Experimental and petrological studies of melt inclusions in phenocrysts from mantle-derived magmas: an overview of techniques, advantages and complications[END_REF] should not be an issue here as the olivine crystals are extremely magnesian-rich and were extracted from type I chondrules in which Fe is in a reduced form. This suggests that experimental homogenization of these inclusions should allow the retrieval of the composition of the parent magma to olivine crystals (Tissandier et al., 2017).

Insight on olivine growth rate and chondrule cooling rate

The Al 2 O 3 /CaO ratio within glass inclusions hosted in chondrule olivines has been investigated in previous studies [START_REF] Varela | Chemistry of glass inclusions in olivines of the CR chondrites Renazzo, Acfer 182 and El Djouf 001[END_REF][START_REF] Varela | A liquid-supported condensation of major minerals in the solar nebula: Evidence from glasses in the Kaba (CV3) chondrite[END_REF]. The majority of glass inclusions record a chondritic ratio in spite of Ca and Al disparities from one inclusion to another. The authors suggest that this chondritic ratio, in association with Al and Ca variability, corresponds to a diffusion process, and propose a rapid crystallization. Data from Varela and coworkers (2009 and reference therein) show that some inclusions also contain high amounts of Na 2 O (up to 15 wt. %). Data show that these Na-rich inclusions are also depleted in CaO and record a different Al 2 O 3 /CaO ratio to the chondritic ratio. However, if diffusion is due to rapid crystallization, Ca and Na should diffuse in the same direction. Hence, if there is an increase in Na, an increase in Ca should also be observed. This rapid crystallization process thus remains questionable.

In the present study, no chemical gradient is observed at the interface between inclusions and olivines hosted in chondrules. This suggests that the olivine overgrowth was not rapid and was therefore not controlled by diffusion processes. Consequently, the last high-temperature cooling episode of the olivine, which probably corresponds to formation of the chondrule, was slow enough that the olivine overgrowth, inside the inclusion, was not controlled by diffusion. Thus, the cooling of the chondrule was slower than Ca and Al diffusion. Furthermore, Ca zoning would be expected within the olivine overgrowth because the Ca content of liquid inclusion increases during crystallization [START_REF] Libourel | Systematics of calcium partitioning between olivine and silicate melt: implications for melt structure and calcium content of magmatic olivines[END_REF][START_REF] Spandler | Diffusion and partition coefficients of minor and trace elements in San Carlos olivine at 1,300°C with some geochemical implications[END_REF]. However, no such variation can be observed in the olivine surrounding the inclusions, neither in the experimental nor the Allende olivines (see Figure 3). Consequently, the cooling rate of the olivine (and likely that of the whole chondrule) must have been at least as slow as the cooling rate used in the experiments (2°C/h) and inevitably slower than the diffusion rate of Al in the liquid and Ca in olivine to erase all chemical profiles. Furthermore, the absence of an incompatible-rich boundary layer at the olivine / glass inclusion interface also allows the inclusion to be considered as having recorded the environment of host-olivine growth. As detailed in the introduction to this paper, if an inclusion has not recorded an incompatible-rich boundary layer, its chemical composition can yield information on the growth environment and / or processes of the host olivine.

As explained earlier in this paper, Ca and Al have different diffusion rates. If we consider that olivines were formed in a molten planetesimal with a chondritic composition (as suggested in [START_REF] Florentin | Origin of Na in glass inclusions hosted in olivine from Allende CV3 and Jbilet Winselwan CM2: implications for chondrule formation[END_REF], then the CaO/Al 2 O 3 ratio should be close to the chondritic ratio as Ca and Al are both incompatible elements. Consequently, if olivine growth had been controlled by diffusion, i.e. that it had undergone rapid crystallization, then the CaO/Al 2 O 3 ratio should be fractionated and thus not correspond to the chondritic ratio. The CaO/Al 2 O 3 ratio remains close to the chondritic ratio (≈ 0.7 when calculated with data from [START_REF] Anders | Abundances of the elements: Meteoritic and solar[END_REF] for both the glass inclusions measured in this study (0.64) and the majority of glass inclusions analyzed in the literature (0.79) (see Florentin et al., 2017 for a compilation). This suggests that there was also no incompatible-rich boundary layer in the inclusion when it was first trapped during olivine formation. Thus, as for chondrules, olivine growth in a molten planetesimal imposes a slow cooling rate. Similarly, the cooling rate of the olivine in the molten planetesimal is necessary low and, in particular, is slower than Ca and Al diffusion rates. It should therefore theoretically be possible to determine the minimum time (i.e. a maximum cooling rate) required to erase Ca zoning on the rim of olivine overgrowth around the inclusion. Diffusion coefficients are not easily calculated. Experiments were carried out by [START_REF] Coogan | Experimental determination of the diffusion coefficient for calcium in olivine between 900°C and 1500°C[END_REF] on both San Carlos olivines and forsterites in order to determine the relationship between fO 2 , temperature T, and diffusion coefficient D for Ca along the three axes of the olivine ([100],

[010], [001]), as a function of temperature and oxygen fugacity fO 2 . The authors proposed the best fits for pre-exponential term Do and E for each axis. The diffusion coefficient D at a given T and fO 2 is then defined by the Arrhenius equation (1): log (D) = log [Do e (-E/RT) ] + 0.31 Δ log fO 2 (1) However, this equation is dependent on many variable parameters. The fO 2 can be determined from previous studies of type I chondrules, rich in magnesian olivine (log (fO 2 ) from -9 to -13 according to [START_REF] Lauretta | Opaque minerals in the matrix of Bishumpur (LL3.1) chondrite: Constraints on the chondrule formation environment[END_REF]. D is strongly dependent on temperature. In this study, we focus on a temperature of 900 °C, close to the end of olivine overgrowth crystallization. This temperature was deduced from the compositions of glass inclusions (Faure and Tissandier, For Coogan and co-workers, E was calculated to be 193 kJ/mol but previous papers showed that for highly magnesian olivines, E could be as high as 400 kJ/mol [START_REF] Hain | 42 Ca Tracer diffusion in Olivin in Abhängigkeit von der Temperatur, dem Sauerstoffpartialdruck und der chemischen Zusammensetzung[END_REF][START_REF] Morioka | Cation diffusion in olivine -II. Ni-Mg, Mn-Mg, Mg and Ca[END_REF]. The variability of E has a critical effect on the final value of D. Thus, it seems difficult to establish a precise D value for our samples, but typical cooling times can be estimated from the available range of D, and the duration of Ca diffusion along 1 µm of olivine, the estimated size of the overgrowth rim around inclusion, can then be calculated, from equation (2):

σ =√Dt (2)
where σ is the typical diffusion distance, set to 1 µm.

Diffusion is a dynamic process, and is fast at high temperature, and slow at low temperature, as suggested by the large range of diffusion coefficients: at low temperature with respect to olivine melt, D ranges from 1.4 10 -17 m²/s at 1200°C [START_REF] Morioka | Cation diffusion in olivine -II. Ni-Mg, Mn-Mg, Mg and Ca[END_REF] to 2.0 10 -20 m²/s at 900°C [START_REF] Coogan | Experimental determination of the diffusion coefficient for calcium in olivine between 900°C and 1500°C[END_REF]. For log(D) values in the range of -20 to -18, Ca diffusion along 1 µm will require between 10 6 and 10 8 s (12 days to 3 years). The longer times are probably overestimated, as our synthetic olivines were formed in 21 days [START_REF] Faure | Contrasted liquid lines of descent revealed by olivine-hosted melt inclusions and the external magma[END_REF] and do not show any limit layer.

We therefore believe that the growth of the natural olivine occurred under conditions close to those of the experimental olivine growth described in Faure and Tissandier, 2014, with a cooling rate of 2°C/h.

To summarize, the olivines appear to have grown slowly via interface attachment processes in a magma from a molten planetesimal and have trapped glass inclusions that recorded information on the olivine surroundings (Figure 8a). Afterwards, an impact event [START_REF] Sanders | The origin of chondrules and chondrites: Debris from lowvelocity impacts between molten planetesimals? Meteorit[END_REF] ejected these olivines into the protoplanetary disk, along with gas, liquid and dust (Figure 8b).

Interactions between gas, liquids and solid resulted in the formation of proto-chondrules (Figure 8c). It is now acknowledged that proto-chondrules underwent several heating episodes (Zanda, 2004 and references therein), as did the olivines they host (Figure 8d). Consequently, heating episodes allowed the olivine overgrowth inside glass inclusions to melt and re-crystallize sequentially. As no incompatible-rich boundary layer is recorded in the present glass inclusions, this suggests that the last cooling episode to have affected the olivines (and therefore the chondrules) was slow enough to have allowed Ca diffusion in olivine.

Conclusion

This study presents the first high resolution (a few tens of nm) chemical 3D-images and depth profiles of glass inclusions. Comparison of synthetic inclusions trapped in slow growing olivine via a magmatic process and of natural glass inclusions from Allende chondrite shows that synthetic and natural samples behave similarly, especially concerning the absence of an incompatible-rich boundary layer at the walls of inclusions. This suggests that both types of olivine (synthetic or from chondrules) were formed by the same process and that their growth was controlled by interface attachment processes. This in turn suggests that olivines, and likely, chondrules, cooled more slowly than the diffusion rate of Al in melt from glass inclusions, allowing olivine overgrowth to form slowly and maintain a metastable equilibrium with its surrounding glass inclusion. Furthermore, as the inclusions display a chondritic CaO/Al 2 O 3 ratio, this means that no diffusion process occurred during growth of olivine. Consequently, the host olivine probably grew slowly, at a cooling rate slower than Ca and Al diffusion rates.

Captions for each figure: In each profile, the transition zone corresponds to the interface between the olivine and the inclusion. In profiles (b) and (c), the transition zone is marked by a chemical variation in incompatible elements. In profile (c), the transition zone is due to dissolution of olivine into the melt at the interface.

Figure 1 :

 1 Figure 1: SEM (BSE) image of glass inclusion Allende 3 inside its host olivine.

Figure 2 :

 2 Figure 2: Microscopic photography of a synthetic melt inclusion and 3D-images of AlO inside the inclusion obtained via SIMS.The color scale for the 3D images corresponds to the number of atoms counted per pixel of the image. For all three images, the scale is in µ m. a: Synthetic inclusion observed in transmitted light. Next to it is the matching surface image (one cycle of acquisition) measured for Al by CAMECA IMS-1280. The atom counting area used to draw the corresponding depth profile is also shown (black square). b. 3D-image of AlO for the synthetic inclusion shown in Figure1a. The bubble appears as a slight decrease in AlO. c: 3D-image of AlO for the same inclusion, slightly rotated, and including the volume corresponding to the area used to draw a depth profile.

Figure 3 :

 3 Figure 3: Depth profiles for Al 2 O 3 /SiO 2 , CaO/SiO 2 and MgO/SiO 2 ratios in two synthetic inclusions. Inclusion 1 was measured with 7F IMS (top) and Inclusion 2 with 1280 IMS (bottom). The MgO/SiO 2 ratio can be observed for inclusion 2. Both Al 2 O 3 /SiO 2 ratio profiles show the typical three segments: a flat segment for the olivine, a smooth slope for the glass inclusion and a steep slope the transition zone.

Figure 4 :

 4 Figure 4: Depth profiles for two inclusions from Allende chondrules. Depth profiles for Al 2 O 3 /SiO 2 , CaO/SiO 2 and MgO/SiO 2 ratios from two glass inclusions hosted in Allende chondrules, Allende 1 (a) and Allende 2 (b.). Both were measured with 1280 ion probe.

Figure 5 :

 5 Figure 5: Sketches of theoretical depth profiles. Sketches of theoretical depth profiles of Al 2 O 3 /SiO 2 (red), CaO/SiO 2 (green) and MgO/SiO 2 (blue) in glassy inclusion (GI) and in olivine. (a.) olivine growth by interface attachment processes, (b.), rapid crystallization of the crystal and rejection of incompatible elements from the crystal in the liquid and (c.) fusion of the host olivine.In each profile, the transition zone corresponds to the interface between the olivine and the inclusion. In profiles (b) and (c), the transition zone is marked by a chemical variation in incompatible elements. In profile (c), the transition zone is due to dissolution of olivine into the melt at the interface.

Figure 6 :

 6 Figure 6: Shapes of Al 2 O 3 /CaO ratio depth profiles near the inclusion-olivine interface for both inclusion types. Shapes of Al 2 O 3 /CaO ratio depth profiles near the inclusion-olivine interface in the two glass inclusions from Allende (All 1 and All 2) presented in Figure 3 and in the two synthetic inclusions (Synth 1 and

  Mg 16 O, 28 Si 16 O, 27 Al 16 O, 40 Ca 16 O for analyses with the CAMECA IMS 7f. For the CAMECA IMS 1280, elements were measured at the masses of 26 Mg, 27 Al, 30 Si, 40 Ca and 24 Mg 16 O, 28 Si 16 O, 27 Al 16 O, 40 Ca 16 O. All masses were measured for

Acknowledgements

We want to thank our editor Sara Russell for her insightful comments. We also wish to thank Alan Rubin, Dominique Hezel and an anonymous reviewer for their thorough reviews and constructive advice. Many thanks to Alice Williams for revising the English. This study was supported by the Programme National de Planétologie: PNP-INSU.

References:

Synth 2) presented in Figure 2. Only the inclusion-olivine interface is represented here to allow comparison of the behavior of Al 2 O 3 /CaO ratios. a: Trapping of glass inclusions in slowly growing olivine within a molten planetesimal. b: Ejection of olivines, gas, liquid and more generally dust during an impact event. c: Interaction between olivine and the surrounding liquid and gas in order to form proto-chondrules. d: After multiple heating episodes of the chondrule, a final phase of slow cooling allows the olivine overgrowth to grow back without forming a chemical boundary layer at the olivine / glass inclusion interface.