Laurent Arcese
email: laurent.arcese@bourges.univ-orleans.fr

Ali Cherry
email: ali.cherry@ensi-bourges.fr

Matthieu Fruchard
email: matthieu.fruchard@bourges.univ-orleans.fr

Antoine Ferreira
email: antoine.ferreira@ensi-bourges.fr

Optimal trajectory for a microrobot navigating in blood vessels

The chemotherapy magnetically controlled under Magnetic Resonance Imaging (MRI) is currently one of the active areas of cancer research. This paper proposes a precise model of a therapeutic microrobot magnetically steered in blood vessels. This modeling approach takes into account the nonnewtonian behavior of blood, as well as wall effect on the blood's profile and robot-to-wall interaction forces. A backstepping approach law is used to ensure a null error between the real trajectory and an optimal reference trajectory deduced from the highly nonlinear model. The strengths and limitations of the overall study are evaluated by simulations.

I. INTRODUCTION

Cancer is a class of diseases in which a group of cells display uncontrolled growth, invasion, and sometimes metastasis. Since complete eradication of cancer cells is imperative, total excision is the best treatment. However, depending on the tumor location and the damage caused to surrounding tissues, surgery is not always possible. In these conditions, chemo or radiotherapy becomes necessary. The drawback of these treatments is that they do not differentiate between healthy and diseased cells. Therefore, the development of techniques that could selectively deliver drug molecules to the diseased site, and among them the chemotherapy magnetically controlled under MRI, is currently a promising area of cancer research [START_REF] Mathieu | Method of propulsion of a ferromagnetic core in the cardiovascular system through magnetic gradients generated by an mri system[END_REF]. Such precision targeting via therapeutic devices will reduce treatment side effects, resulting in better patient compliance. This novel technique relies on generating thrust on a magnetic device, containing therapeutical solution, using the magnetic force related to the gradients of the MRI magnetic field. Not only does MRI device provide propelling energy, but accurate observation capability, thanks to control and imaging multiplexing approach [START_REF] Tamaz | Real-time mri-based control of a ferromagnetic core for endovascular navigation[END_REF].

Because the tiniest capillaries are in the 5-6 micron range, the robot should be in the range of 10-100 nanometers [START_REF] Bawa | Nanoparticle-based therapeutics in humans : A survey[END_REF] to avoid embolization and to drive the drugs as close as possible to the tumor. Analysis on the force capabilities of clinical MRI systems shows that it is impossible to steer nanoscale robots in arteries. An optimal ratio between robot's and vessel's radii, evocated in [START_REF] Mathieu | Method of propulsion of a ferromagnetic core in the cardiovascular system through magnetic gradients generated by an mri system[END_REF], should be respected. A way to overcome the MRI limitations is to make the radius of the robot decrease as it goes in smaller vasculature. This can be done if the robot is a polymer binded aggregate of magnetic particles, loosing substance with time.

This paper proposes a precise model of a robot in blood vessels including wall effects (parabolic profile of blood flow, pulsatile vessel walls), wall interactions (Van der Waals, electrostatic and contact forces) and non-Newtonian behavior of blood (Section II). This model points out the importance of wall interactions, often neglected, in the determination of an optimal trajectory (Section III) which minimize significantly the control efforts. A backstepping control approach [8] [6] is used to ensure a stabilization along any desired trajectory. The performance and the stability of the closed loop system with respect to noise measurement and parameters variations are illustrated by simulations (Section IV). Finally, we sum up results and discuss prospects (Section V).

II. MODELING

The purpose of this section is to present a 2D highly nonlinear model for a polymer binded aggregate of ferromagnetic (NdFeB) particles immersed in blood vessel. The model encompasses the different forces that affect the robot's motion as well as its interaction with the vessel wall. The translational of the robot is expressed by:

m d⃗ v dt = ⃗ F m + ⃗ F d + ⃗ W a + ⃗ F c + ⃗ F vdw + ⃗ F elec (1)
where ⃗ v is the translational velocity of the robot and m its mass. ⃗ F m , ⃗ F d , ⃗ W a , ⃗ F c , ⃗ F vdw and ⃗ F elec respectively denote the magnetic force produced by the MRI gradient coils, blood hydrodynamic drag force, apparent weight, the robot-to-wall contact force, the Van der Waals and the electrostatic force.

In the rest of this paper, we assume that the orientation of the robot does not change due to the magnetic torque which tends to align the magnetization of the robot along the external field. We also assume that the robot is large enough to neglect the effect of Brownian motion.

A. Magnetic force

The gradient coils of the MRI system provide magnetic gradients which produce a magnetic force F m on the robot:

⃗ F m = τ m µ 0 V ( ⃗ M .∇) ⃗ H ( 2
)
where V is the robot's total volume, τ m = Vm V with V m the ferromagnetic volume, ⃗ M is the magnetization of the material, µ 0 is the permeability of free space, ⃗ H is the external magnetic field, and ∇ is the gradient operator.

B. Hydrodynamic drag force

The hydrodynamic drag force ⃗ F d exerting on a spherical body is expressed as:

⃗ F d = - 1 2 ρ f [ ∥(⃗ v -⃗ v f )∥ β ] 2 AC d (⃗ v -⃗ v f ) ∥(⃗ v -⃗ v f )∥ (3)
where ⃗ v -⃗ v f denotes the relative velocity of the robot with respect to the fluid, A is the frontal area of the core and ρ f is the density of the fluid.

For endovascular applications, influence of the vessel walls on the velocity of the robot has to be taken into account. In general, this wall effect is expressed as [START_REF] Kehlenbeck | Empirical relationships for the terminal settling velocity of spheres un cylindrical columns[END_REF]:

β = 1 + ( λ λ0
) α0

-λ α0

with ratio λ = 2r/D and D denoting the vessel diameter (in meter). Parameters α 0 and λ 0 are functions of Reynolds number, but are commonly set to 1.5 and 0.29, respectively. The drag coefficient C d , which is a function of the Reynolds number, is given by [START_REF] White | Viscous Fluid Flow[END_REF]:

C d = 24 Re + 6 1 + √ Re + 0.4
with Re denoting the Reynolds number:

Re = 2rρ f |⃗ v -⃗ v f | βη
In the case of blood, which exhibits a non-newtonian behavior, the fluid's viscosity η is a function of vessel diameter d (in micron) and hematocrit rate h d according to the following empirical relations [START_REF] Pries | Biophysical aspects of blood flow in the microvasculature[END_REF]:

η = η plasma d 2 (d -1.1) 2 [ 1 + (η 0.45 -1)d 2 (d -1.1) 2 (1 -h d ) c -1 (1 -0.45) c -1 ]
with parameters η plasma and η 0.45 denoting respectively the plasma's viscosity and the relative apparent blood viscosity for a fixed discharge hematocrit of 0.45, given by:

η 0.45 = 6 e -0.085d + 3.2 -2.44 e -0.06d 0.645
The shape of the viscosity dependance on hematocrit is:

c = 10

)

Wall effects on the fluid in a vessel traditionally result in a parabolic profile of blood flow (see Figure 1). Besides, to fully take into account pulsatile flow caused by heart pumping in arteries, one has to consider a periodic deformation of the vessel's diameter D(t) synchronized with the pulsative blood velocity v f (t).

C. Apparent weight

In addition to the magnetic and hydrodynamic forces, apparent weight (combined action of weight and buoyancy) is acting on the spherical robot:

⃗ W a = V (ρ -ρ f )⃗ g ( 4
)
where ρ = τ m ρ m + (1 -τ m )ρ poly with ρ m and ρ poly the magnetic material's and polymer's densities.

D. Contact force

The normal and tangential interactions between the robot and the wall are assumed to be expressed by a Voigt model (see figure 2), where the spring constant K and the decay coefficient of dashpot b are functions of Young's modulus E and Poisson's ratio σ [START_REF] Tsuji | Lagragian numerical simulation of plug flow of cohesionless particles in a horizontal pipe[END_REF].

The normal component of the contact force ⃗ F cn acting on the robot is given by:

⃗ F cn = (K n δ 3/2 + b n δ)⃗ n ( 5
)
where δ is the elastic deformation of the wall at the contact point and ⃗ n is the normal unit vector pointing from the robot to the contact surface.

The tangential component of the contact force takes part when the robot is rotating or in case of oblique collision with the wall:

⃗ F ct = (K t ζ + b t ζ) ⃗ t (6)
where ζ is the displacement in the tangential direction and ⃗ t is the tangential unit vector.

E. Van der Waals and electrostatic forces

When the robot and the wall are not in contact, they interact each other through Van der Waals and electrostatic forces. The Van der Waals potential between the robot and the wall is given by [START_REF] Iimura | Simulation of entrainment of agglomerates from plate surfaces by shear flows[END_REF]:

⃗ V vdw = - A h 6 ( 1 h + 1 2 + h + ln h 2 + h ) ⃗ n (7)
where A h is the Hamaker constant and h is the distance between the robot and the wall. Then, the Van der Waals interaction force is given by differenciating [START_REF] Kehlenbeck | Empirical relationships for the terminal settling velocity of spheres un cylindrical columns[END_REF]:

⃗ F vdw = -∇ ⃗ V vdw (8)
The electrostatic force between the robot and the wall considered as an uncharged surface is given by [START_REF] Hays | Electrostatic adhesion of non-uniformly charged dielectric sphere[END_REF]:

⃗ F elec = q 2 4πϵϵ 0 (r + h) 2 ⃗ n (9)
with q the robot charge, ϵ the dielectric density of the medium in which the interaction occurs and ϵ 0 the vacuum permittivity. [START_REF]Role of Electrostatics in Adhesion, in Fundamentals of Adhesion[END_REF] gives the expression of the maximum allowable charge for a spherical body of radius r:

q(µC) = S × Q = 4πr 2 × 30(100r) -0.3

III. OPTIMAL TRAJECTORY AND CONTROL APPROACH

Previous forces balance gives us sufficient informations to plan an optimal trajectory. At least two positions A and B, shown on the Figure 1, should be taken into account. In the first one, the robot is in a vertical vessel and the magnetic force ⃗ F m should counter both contributions of the robot's apparent weight ⃗ W a and the drag force ⃗ F d when blood is flowing back (Curve A of Figure 3). The drag force decreases when the robot approaches the wall due to the parabolic profile of velocity. Thus the reference trajectory should be near the wall.

In the second case, the robot is in a horizontal vessel and the magnetic force should counter contributions of the robot's apparent weight, electrostatic and Van der Waals forces, ⃗ F elec and ⃗ F vdw (Curve B of Figure 3). This case shows that near the wall, ⃗ F elec and ⃗ F vdw , which point to the wall, are dominant and the magnetic force is no more sufficient to counter it. Nevertheless there is an optimal position where the sum of the two forces compensates perfectly the robot's weight. Moreover, the curve A shows that at this point, the magnetic force exceeds the drag force.

From these observations, we define an optimal path as an arc passing through the point C. A backstepping control approach have been developed in order to ensure a null error between real and desired trajectory. The reader may refer to [START_REF] Arcese | Nonlinear modeling and robust controller-observer for a magnetic microrobot in a fluidic environment using mri gradients[END_REF] for details on the control approach.

IV. SIMULATIONS

Simulations are performed by taking into account the limitations of a clinical MRI system. The pulsative blood's velocity is modeled by an affine combination of a timevarying periodic flow with a spatial parabolic shape. So as to simplify the analytical expression, but with no loss of generality, we only consider the first terms in the timevarying Fourier series of the physiological pulse. In the case of an artery, such an approximation leads to:

v f (t) = 0.025(1 + 1.15 sin 2πt) × [ 1 - ( D/2 -h r ) 2 ]
Our studies assume the presence of minor bifurcations (see Figure 1). The developed controller must be sufficiently robust to compensate this effect which could be considered as a disturbance. Major bifurcations will require a further study of velocity's field profile.

In the following, the performances and stability of the controller with respect to noise measurement, parameters variations and uncertainties are illustrated by a simulation, whose parameters are given in Table I. The simulation is performed by assuming that the blood's viscosity and permittivity are affected by uncertainties of 100% of their nominal values and the vessel's diameter by uncertainties of 10%. Besides, a white gaussian noise about 10% of the measured signal is applied on the position measurement.

This simulation shows that after a transient phase (see figure 4), the position tracking performance is robust enough to model's error, but proved quite sensitive to output noise, though remains stable. Despite these uncertainties and disturbances, Figure 5 shows that the control inputs do not reach saturation. This is due to the choice of an optimal reference trajectory deduced from the analysis of Section III.

V. CONCLUSION

In this paper, we have developed a highly nonlinear model for a MRI guided microrobot in blood vessels with minor bifurcations. This model takes into account the nonnewtonian behavior of blood, as well as wall effects and interactions. It makes it possible to hence deduce an optimal trajectory. A backstepping approach have been used to follow desired trajectory. Parameters uncertainties and noise effects have been illustrated by simulation. It appears that the system is robust to uncertain physiological parameters, but proved quite sensitive to output noise, though remains stable.

Actually, we are working on the estimation of the blood's velocity and frequency, assumed to be known in the present simulation. The modeling of the impact of major bifurcations on the blood's velocity profile is also underway.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Scheme of a blood vessel with minor bifurcations

Fig. 3 .

 3 Fig. 3. Forces balances in vertical (A) and horizontal (B) artery

Fig. 4 .Fig. 5 .

 45 Fig. 4. XZ trajectory : reference trajectory (red dotted) and real trajectory (blue solid line)

TABLE I SIMULATIONS

 I

		DATA
	Plasma's viscosity	η plasma	5 × 10 -3 [P a.s]
	Blood's density	ρ f	1060 [kg.m -3 ]
	Robot's density	ρm	8000 [kg.m -3 ]
	Robot's radius	r	300 [µm]
	Vessel's diameter	D	3 [mm]
	Polymer's density	ρ poly	1500 [kg.m -3 ]
	Ferromagnetic ratio	τm	0.8
	Magnetization	M	1.95 × 10 6 [A.m -1 ]
	Hematocrit	h d	0.45
	Robot's Young's modulus	Ep	10 9 [P a]
	Wall's Young's modulus	Ew	0.75 × 10 6 [P a]
	Robot's Poisson's ratio	σp	0.27
	Wall's Poisson's ratio	σw	0.30
	Hamaker constant	A h	4 × 10 -19 [J]
	Blood's dielectric density	ε	77 [C 2 .N -1 .m -2 ]
	Inputs saturations	u i,max	45 [mT.m -1 ]

This work was supported by European Union's 7th Framework Program and its research area ICT-2007.3.6 Micro/nanosystems under the project NANOMA (Nano-Actuactors and Nano-Sensors for Medical Applications).


Laurent Arcese
email: laurent.arcese@bourges.univ-orleans.fr

Ali Cherry
email: ali.cherry@ensi-bourges.fr

Matthieu Fruchard
email: matthieu.fruchard@bourges.univ-orleans.fr

Antoine Ferreira
email: antoine.ferreira@ensi-bourges.fr

Dynamic behavior investigation for trajectory control of a microrobot in blood vessels

This paper reports modeling and control of a microsized polymer aggregate of magnetic particles inside an artery, using a MRI device for supplying propulsion in order to achieve targeted chemotherapy. Non-newtonian behaviour of blood is took into account, as well as wall effects and interactions, resulting in a highly nonlinear model. A backstepping approach is synthesized to ensure Lyapunov stability along a pre-planned trajectory inherited from the model, with robustness concerns.

I. INTRODUCTION

Autonomous microrobots designed to perform targeted therapy by navigating in the cardiovascular system are a prolific research area for minimally invasive surgery purposes [START_REF] Abbott | Robotics in the small: Part 1 microrobotics[END_REF]. Untethered robots are particularly well suited either for reaching remote parts of the body which remain inaccessible to present tools without operating, or for targeting therapy in such places. In the case of chemotherapy, such robots carrying drugs can both improve the therapy's efficiency, reduce drug dose and in turn side effects. On such a scale, it proves difficult to embed actuators sufficiently powerful to propel the robot, especially when swimming against the tide. Therefore, the great majority of swimming approaches considers magnetic fields to wirelessly transmit power to the robot [START_REF] Gillies | Magnetic manipulation instrumentation for medical physics research[END_REF], [START_REF] Honda | Micro swimming mechanisms propelled by external magnetic fields[END_REF].

Bead pulling has first been studied using magnets [START_REF] Quate | Goniometric motion controller for the superconducting coil in amagnetic stereoaxis system[END_REF] or superconducting magnets [START_REF] Takeda | Development of magnetically targeted drug delivery system using superconducting magnet[END_REF]. The last innovation in this domain has been provided by Martel [START_REF] Mathieu | Method of propulsion of a ferromagnetic core in the cardiovascular system through magnetic gradients generated by an mri system[END_REF], where the basic idea is to use the gradients coils of a clinical MRI to pull the robot. Since medical applications are the underlying theme of these works, using a MRI device to actuate the robot combines several advantages. MRI devices are widely implanted and provide magnetic gradients that can be used for propelling purposes. Besides, the MRI imager can provide fine observation of the scene for navigation, using multiplexing for controlling and imaging [START_REF] Tamaz | Real-time mri-based control of a ferromagnetic core for endovascular navigation[END_REF].

What is the optimal size for such a robot? On the one hand -from a therapeutic viewpoint-the smaller the better, so as to avoid embolization far from the tumor and to drive the drugs as close as possible to the tumor. On the other hand, since magnetic forces are volumic, they decrease in cubic with size. Besides, partial vessel occlusion by the robot results in an optimal ratio between robot's and vessel's radii evocated in [START_REF] Mathieu | Method of propulsion of a ferromagnetic core in the cardiovascular system through magnetic gradients generated by an mri system[END_REF]. Along the reference trajectory, from injection point to tumor's area, vessels radius decreases. It is thus necessary to make the radius of the robot decrease as it goes in smaller vasculature. This can be done if the robot is a polymer binded aggregate of magnetic particles, loosing substance with time.

This paper deals with both modeling and controlling such a robot along a reference trajectory. This contribution mainly relies on defining a precise forces balance (Section II), which includes wall effects (parabolic profile of blood flow, pulsatile arterial walls and effect of the ratio of robot's on vessel's radii), wall interactions (Van der Waals, contact and electrostatic forces) and non-newtonian behaviour of blood. The goal is not only to describe the robot's behaviour accurately, but also to make a link with off-line approaches [START_REF] Vartholomeos | Simulation platform for selfassembly structures in mri-based nanorobotic drug delivery systems[END_REF] without their time-consuming drawbacks. An optimal reference trajectory is hence deduced from the nonlinear model (Section III), in order to minimize the control efforts.

A state space representation is then inherited from the forces balance model (Section IV), so as to synthesize a backstepping control law (Section V). To our knowledge, most of feedback controller schemes designed in the scope of microrobotics facing drag are PID approaches [START_REF] Tamaz | Real-time mri-based control of a ferromagnetic core for endovascular navigation[END_REF], [START_REF] Mathieu | In-vivo validation of a propulsion method for untethered medical microrobots using a clinical magnetic resonance imaging system[END_REF]. Authors report instabilities and important oscillations around the equilibrium, especially when the blood stream is modeled as a pulsatile flow. They also indicate a lack of robustness to noise and unmodeled dynamics. In this paper, the control design is based on the backstepping approach recently developed in [START_REF] Arcese | Nonlinear modeling and robust controller-observer for a magnetic microrobot in a fluidic environment using mri gradients[END_REF], which ensures Lyapunov stability. Since backstepping can be regarded as a pre-compensation of the model's nonlinearities, working with a more accurate model is a real improvement of our previous work. Some simulations results illustrate the strengths and limitations of the overall study (Section VI). Lastly we sum up results and discuss prospects (Section VII).

II. MODELING: FORCES BALANCE

The purpose of this section is to present a 2D state space model of a ferromagnetic microrobot of radius r immersed in blood vessel. The model encompasses the different forces that affect the robot's motion as well as its interaction with the vessel wall. The translational and rotational motions of the robot are expressed by:

m d v dt = F m + F d + W a + F c + F vdw + F elec J d w dt = T m + T d + M c (1)
where v and w are, respectively, the translational and rotational velocity of the robot, m and J its mass and moment of inertia. F m , F d , W a , F c , F vdw and F elec respectively denote the magnetic force produced by the MRI gradient coils, blood hydrodynamic drag force, apparent weight, the robot-to-wall contact force, the Van der Waals force and the electrostatic force. T m , T d and M c denote respectively the magnetic torque, hydrodynamic drag torque and the robotto-wall contact moment.

In the rest of this paper, we assume that the orientation of the robot does not change due to the magnetic torque which tends to align the magnetization of the robot along the external field. T d and M c are much smaller order than T m .

To determine whether statistical mechanics or continuum mechanics formulations of fluid dynamics should be used, we refer to the dimensionless Knudsen number K n = k b T P πr 3 2 5/2 , where T , k b and P respectively denote temperature, the Boltzmann constant and pressure. In our case, the low value K n ≈ 10 -13 (≪ 10 -3 ) ensures that the robot is large enough to neglect the effect of Brownian motion: the robot's motion is subjected to generic Navier-Stokes equations.

A. Magnetic force

The gradient coils of the MRI system provide magnetic gradients which produce a magnetic force F m on the robot:

F m = τ m µ 0 V ( M .∇) H (2)
where V is the robot's total volume, τ m = Vm V with V m the ferromagnetic volume, M is the magnetization of the material, µ 0 is the permeability of free space, H is the external magnetic field and ∇ is the gradient operator.

B. Hydrodynamic drag force

The hydrodynamic drag force F d exerting on a spherical body in an infinite extent of fluid is expressed as:

F d = - 1 2 ρ f ( v -v f ) 2 ∞ AC d∞ (3)
where v -v f denotes the relative velocity of the robot with respect to the fluid, A is the frontal area of the core and ρ f is the density of the fluid. C d∞ is the drag coefficient, a dimensionless quantity used to quantify the drag or resistance of an object in a fluid.

C d∞ is a function of the Reynolds number Re ∞ = 2rρ f | v-v f | ∞ η
; its expression for a spherical body in laminar flow is given by [START_REF] White | Viscous Fluid Flow[END_REF]:

C d∞ = 24 Re ∞ + 6 1 + √ Re ∞ + 0. 4
1) Non-Newtonian fluid: In case of blood, which exhibits a non-newtonian behavior, the fluid's viscosity η is a function of vessel diameter d (in micron) and hematocrit rate h d according to the following empirical relations [START_REF] Pries | Biophysical aspects of blood flow in the microvasculature[END_REF]:

η = η plasma d 2 (d -1.1) 2 1 + (η 0.45 -1)d 2 (d -1.1) 2 (1 -h d ) c -1 (1 -0.45) c -1 B A Output Blood Profile Input Blood Profile W a F c F m W a F e + F vdw F d F m i θ k C r e f
Fig. 1. Scheme of a blood vessel with minor bifurcations with parameters η plasma and η 0.45 denoting respectively the plasma's viscosity and the relative apparent blood viscosity for a fixed discharge hematocrit of 0.45, given by: η 0.45 = 6 e -0.085d + 3.2 -2.44 e -0.06d 0.645

The shape of the viscosity dependance on hematocrit is: 2) Wall effects: For endovascular applications, influence of the vessel walls on the velocity of the robot has to be taken into account. In general, this wall effect is expressed as a ratio between the terminal relative velocity of the robot ( v -v f ) t and its velocity ( v -v f ) ∞ in an infinite extent of fluid [START_REF] Kehlenbeck | Empirical relationships for the terminal settling velocity of spheres un cylindrical columns[END_REF]:

c = 10
| v -v f | t | v -v f | ∞ = 1 -λ α0 1 + λ λ0 α0 (4)
with ratio λ = 2r/D and D denoting the vessel diameter (in meter). Parameters α 0 and λ 0 are functions of Reynolds number, but are commonly set to 1.5 and 0.29, respectively. Thus, equation ( 3) is corrected as follows:

F d = - 1 2 ρ f     1 + λ λ0 α0 1 -λ α0   ( v -v f ) t   2 AC dt (5)
Wall effects on the fluid in the vessel traditionally result in a parabolic profile of blood flow (see Figure 1). Besides, to fully take into account pulsatile flow caused by heart pumping in arteries, one has to consider a periodic deformation of the vessel's diameter D(t) synchronized with v f (t).

C. Apparent weight

In addition to the magnetic and hydrodynamic forces, apparent weight (combined action of weight and buoyancy) is acting on the spherical robot:

W a = V (ρ -ρ f ) g (6)
where ρ = τ m ρ m + (1 -τ m )ρ poly with ρ m and ρ poly the magnetic material's and polymer's densities.

D. Contact force

The normal and tangential interaction between the robot and the wall are assumed to be expressed by a Voigt model with the spring constant K and the decay coefficient of dashpot b, as shown in Figure 2, where the indexes n and t indicate respectively the normal and tangential normalized vectors.

The normal component of the contact force F cn acting on the robot is given by the sum of a non-linear elastic component and a linear viscous damping: [START_REF] Tijskens | Discret element modelling for process simulation in agriculture[END_REF]:

F cn = (K n δ 3/2 + b n δ) n ( 7
)
where δ is the elastic deformation of the wall at the contact point and n is the normal unit vector pointing from the robot to the contact surface. The stiffness K n can be calculated by Hertzian contact theory [START_REF] Tsuji | Lagragian numerical simulation of plug flow of cohesionless particles in a horizontal pipe[END_REF] as follows:

K n = 4 √ r 3 1-σ 2 p Ep + 1-σ 2 w
Ew where E p and E w are the Young's modulus of robot and wall, and σ p , σ w are respectively the Poisson's ratios of robot and wall. The damping coefficient b n is deduced from K n using b n = 2 √ mK n . The tangential component of the contact force takes part when the robot is rotating or in case of oblique collision with the wall:

F ct = (K t ζ + b t ζ) t ( 8
)
where ζ is the displacement in the tangential direction and t is the tangential unit vector. The stiffness K t and damping coefficient b t are given by:

K t = 8 √ rG p 2 -σ p δ 1/2 , b t = 2 mK t
where G p is the shear modulus related to the Young's modulus E p and Poisson ratio σ p , i.e. G p = Ep 2(1+σp) .

E. Van der Waals and electrostatic forces

When the robot and the wall are not in contact, they interact each other through Van der Waals and electrostatic forces. These two interaction forces have different dominant regimes. In fact, when the robot is close to the wall, Van der Waals force is dominant. As the robot move away from the wall, the Van der Waals force rapidly decrease and the electrostatic force becomes dominant.

The Van der Waals potential between the robot and the wall is given by [START_REF] Iimura | Simulation of entrainment of agglomerates from plate surfaces by shear flows[END_REF]:

V vdw = - A h 6 1 h + 1 2 + h + ln h 2 + h n (9)
where A h is the Hamaker constant and h is the distance between robot and wall. Then, the Van der Waals interaction force is given by differenciating ( 9) F vdw = -∇ V vdw . The electrostatic force between the robot and the wall considered as an uncharged surface is given by [START_REF] Hays | Electrostatic adhesion of non-uniformly charged dielectric sphere[END_REF]:

F elec = q 2 4πǫǫ 0 (r + h) 2 n ( 10
)
with q the robot charge, ǫ the dielectric density of the medium in which the interaction occurs and ǫ 0 the vacuum permittivity. [START_REF]Role of Electrostatics in Adhesion, in Fundamentals of Adhesion[END_REF] gives the expression of the maximum allowable charge Q = 30 × (100r) -0.3 for a spherical body of radius r.

III. MODELING: OPTIMAL TRAJECTORY

Previous forces balance gives us sufficient informations to plan an optimal trajectory. At least two cases A and B, shown on the Figure 1, should be taken into account. In the first one, the robot is in a vertical vessel and the magnetic force F m should counter both contributions of the robot's apparent weight W a and the drag force F d when blood is flowing back (Curve A of Figure 3). The drag force decreases when the robot approaches the wall due to the parabolic profile of velocity. Thus the reference trajectory should be as near as possible to wall. In the second case, the robot is in a horizontal vessel and the magnetic force should counter contributions of the robot's apparent weight, electrostatic and Van der Waals forces, F elec and F vdw (Curve B of Figure 3). This case shows that near the wall, F elec and F vdw , which point to the wall, are dominant and the magnetic force is no more sufficient to counter it. Nevertheless there is an optimal position where the sum of the two forces compensates perfectly the robot's weight. Moreover, the curve A shows that at this point, the magnetic force exceeds the drag force. From these observations, we define an optimal path as an arc passing through the point C.

IV. MODELING: STATE SPACE REPRESENTATION

Let (x, z) denote the position of the robot in the vessel with respect to a given frame F (O, i, k). The state model is established from differential equation (1) defining the robot's dynamic behavior, projected on i and k axes: where indexes x (respectively z) denote projections on i ( k).

mẍ = F mx + F dx + F cnx + F ctx + F vdwx mz = F mz + F dz + W a + F cnz + F ctz + F vdwz ( 11
)
Let x 1 , x 2 , x 3 , x 4 denote respectively the particule's position and velocity along i axis, and the position and velocity along k axis. Assuming that positions x 1 and x 3 can be measured thanks to the MRI imaging system, let y denote the state's measure. Using expressions of forces given by ( 2), ( 5), ( 6), ( 7), ( 8), [START_REF] Kehlenbeck | Empirical relationships for the terminal settling velocity of spheres un cylindrical columns[END_REF] and Van der Waals interaction derivated from potential (9), and adequate projection of local frame along the geometry of the vessel 1 , system (11) can be written in the form:

(S)                ẋ1 = x 2 ẋ2 = f 2 (x 1 , x 2 , x 3 , x 4 ) + au 1 ẋ3 = x 4 ẋ4 = f 4 (x 1 , x 2 , x 3 , x 4 ) + au 2 y = x 1 x 3 (12)
where control inputs u 1 = ∇B x and u 2 = ∇B z are the magnetic gradients, parameter a = τmM ρ , and f i are very highly nonlinear functions of the full state.

V. CONTROL APPROACH

In this section, we present a control design for Lyapunovstabilizing trajectories for system (S). Since determination of Lyapunov functions is generically a challenging issue, it is preferable to use Lyapunov control functions in a backstepping control approach [START_REF] Krstić | Nonlinear and Adaptative Control Design[END_REF], [START_REF] Kanellakopoulos | A toolkit for nonlinear feedback design[END_REF]. Since this design requires a triangular form for the control system, we propose the following change of coordinates:

X = x 1 x 3 ; Y = x 2 x 4 ; U = u 1 u 2 (13)
Thus, we obtain a new system in a triangular form:

(S ′ )    Ẋ = Y Ẏ = F (X, Y ) + U y = X (14)
1 The drag force is not linearly distributed between i and k axis. Denoting X ref , Ẋref and Ẍref respectively the desired reference trajectory, velocity and acceleration, the control law for system (S ′ ) can be expressed as [START_REF] Arcese | Nonlinear modeling and robust controller-observer for a magnetic microrobot in a fluidic environment using mri gradients[END_REF]:

U = Ẍref + (k 1 + k 2 ) Ẋref + (1 + k 1 k 2 )X ref -(1 + k 1 k 2 )X -(k 1 + k 2 )Y -F (X, Y ) (15)
where k 1 and k 2 are the strictly positive backstepping gains.

VI. SIMULATIONS

The robot is made of τ m = 80% NdFeB particles which has a combination of very high remanence and coercivity, and 20% of binding polymer. Simulations are performed by taking into account the limitations of a clinical MRI system. In order not to exceed the capacity of existing MRI systems, the applied control law ( 15) is now corrected as ui k(t) , with k(t) = max 1, ui ui,max .

A pulsatile flow is included by imposing a time-varying velocity. As a first approximation of a physiological pulse, we use a time-sinusoidal profile with spatial parabolic form. In the case of artery, such an approximation leads to:

v f (t) = 0.035(1 + 1.15 sin 2πt) × 1 - D/2 -h r 2
Our studies assume the presence of minor bifurcations (see Figure 1). This geometry leads to only slight change in the velocity's profile and amplitude (see Figure 4). The developed controller must be sufficiently robust to compensate this effect which could be considered as a disturbance. Major bifurcations will require a further study of velocity's field profile.

The purpose of the control is to ensure a null error between desired and a pre-planned trajectory. In the following, the performances and stability of the controller with respect to I. A first simulation is led so as to pinpoint robustness to model errors, and the second one to study sensitivity to output noise.

Simulation of Figures 5,6 is performed by assuming that the blood's viscosity and permittivity are affected by uncertainties of 100% of their nominal values and the vessel's diameter by uncertainties of 10%. This simulation shows that after a transient phase, the position tracking performance is robust enough to model's error (subfigure of Figure 5). Despite these uncertainties, Figure 6 shows that the control inputs don't reach saturation. This is due to the choice of an optimal reference trajectory deduced from the analysis of Section III.

In a second simulation (Figures 7 and8), we assume a white gaussian noise is applied on the position measurement. This noise is about 10% of the measured signal. Figure 7 illustrate that the controller is quite stable despite the noise, and does not even reach saturations (Figure 8). Nevertheless, it's clear that the system is sensitive to noise. We can expect that this drawback effect can be reduced using an observer. In addition to estimate some parameters like blood's velocity (assumed to be known in this paper), an observer will tend to smooth the output signals.

VII. CONCLUSION

In this paper, we have presented a highly nonlinear model for a MRI guided microrobot in blood vessels. This model takes into account the non-newtonian behavior of blood, as well as wall effects and interactions. It makes it possible to hence deduce an optimal trajectory. Besides, we have developed a nonlinear control law based on the backstepping approach. Parameters uncertainties and noise effects have been illustrated by simulations. It appears that the system is robust to uncertain physiological parameters, but proved quite sensitive to output noise, though remains stable.

As stated in Section VI, system robustness to noise can be improved by implementing a high gain observer based on [START_REF] Arcese | Nonlinear modeling and robust controller-observer for a magnetic microrobot in a fluidic environment using mri gradients[END_REF] but fitted for this more complete and complex model. We are also working on estimating the blood's velocity and frequency, assumed to be known at the moment, using Kazantis-Kravaris Luenberger observers [START_REF] Andrieu | On the existence of a kazantis-kravaris / luenberger observer[END_REF]. Additionnaly, the modeling of impact of major bifurcations on the blood's velocity profile is underway.

Fig. 2 .

 2 Fig. 2. Model of contact forces: robot-to-wall

Fig. 3 .

 3 Fig. 3. Forces balances in vertical (A) and horizontal (B) artery

Fig. 4 .

 4 Fig. 4. Comsol results : Velocity profile in a blood vessel with minor bifurcations

Fig. 5 .

 5 Fig. 5. XZ trajectory : reference trajectory (red dotted) and real trajectory (blue solid line)

Fig. 6 .

 6 Fig. 6. Control input : magnetic field gradients on i-axis (dotted) and on k-axis (solid line)

Fig. 7 .Fig. 8 .

 78 Fig. 7. XZ trajectory : reference trajectory (red dotted) and real trajectory (blue solid line)

TABLE I SIMULATIONS

 I

	Robot's Poisson's ratio	σp	0.27
	Wall's Poisson's ratio	σw	0.30
	Hamaker constant A Controller gains k 1	15
		k 2	30

DATA Plasma's viscosity η plasma 5 × 10 -3 [P a.s] Blood's density ρ f 1060 [kg.m -3 ] Robot's density ρm 8000 [kg.m -3 ] Robot's radius r 300 [µm] Vessel's diameter D 3 [mm] Polymer's density ρ poly 1500 [kg.m -3 ] Ferromagnetic ratio τm 0.8 Magnetization M 1.95 × 10 6 [A.m -1 ] Hematocrit h d 0.45 Robot's Young's modulus Ep 10 9 [P a] Wall's Young's modulus Ew 0.75 × 10 6 [P a] h 4 × 10 -19 [J] Blood's dielectric density ε 77 [C 2 .N -1 .m -2 ] Initial condition on x x 0 (0, 0, 0, 0) T Inputs saturations u i,max 45 [mT.m -1 ]

VIII. ACKNOWLEDGMENTS

This work was supported by European Union's 7th Framework Program and its research area ICT-2007.3.6 Micro/nanosystems under the project NANOMA (Nano-Actuactors and Nano-Sensors for Medical Applications). The authors thank Christophe Boulnois for Comsol helpful advice.


Francis Comets
email: comets@math.jussieu.fr

François Simenhaus
email: simenhaus@math.jussieu.fr

Random walk delayed on percolation clusters

Keywords: Random walk on percolation clusters Random walk in random environment, subcritical percolation, anomalous transport, anomalous diffusion, environment seen from the particle, coupling AMS 2000 subject classifications: 60K37; 60F15, 60K35, 82D30

We study a continuous time random walk on the d-dimensional lattice, subject to a drift and an attraction to large clusters of a subcritical Bernoulli site percolation. We find two distinct regimes: a ballistic one, and a subballistic one taking place when the attraction is strong enough. We identify the speed in the former case, and the algebraic rate of escape in the latter case. Finally, we discuss the diffusive behavior in the case of zero drift and weak attraction.

clusters size has an exponential tail: For any p < p c , there exists ξ = ξ(p) > 0 such that for all x, lim n→∞ 1 n ln P(C x n) = -ξ .

In this paper, we fix p < p c . Let ℓ = (ℓ k ; 1 k d) be a unit vector, λ and β two non-negative number. For every environment ω, let P ω be the law of the continuous time Markov chain Y = (Y t ) t 0 on Z d starting at 0 with generator L given for continuous bounded functions f by

Lf (x) = K e∼0 e λℓ•e-βCx f (x + e) -f (x) ,
where we chose the normalizing constant K as K = e∼0 e λℓ•e -1 for simplicity. Given ω, define the measure µ on Z d by µ(x) = e 2λℓ•x+βCx .

(

) 1
The random measure µ combines a shift in the direction ℓ together with an attraction to large clusters. Observe that the process Y admits µ as invariant, reversible measure. Markov processes having µ as invariant measure are of natural interest in the context of random walks in random environment. They describe random walks which have a tendency to live on large clusters, the attraction becoming stronger as β is increased. The isotropic case, λ = 0, has been considered in [START_REF] Popov | Random walk attracted by percolation clusters[END_REF] with a different, discrete-time dynamics. There, the authors proved that the walk is diffusive for small β, and subdiffusive for large β.

The investigation of slowdowns in the anisotropic case is then natural. In [START_REF] Shen | Asymptotic properties of certain anisotropic walks in random media[END_REF], a random resistor network is considered with a invariant reversible measure of the form C(x, ω)e 2λℓ•x where the random field (C(x, ω); x ∈ Z d ) is stationary ergodic and bounded away from 0 and +∞: in this case, the random walks in random environment is ballistic for all positive λ.

The study of a general dynamics in the presence of a drift contains many difficult questions, and the advantage of the particular process Y considered here is that we can push the analysis farther. We could as well handle the discrete time analogous of Y , i.e. the random walks in random environment with geometric holding times instead of exponential ones, which falls in the class of marginally nestling walks in the standard classification (e.g., [START_REF] Zeitouni | Random walks in random environment[END_REF]). The Markov process Y can also be described with it skeleton and its jump rates. The skeleton X = (X n ) n∈N is defined as the sequence of distinct consecutive locations visited by Y . Then, X is a discrete time Markov chain with transition probabilities P , given for x ∈ Z d and e ∼ 0 by

∀x ∈ Z d , ∀e ∼ 0, P (X n+1 = x + e|X n = x) = e λℓ•e
e ′ ∼0 e λℓ•e ′ =: p e , and P (X n+1 = y|X n = x) = 0 if y is not a nearest neighbor of x. The jump rate of Y at site x is equal to exp -βC x , and the holding times are independent, exponentially distributed with mean exp βC x . The Markov chain X is quite simple, it is the random walk on Z d with drift

d(λ) = 1 d k=1 cosh(λℓ k ) sinh(λℓ k ) 1 k d . (2)
It is plain that for the random walk,

X n n -→ d(λ) P -a.s., (3)
so directional transience is clear, and the law of large number for Y boils down to studying the clock process which takes care of the jump times. As can be seen from formula [START_REF] Durrett | Multidimensional random walks in random environments with subclassical limiting behavior[END_REF], the process considered here is a generalization of the so-called random walk in a random scenery, or the random walk subordinated to a renewal process, which are used as effective models for anomalous diffusions. The difference is essentially that the environment (i.e., the field of jump rates) has here some space correlations, which are short-range. It is also related to the trap model considered in the analysis of the aging phenomenon introduced in [START_REF] Bouchaud | Weak ergodicity breaking and aging in disordered systems[END_REF]: the aging of this model has been studied in details, see [START_REF] Ben Arous | Dynamics of Trap Models[END_REF] for a recent review. For a fixed ω, P ω is called the quenched law and we define the annealed law P by

P = P × P ω .
Of course, statements which hold P -a.s., equivalently hold P ω -a.s. for P-a.e. environment. Finally, we stress that we assume d 1 in this paper. The case d = 1 is special since the critical threshold p c (1) = 1. Moreover, specific techniques are available in one dimension, e.g. [START_REF] Zeitouni | Random walks in random environment[END_REF] for a survey, however we will stick as much as possible to techniques applying for all d.

Our first result is the law of large numbers.

Theorem 1. (Law of large numbers) For any λ 0 and any β 0,

Y t t ----→ t→+∞ v(λ, β), P -a.s., where v(λ, β) = Ee βC 0 -1 d(λ) . (4)
In particular, v(λ, β) = 0 if β > ξ or λ = 0 though v(λ, β) • ℓ > 0 if β < ξ and λ = 0.

As in the case λ = 0 considered in [START_REF] Popov | Random walk attracted by percolation clusters[END_REF], slowdowns occur for large disorder intensity β, when the walk gets trapped on large percolation clusters. This behavior is reminiscent of the biased random walk on the supercritical percolation infinite cluster [START_REF] Sznitman | On the anisotropic walk on the supercritical percolation cluster[END_REF], [START_REF] Berger | The speed of biased random walk on percolation clusters[END_REF] where ballistic or subballistic regimes take place according to the parameters values. The slowdowns in our paper have a similar nature to those in some one dimensional random walks in random environment, see [START_REF] Solomon | Random walks in a random environment[END_REF], [START_REF] Kesten | A limit law for random walk in a random environment[END_REF] and [START_REF] Sinai | The limit behavior of a one-dimensional random walk in a random environment[END_REF]. Moreover, as in the one dimensional case, we obtain here explicit values for the rate of escape, a rather unusual fact in larger dimension. More drastic (logarithmic) slowdowns were also found for an unbiased walker in a moon craters landscape in [START_REF] Durrett | Multidimensional random walks in random environments with subclassical limiting behavior[END_REF], [START_REF] Bramson | Random walk in random environment: a counterexample?[END_REF], or diffusions in random potentials [START_REF] Mathieu | Zero white noise limit through Dirichlet forms, with application to diffusions in a random medium[END_REF], but in these models the behavior at small disorder is qualitatively different from the behavior without disorder.

The next result contains extra information on the subballistic behavior.

β 2ξ + 1 2 -1 P -a.s.
Hence, the spread of the random walks in random environment scales algebraically with time in all cases. Note that in the isotropic case λ = 0, the slowdown is larger for d = 1 than for d 2. This will appear in the proof as a consequence of the strong recurrence of the simple random walk X in the one-dimensional case. Note that our results are only in the logarithmic scale, though the scaling limit has been obtained for the isotropic trap model, in dimension d = 1 (e.g., [START_REF] Ben Arous | Dynamics of Trap Models[END_REF]), and d 2 [START_REF] Ben Arous | Scaling limit for trap models on Z d[END_REF] with limit given, if the disorder is strong, by the time change of a Brownian motion by the inverse of a stable subordinator (fractional kinetics). Though we believe that the scaling limit of our model without drift (λ = 0) is the same, we could not get finer results because of the presence of correlations in the medium. Moreover, the case of a drift λ = 0 has not been considered in the literature, except for d = 1 with renormalization group arguments [START_REF] Monthus | Nonlinear response of the trap model in the aging regime: Exact results in the strong-disorder limit[END_REF].

To complete the picture, we end by the diffusive case. (Recall that β < ξ is sufficient for E(e βC 0 ) < ∞.) Theorem 3. (Diffusive case regime) Assume λ = 0, and E(e βC 0 ) < ∞. Then, we have a quenched invariance principle for the rescaled process Z ǫ = (Z ǫ t ) t 0 , Z ǫ t = ǫ 1/2 Y ǫ -1 t : For almost every ω, as ǫ ց 0, the family of processes Z ǫ converges in law under P ω in the Skorohod topology to the d-dimensional Brownian motion with diffusion matrix

Σ = d × E(e βC 0 ) -1 I d . Moreover, lim sup t→+∞ ln |Y t | ln t = 1 2 a.s. (5)
For the proof of our results we will take the point of view of the environment seen from the walker. It turns out that the "static" environmental distribution is invariant for the dynamics. Hence the environment is always at equilibrium.

The paper is organized as follows. In the next section, we introduce the basic ingredients for our analysis and we prove the law of large numbers of Theorem 1. The last section is devoted to the subballistic regime and contains the proofs of Theorem 2 and 3.

Preliminaries and the proof of Theorem 1

For x ∈ Z d , T x will denote the space shift with vector x. We will consider also the time shift θ.

Skeleton and clock process of Y . The sequence (S n ; n 0) of jump times of the Markov process Y with right-continuous paths is defined by

S 0 = 0 < S 1 < S 2 < . . ., Y t = Y Sn for t ∈ [S n , S n+1 ), Y S n+1 = Y Sn .
The skeleton of Y is the sequence X given by X n = Y Sn , n 0. As mentioned above, the skeleton X of Y is the simple random walk with drift. For any x in Z d , the jump rate of (Y t ) t 0 at x is e -βCx . Hence the time S n of the n-th jump is the sum of n independent random variables with exponential distribution with mean e βC X i , i = 1, . . . n. This means that the sequence E = (E i ) i∈N , with E i = e -βC X i (S i+1 -S i ), is, under the quenched law and then also under the annealed law, a sequence of i.i.d. exponential variables with mean 1, with E and X independent. The law of this sequence will be denoted by Q (Q = Exp(1) ⊗N , with Exp(1) the mean 1, exponential law). For any n in N, the time S n of the n-th jump is given by

S n = n-1 i=0 E i e βC X i . ( 6
)
This sequence can be view as a step function

S t := S [t]
, where [•] is the integer part, and we also define its generalized inverse S -1 : for any t 0,

S -1 (t) = n ⇐⇒ S n t < S n+1 .
We observe that S n → ∞ as n → ∞ P ω -a.s. for all ω, making the function S -1 defined on the whole of R + . Then, P ω -a.s.,

X S -1 (t) = Y (t) , ∀t 0 . (7)
and therefore, the process S -1 is called the clock process. Conversely, let E, X and ω be independent, with distribution Q, P and P respectively, defined on some new probability space. Then, fixing λ and viewing β as a parameter, by ( 6) and [START_REF] Bramson | Random walk in random environment: a counterexample?[END_REF] we construct, on this new probability space, a coupling of the processes Y = Y (β) for all β ∈ R. The coupling has the properties that the skeleton is the same for all β, and that the clock processes are such that for β β ′ and t 0,

S -1 (β; t) S -1 (β ′ ; t). ( 8
)
The environment seen from the walker. Depending on the time being discrete or continuous, we consider the processes ( ω n ) n∈N and (ω t ) t 0 defined by

ω n = T Xn ω , ωt = T Yt ω = ω S -1 (t)
for n 0, t 0. We start with the case of discrete time.

Lemma 1. Under P , ( ω i ) i∈N is a stationary ergodic Markov chain. The same holds for

( ω i , E i ) i∈N .
Proof of Lemma 1. As (E i ) i∈N is an i.i.d. sequence of variables independent of ω, it is enough to prove Lemma 1 for the process ( ω i ) i∈N . Under P (resp P ω ) ( ω i ) i∈N is markovian with transition kernel R defined for any bounded function f by

Rf (ω) = e∼0 p e f (T e ω) ∀ω ∈ Ω,
and initial distribution P (resp δ ω ). The transitions of ( ω i ) i∈N does not depend on ω like those of X and, in this sense, the sequence is itself a random walk. Since P is invariant by translation,

E[f ( ω 1 )] = e∼0 p e f (T e ω)dP = e∼0 p e f (T e ω)dP = E[f (ω)],
showing that P is an invariant measure for ( ω i ) i∈N . We will use F to denote the product σ-field on Ω N , and for any k 0, F k will denote the σ-field generated by the k first coordinates. Note that θ is measurable and preserves the law of ω under P . We have to prove that the invariant σ-field Σ := {A ∈ F, 1 A ( ω) = 1 A (θ ω), P -a.s.} is trivial. Let Y be a Σ-measurable bounded random variable on Ω N , we have to show that it is P -a.s. constant.

Define for all ω in Ω, h Y (ω

) := E ω [Y ]
. We will study this function with standard arguments e.g. chapter 17.1.1 of [START_REF] Meyn | Markov chains and stochastic stability[END_REF]. Using Markov property and the θ-invariance of Y , we can show that,

h Y ( ω k ) = E[Y |F k ] ∀k ∈ N, P -a.s. ( 9
)
As a consequence, under P , (h Y ( ω k )) k 0 is both a stationary process and an a.s. convergent martingale, and hence it is a.s. constant. In particular,

Y = h Y ( ω 0 ) P -a.s.,
what means that Y can be consider as a function of the first coordinate alone. The next step is to show that h Y is P-a.s. harmonic, that is

Rh Y ( ω 0 ) = h Y ( ω 0 ), P -a.s.
It is a consequence of the following computation,

Rh Y ( ω 0 ) = E[h Y ( ω 1 )|F 0 ] P -a.s. = E[E[Y |F 1 ]|F 0 ] P -a.s. = h Y ( ω 0 ) P -a.s.,
where the second equality is true because of (9). We will now show that Y is invariant by translation in space. By invariance of P and harmonicity of h Y , it is true that

e∼0 p e (Y -Y • T e ) 2 dP = 0.
For every e neighbour of 0, p e > 0, and the previous equation implies that, P almost surely Y = Y • T e for any e ∼ 0. Together with the ergodicity of this shows that Y is P -a.s. constant, and completes the proof.

As a consequence of Lemma 1 and Birkhoff's ergodic theorem, for any function f in L 1 (Ω N ) (or f non negative),

1 n n-1 k=0 f (θ k ω) n→+∞ -----→ E[f ] P -a.s.
Now, we turn to the time continuous case, and we consider the empirical distribution 1 t t 0 δ ωs ds of the environment seen from the walker up to time t. Our next result is a law of large numbers for this random probability measure. For small β, the empirical distribution converges to some limit P 0 , which is then an invariant measure for (ω t ) t 0 .

Corollary 1. If β < ξ then P -almost surely, the empirical distribution of the environment seen from the walker, 1 t t 0 δ ωs ds, converges weakly to P 0 defined by dP 0 = e βC E[e βC ] dP.

Proof of Corollary 1. We need to show that t -1 0 f (ω s )ds → f dP 0 as t → ∞, for all real bounded continuous function f on Ω. Since e βC 0 is integrable when β < ξ, this follows from the convergence along the sequence t = S n , n → ∞. By [START_REF] Durrett | Multidimensional random walks in random environments with subclassical limiting behavior[END_REF], this is equivalent to

n -1 n-1 i=0 E i e βC X i f ( ω i ) n -1 n-1 i=0 E i e βC X i -→ Ω f dP 0 , n → ∞.
We first study the P -almost sure convergence of the denominator, i.e. of n -1 S n . Define the real function g on (R N , Ω N )

g : ((E i ) i∈N , ( ω i ) i∈N ) → E 0 e βC 0 (e ω 0 )
and note that C Xn = C 0 ( ω n ). Applying Lemma 1 and the ergodic theorem to ( ω, E) and to the non negative function g, we obtain that n -1 S n converges P -almost surely to E[e βC 0 ]. The numerator can be studied with the same arguments, and we obtain the claim since for β < ξ both limits are finite.

With this in hand, we can easily complete the Proof of Theorem 1. Write

Y t t = X S -1 (t) S -1 (t))
S -1 (t) S(S -1 (t))

S(S -1 (t))

t .

Recall from (3) that the first factor in the right-hand side converge almost surely to d(λ) as t → ∞. In the proof of Corollary 1 we have shown that S(S -1 (t))/S -1 (t) → E[e βC 0 ] a.s. for β < ξ, but clearly the result remains true for all β (the limit is infinite for β > ξ).

For the last factor in the right-hand side we simply observe that

S(S -1 (t)) S(S -1 (t) + 1) S(S -1 (t)) t 1 , (10)
yielding that S(S -1 (t))/t converges P -almost surely to 1 if E[e βC 0 ] < ∞: in this case, we then conclude that Y t /t converges P -almost surely to v(λ, β) given by ( 4).

In the case E[e βC 0 ] = ∞, we just use the right inequality in [START_REF] Menshikov | Coincidence of critical points in percolation problems. (Russian)[END_REF] to obtain the P -almost surely convergence of Y t /t to v(λ, β) = 0.

Subballistic regime, and the proofs of Theorem 2 and 3

We start with a few auxiliary results. Lemma 2. Assume d 2 or λ > 0. Then, for any ǫ > 0, there exists α > 0 such that P -almost surely, we eventually have

♯ i n, C X i > ( 1 ξ -ǫ) ln n n α
with the notation ♯A for the cardinality of a set A.

Proof of Lemma 2. Define the range R n as the number of points by (X i ) i∈N during the first n steps. For λ > 0, there exists a constant c 1 > 0 such that P -almost surely eventually R n > c 1 n. For λ = 0 and d 2, it is well known (see chapter 21 of [START_REF]Revész Random walk in random and non-random environments[END_REF]) that there exists a constant c 2 such that P -almost surely eventually R n > c 2 n ln n (when d 3, the walk is transient and the correct order of R n is n). In all cases, there exists a constant c 3 > 0 such that under the assumptions of Lemma 2, we have P -almost surely, eventually, R n > c 3 n ln n . For a fixed n in N, we define recursively the time T n i by T n 0 = 0,

T n i = inf{T n i-1 < k n, |X k -X T n j | > 2( 1 ξ -ǫ) ln n, ∀j < i} ∀i 1, inf ∅ = +∞.
Note that the balls with center X T n j and radius (ξ -1ǫ) ln n are pairwise disjoint, and define K n the number of such balls, i.e.

K n = max{i 0 : T n i < +∞} As the cardinality of those ball is c 4 ln d n (for some c 4 > 0), it follows from the previous discussion on the range that P -almost surely, eventually, K n > c n ln d+1 n , where c denotes a positive constant. From now on we fix a path (X i ) i 0 such that K n > c n ln d+1 n eventually. In the rest of the proof, we take n large enough so that the inequality holds. Then,

P ♯ i K n , C X T n i ( 1 ξ -ǫ) ln n K n -n α = P ∃I ⊂ {1, . . . K n }, ♯I = K n -[n α ] : ∀i ∈ I, C X T n i ( 1 ξ -ǫ) ln n I⊂{1,...Kn},♯I=Kn-[n α ] P ∀i ∈ I, C X T n i ( 1 ξ -ǫ) ln n
For all j such that 0 j K nn α , B n i denotes the ball with center X T n i and radius

( 1 ξ -ǫ) ln n. The event {C X T n i ( 1 ξ -ǫ) ln n)} is σ{ω x , x ∈ B n i } measurable.
As the balls B n i are disjoint and the environment is i.i.d.,

P ♯ i K n , C X T n i ( 1 ξ -ǫ) ln n K n -n α K n [n α ] 1 -P(C 0 > ( 1 ξ -ǫ) ln n) Kn-[n α ] c 5 n n α 1 -n -(1-ǫξ)+o(1) c n ln d+1 n -n α ,
for some suitable constant c 5 > 0. We now choose α < min(1, ǫξ), so that

n P ♯ i K n , C X T i ( 1 ξ -ǫ) ln n K n -n α < ∞
We conclude using Borel-Cantelli's lemma. Proof of Lemma 3. Let η be a positive real number. With := η/β, from Lemma 2, there exists α > 0 such that P ⊗ P-almost surely, there exists a natural number N = N (X, ω) such that for n > N , the set

I = {i n, C X i > ( 1 ξ -ǫ) ln n} has cardinality ♯I n α . For n > N , Q(S n < n β/ξ-η ) Q(E i e βC X i < n β/ξ-η , i ∈ I) Q(E 1 e βC X 1 < n β/ξ-η ) n α Q(E 1 < n βǫ-η ) n α = (1 -e -1 ) n α .
From previous inequality, we obtain that Q(S n < n β/ξ-η ) is the general term of a convergent series and we can use Borel-Cantelli's Lemma to conclude. Proof of Lemma 4. For any α in (0, 1), by subadditivity we have (u + v) α u α + v α for all positive u, v, and then

S α n n i=1 E α i e αβC X i . Now, define the function f α f α : (R N , Ω N ) → R ((E i ) i∈N , ( ω i ) i∈N ) → E α 0 e αβC 0 (e ω 0 ) .
Applying Lemma 1 and the ergodic theorem to ( ω, E) with the non negative function f α , we obtain that for any α such that αβ < ξ,

lim sup n→+∞ S α n n lim n→+∞ n i=1 E α i e αβC X i n = E Q (E α 1 ) × E(e αβC 0 ) < ∞ almost surely. Therefore, lim sup n→+∞ ln S n ln n < 1 α .
Since α is arbitrary in (0, ξ/β), the proof is complete.

The two following lemmas deal with the one dimensional case. Notice that when d = 1, for all n > 0,

P(C n) = p n-1 k=0 p k p n-1-k = np n ,
and as a consequence ξ =ln p. Proof of Lemma 5. Here we need to relabel our sequence of exponential variables (E i ; i 0). For y ∈ Z, k ∈ N, define E y,k by E y,k = E i with i such that X i = y, ♯{j : 0 j i, X j = y} = k , i.e. the exponential corresponding to the k-th passage at y. These new variables are a.s. well defined when d = 1 and λ = 0, and it is not difficult to see that the sequence (E y,k ) y∈Z,k∈N is i.i.d. with mean 1 exponential distribution, and independent of X and of ω. The number of visits of the walk to a site y at time n will be denoted by θ(n, y). We can rewrite S n in the following way,

S n = n-1 i=0 e βC X i E i = y∈Z e βCy   θ(n,y)-1 k=0 E y,k   . (11)
Notice that for any η > 0, Pa.s. for n large enough, θ(n, y) = 0 ∀y > n

1 2 +η (12)
(see for example Theorem 5.7 p.44 in [START_REF]Revész Random walk in random and non-random environments[END_REF]). As a consequence, we obtain that for any positive α < 1, P -almost surely for n large enough,

S α n n 1 2 +η y=-n -1 2 +η e αβCy   θ(n,y)-1 k=0 E y,k   α .
Here and below, the sum b y=a with real numbers a < b, ranges over all y ∈ Z with a y b. Notice now that for any ν > 0, Pa.s. for n large enough, sup{θ(n, y), y ∈ Z} < n 1 2 +ν [START_REF] Monthus | Nonlinear response of the trap model in the aging regime: Exact results in the strong-disorder limit[END_REF] (see for example Theorem 11.3 p118 in [START_REF]Revész Random walk in random and non-random environments[END_REF]), and we obtain for such n,

Using the large deviation upper bound similarly to the lines below [START_REF] Popov | Random walk attracted by percolation clusters[END_REF], we obtain from the last inequality that P -almost surely, for n large enough,

S n 1 2 n 1 2 + β 2ξ -νβ-η .
Since ν and η can be chosen arbitrary small, this last inequality ends the proof.

Proof of Theorem 2. We first assume that β > ξ. From Lemma 3 and Lemma 4, we know that under assumptions of parts 1 or 2 of Theorem To conclude in the case β > ξ, note that under assumptions of part 1, ln |Xn| ln n converges P -almost surely to 1 and under assumption of part 2 and 3, P -almost surely, lim sup n→+∞ ln |Xn| ln n = 1 2 by the law of iterated logarithm. To extend the results to the border case β = ξ, we use the property (8) of the coupling, which implies that the long-time limit of ln |Yt| ln t is non-increasing in β. This completes the proof of part 1 with β = ξ. For the other parts, we use the property [START_REF] Ben Arous | Scaling limit for trap models on Z d[END_REF], that we will prove independently below. Again, the results claimed for β = ξ in parts 2 and 3 follow from the monotonicity of the coupling.

Proof of Theorem 3. First observe that when λ = 0, Since the martingale Z ǫ has jumps of size ǫ -1/2 tending to 0 and since its bracket converges to a deterministic limit, it is well known (e.g. Theorem VIII-3.11 in [START_REF] Jacod | Limit theorems for stochastic processes[END_REF]) that the sequence (Z ǫ , ǫ > 0) converges to the centered Gaussian process with variance tΣ, yielding the desired invariance principle under P ω . We now prove [START_REF] Ben Arous | Scaling limit for trap models on Z d[END_REF]. Since λ = 0 we have lim sup n ln |X n |/ ln n = 1/2, P -a.s., and since Ee βC 0 < ∞ it holds a.s. lim t ln S -1 (t)/ ln t = 1. This implies the claim.

f (Y t ) - t 0 (2d) -1 e -βC Ys e∼0 f (Y s + e) -f (Y s ) ds
Concluding remarks: (i) Part 2 of Theorem 2 deals with the upper limit in the subdiffusive case λ = 0, β > ξ. We comment here on the lower limit. In dimension d 3, n -1/2 (ii) A natural question is: What does the environment seen from the walker look like in the subballistic case? In fact, the prominent feature is that the size of surrounding cluster is essentially the largest one which was visited so far. Consider for instance the case of positive λ. One can prove that, for β > ξ and ǫ > 0, 1 t s ∈ [0, t] : (ln t) -1 C Ys ∈ [β -1ǫ, β -1 + ǫ] -→ 1 P -a.s. as t ր ∞.

(iii) We end the paper with a short comment on the case when the environment is a general random field, not necessarily coming from site percolation. It is easy to check that Theorems 1 and 3, together with their proofs, remain valid for a stationary, ergodic random field (C x , x ∈ Z d ). On the contrary, our proof of Theorem 2 uses some independence property specific to the percolation model.

Theorem 2 .

 2 (Subballistic regime) Let β ξ.

Lemma 3 .

 3 Assume β > ξ. For d 2 or λ > 0, we have lim inf n ln Sn ln n β ξ , P -almost surely.

Lemma 4 .

 4 Assume β > ξ. For d 1 and λ 0, we have P -almost surely, lim sup n ln Sn ln n β ξ .

Lemma 5 .

 5 Assume β > ξ. For d = 1 and λ = 0, we have P -almost surely, lim sup n

E

  y,k ) α .

0 ee

 0 is a P ω -martingale for f continuous and bounded. Then, for all ω, the process Y is a square integrable martingale under the quenched law P ω . Its bracket is the unique process Y taking its values in the space of nonnegative symmetric d × d matrices such that Y t Y * t -Y t is a martingale and Y 0 = 0. We easily computeY t = t -βC Ys ds × d -1 I dBy Corollary 1, we see that the bracket Z ǫ is such that, for all t 0, -βC Ys ds × d -1 I d -→ tΣ as ǫ ց 0 P -a.s., and then in P ω -probability for a.e. ω. Let us fix such an ω, and use the law P ω .

  |X [ns] | converges to a transient Bessel process, and it is not difficult to see that In dimension d 2, X is recurrent, and then lim inf t |Y t | = 0 and lim inf

	lim sup t→∞	ln |Y t | ln t	= lim t→∞	ln |Y t | ln t	= ξ/(2β)

t→∞ ln |Y t | ln t = -∞

Acknowledgement: We thank Marina Vachkovskaia for stimulating discussions on the model.

* Partially supported by CNRS (UMR 7599 "Probabilités et

For any y in Z and n in N, we define u y,n = 1 n 1 2 +ν n 1 2 +ν k=0 E y,k . Fix µ > 0, according to the large deviation principle for i.i.d. sequences, there exists I µ > 0 such that, for any y in Z and any n in N,

Using the independence of the (E y,k ) y∈Z,k∈N , it is easy to check that Q(∃y ∈ [-n

is the general term of a convergent series and using Borel-Cantelli's lemma we obtain that Q-almost surely, for n large enough and for any -n

From the ergodicity of the environment, it is true that P-almost surely,

Using now ( 14),( 15) and ( 16), we obtain that for any α < ξ/β, there exists M < +∞ such that, P -almost surely for n large enough,

Since the last inequality is true for η and µ arbitrary small and α arbitrary close to ξ/β, the proof is complete. Proof of Lemma 6. Let η and ν be two positive real numbers. As a consequence of ( 12) and ( 13), P -almost surely, for n large enough, at least n 2 ( 1 ξ -ν) ln n points (y i ) i αn in O n such that the intervals (I i ) i αn centered in (y i ) i αn and of length 1 2 ( 1 ξν) are disjoint. If all sites of an interval are open, it will be said open, otherwise it will be said closed. Using the fact that the (I i ) i αn are disjoint, we obtain that,

As a consequence of Borel-Cantelli's lemma we obtain that P -almost surely, for n large enough, there exists at least one site visited more than n 1 2 -η times and that belongs to a cluster of size greater than 1 2 ( 1 ξν) ln n, we will note this site y n and therefore,

-νβ E e yn,i .


Jean-Yves Blaise

Iwona Dudek

Concentric time: enabling context + focus visual analysis of architectural changes

In order to acquire and share a better understanding architectural changes, researchers face the challenge of modelling and representing events (cause and consequences) occurring in space and time, and for which assessments of doubts are vital. This contribution introduces a visualisation designed to facilitate reasoning tasks, in which a focus view on evidence about what happens to artefact λ at time t is a complemented with a context view where successive spatial configurations of neighbouring artefacts, durations of changes, and punctual events are correlated and tagged with uncertainty markers.

Introduction

Analysing, representing, cross-examining what is really known about transformations occurring during the lifetime of historic artefacts is a research topic where issues specific to historic sciences often meet models and formalisms developed in computer science. Although largely outnumbered by self-evident applications of computer graphics, fruitful knowledge representation-oriented research efforts do exist on the modelling of such spatial dynamics.

But gaining insight (in the sense of [START_REF] Spence | Information vizualisation[END_REF]) into architectural changes is not only about circumscribing the relevant pieces of knowledge -i.e. abstracting relevant properties of the reality observed, and foreseeing possible processing. In historic sciences particularly, it is also basically about giving users means to shed a new light on their knowledge, to draw individual, self-achieved, maybe temporary conclusions -i.e. graphics as a discovery tool [START_REF] Bertin | Sémiologie graphique[END_REF]. And indeed an effective way to support reasoning is to capitalize on the user's capacity to use vision to think [START_REF] Card | Information vizualisation: Using vision to think[END_REF].

Now what do we here need to reason about? Architectural changes occur in space, and in time -analysing changes will mean handling both dimensions. Changes may have implications on the shape of an artefact or not, anyway the evidence is liable to be distributed in space/time slots. Of course, visual solutions that help distributing and analysing information in space exist, either in 2D or 3D [START_REF] Sanders | Models in Spatial Analysis London[END_REF]. In parallel, other visual solutions may help us to distribute events in time, using the basic timeline paradigm, or more sophisticated concepts [START_REF] Dürstener | Visualising Time[END_REF] [START_REF] Jensen | Semantic Timeline Tools for History and Criticism[END_REF]. But architectural changes do not occur in space, on one hand, and in time, on another hand. If we want graphics to help us perform reasoning tasks, we need them to combine a spatial representation and a chronology representation 1 . So-called "time sliders" (nested inside 2D or 3D interactive graphics) could be seen as an answer: they allow users to select a time slot and observe the corresponding spatial configuration. But what the user sees is one period at a time: in other words, a focus view. And beyond merely pointing out facts -artefact λ changes at time t -amplifying cognition [START_REF] Kienreich | Information and Knowledge Visualisation: an oblique view[END_REF] about architectural changes implies a teleological approach through which we connect facts about an artefact to a context -previous and following changes, things that happen in the neighbourhood, events that occur during that same period of time, etc. Besides a focus view, we need a context view where successive spatial configurations, durations of changes, punctual events, etc., are correlated. Finally, as usual in historic sciences, the evidence we base on is uncertain: efficient graphics should help visualising that uncertainty 2 .

In previous contributions, we focused on modelling issues, starting from two ends: the spatial bias (modelling the lifetime of artefacts as chains linking version to one another, with identification/classification issues) and the temporal bias (modelling key moments in the evolution of artefacts and processes of transformation in a cause + consequence approach). As a result, we developed solutions ranging from knowledge visualisation to information systems where either spatial interfaces or timeline-like graphics helped visualise architectural changes [8][9]. But at the end of the day, although our KR effort did integrate time and space, means to perform visual analysis privileged either time or space. In this research we present an experimental visualisation-called concentric time -designed according to four simple ideas:

• a combination of spatial features and chronology allowing comparisons enforced within the eyespan [START_REF] Tufte | Visual explanations[END_REF],

• an application of the context+focus principle, • a support for uncertainty assessment (fuzziness, impreciseness, lacks, etc),

• a simple (in the sense of [START_REF] Maeda | No simplicity without complexity[END_REF]) visualisation minimising the learning curve and the decoding effort, facilitating information discovery (including by a public of non-experts). It is applied and evaluated on the market square in Krakow, Poland, a 200*200m urban space where 24 artefacts have been built, modified, and for most of them destroyed, over a period of 750 years. Primarily designed as a visualisation, it also acts as an interface. Section 2 details where we start from. Section 3 presents the concentric time visualisation itself. In section 4, its evaluation is discussed. Finally, in sections 5 and 6, we point out some of its limits, and further possible investigations.

Combining spatial information and a chronology inside one unique representation is far from being a new problem. In a way, it is a feature of Roman "ribbon maps", where localities across the empire are connected through lines that mean travel durations, and not actual distances. Far before the "computer era", XIXth century scientists developed inventive and sometimes brilliant solutions like Minard's "figurative maps", to this day still regarded as exemplary [START_REF] Friendly | Re-visions of Minard[END_REF]. John Snow's analysis of the 1854 London cholera epidemic, as masterly analysed by E.R Tufte [START_REF] Tufte | Visual explanations[END_REF], shows how reasoning both on time and space, in a casual manner, can been vital. Yet we today expect from graphics features that paper-based solutions cannot offer (updatability, interactive browsing, interface capabilities, etc.). Let us still remember to counterbalance the natural verbosity and jumble-hungriness of computer-based solutions with the clarity of mind of the above mentioned precursors.

In this section we briefly present some of the prominent solutions to combine space and time explicitly, at visualisation level. Our intent is not exhaustiveness -unreachable here (see for instance [START_REF] Sanders | Models in Spatial Analysis London[END_REF]) -but to say basically where we start from.

Originating from geosciences for the former, from engineering sciences for the latter, GIS platforms and CAD tools can allow users to display at a given position p successive versions of an object using the "layer" concept (although this concept may take a different name depending on the actual software). And indeed, the layering of time-related info may be efficient in terms of information delivery. Each version of the object can correspond to a given time slot, and so the job is done -space and time are present. Well apparently the job is done, but apparently only: when giving a closer look time does not really exist here as such (with for instance vital parameters such as durations of changes not present). Versioning offers a focus view, not a context view: what you get to see are moments in an object's evolution, hardly its whole history. A number of interesting research works introduce, at modelling level, ideas to help overcoming weaknesses inherent to the abovementioned solutions, like [START_REF] Rodier | Social features, Spatial features, Time features : An urban archaeological data model[END_REF]. Nevertheless, GIS-inspired solutions have found a wide audience among archaeologists, with a recurrent trend to try and implement "4D" information systems about the history of a site [START_REF] Johnson | Mapping the fourth dimension: the TimeMap project[END_REF] [START_REF] Van Ruymbeke | Development and use of a 4D GIS to support the conservation of the Calakmul site[END_REF]. But GIS and CAD tools are basically concerned with distributing things in space. Consequently, x and y axes of the representation are requisitioned for spatial data -and time has to manage with what is left. Time geography may have started from the same observation -the z axis is here used to superimpose on a basic x,y cartography movements across a territory [START_REF] Yu | Revisiting Hägerstrand's Time-Geographic Framework for Individual Activities in the Age of Instant Access[END_REF]. But the primary goal of time geography is to record and represent movements occurring in space over short periods of time. In our case artefacts do not "move" (or rarely…): they just change, and over long periods of time. Time geography is a promising concept -its adequacy to represent architectural changes remains questionable.

Widespread in information sciences, timelines, time bars, time charts are used to represent dates, periods, etc. A wide attention is notably put on the issue of visualising time-related phenomena in the field of information visualisation (see for instance [START_REF] Spence | Information vizualisation[END_REF], [START_REF] Tufte | Visual explanations[END_REF]). Such graphics do help reading a chronology, with events, durations, and possibly uncertainty assessments. A number of generic solutions have been developed in the past years, with some convincing results like [START_REF] Sabol | Visualizing Temporal-Semantic Relations in Dynamic Information Landscapes[END_REF]. But because time is represented in space -typically with a line or a spiral -spatial features tend to be scattered here and there, causing visual discomfort rather than helping to understand spatial changes. The chronographs experience [START_REF] Blaise | Profiling artefact changes: a methodological proposal for the classification and visualisation of architectural transformations[END_REF] taught us that although a timeline-like visualisation does greatly enhance reasoning about architectural changes, it is not well suited to spatial reasoning.

With computer implementations, interactive time sliders have been successfully introduced in many research and/or communication works [START_REF] Blaise | From artefact representation to information visualisation: genesis of informative modelling[END_REF]. As mentioned before, time sliders allow users to choose time slots they are interested in, and investigate the corresponding spatial configuration. But paradoxically, time sliders are mainly efficient in helping us to understand spatial changes. In fact, they are fundamentally a mean to browse through versions of objects using time as selection criteria, and in no way a context+focus visualisation. In a recent experiment we added to this time slider mechanism a visualisation of densities of changes: this context view however only focuses on time aspects -spatial changes remaining readable only one at a time.

As alternative, visual metaphors can be proposed that link spatial changes and a chronology -T.Ohta presents some nice concepts in [START_REF] Ohta | Diagram design course, In Informational Diagram Collection[END_REF]. In previous works, we developed a 3D metaphor called Infosphere, and a 2D metaphor called "ladybug race" that were equally "enjoyable". But neither the former (no durations, 3D navigation disconcerting for some users) nor the latter (no neighbourhood relations) did solve the problem. And so at the end of the day we are left with a number of solutions that correspond only in part to what we expect from a visualisation of dynamics of change.

The concentric time visualisation can be seen as coming to an arrangement between some of these solutions. It capitalizes on basic, well-known, well understood representations (cartography, timelines). It differs from them by introducing a specific, con-text+focus, graphic combination of spatial features and chronology. It is designed to help researchers reconsider their knowledge and the doubts that travel with it, and puts space and time on equal terms.

The concentric time visualisation: principles, implementation

To start with, let us introduce some of the terms that we will use in this section:

• the evolution of an artefact is the time span separating its creation from its extinction (i.e. full physical removal, including of sub-structures, [START_REF] Blaise | Profiling artefact changes: a methodological proposal for the classification and visualisation of architectural transformations[END_REF]); • morphological changes imply transformations of the artefact's "shape",

• recurrent changes may concern upkeep, ownership, function, etc.,

• contours are simplified representations of an artefact's plan,

• uncertainty may concern the dating of events, the very existence of events, or the physical reality of an artefact. In all cases it is a consequence of the evidence we base on. It is represented in this experiment by a simplified lexical scale (known, uncertain, hypothetical). Figures below present the visualisation's construction as a series of steps introducing spatial or temporal information -these steps should not be interpreted as "successive actions inside the construction procedure". The left part of figures illustrates the idea, the right part what it looks like when applied to our real case (for enhanced readability of the paper output, original contrast and line thicknesses of the computer output are accentuated). draw Fig. 1. A 2D map is displayed inside a red circle, showing simplified outlines of the artefacts under scrutiny. These contours are there to get an idea of where the artefacts were located, but they are not exact contours 3 . With no reasonable justification, an old map acts as a background. d,j). When the dating of a period is uncertain, coloured rectangles as well as circles are filled in white, and only outlined in colour (i, j). When the dating of a period is hypothetical, the outlined is dashed (e, g). When the very existence of an artefact is uncertain, the size of coloured rectangles is modified (halfwidth, g). Visible on the right figure, an arc marks an artefact's destruction (c). It is withdrawn from the representation when zooming in, in order to lower the overall graphic weight. Timelines for artefacts that remain up to now are continued up to the central circle (a'). Capitalizing on the graphics' radial structure, dotted grey lines connect the timeline to a position on the central circle; white lines extend short timelines backwards in time (b). Depending on what really happened, the actual geometry of the new contour may be left unchanged from version to version (a -adding a new storey for instance does not change the artefact's contour). However, the visualisation acts also an interface, each contour is connected to queries that are specific to it -allowing users to retrieve info (basic data, 2D/3D content, bibliography, etc.) about what really happened (this point is re-discussed in the perspective section). The timelines radiate around the central circle, but the contour's orientation to the north is maintained whatever angle the timeline is drawn at. Fig. 5. Finally, in order to allow basic spatial analysis, each of the contours can be projected inside the central circle on user selection. Because contours are drawn with a level of transparency, this mechanism enables user-monitored comparisons of contours at two levels: -Left, Neighbourhood analysis -in order to uncover the likelihood of event propagation -fires for instance, a plague at those periods.

-Right, Transformation analysis -shows the density of use of the ground over time, or helps spotting the period of emergence of a given component, and check on neighbouring artefacts whether "there is room for the added component" at that moment in history. When contours are projected, the central circle's content is simplified (background image is hidden as well as simplified contours) so as to avoid visual overload.

The concentric time visualisation has been applied to the evolution of Krakow's Main Square (24 artefacts, 79 "contours", 391 "events" over 750 years). Although we already knew quite well this site, the visualisation did help us spot some interesting patterns -like an unexpected pattern of uncertainty on XIXth century destructions in the north-west corner. The following figures illustrate on two examples how the visualization may support reasoning tasks. Fig. 6. Left, a focus view helps analysing causal relations between events and changes on neighbouring objects note parallel recurrent changes in both artefacts (red dots) independently of morphological changes. Right, a focus view allows comparisons of patterns between objects: compare densities of changes for (a) , (b) to this of (c) (a prison). The implementation is a fairly primitive one, combining two main levels:

• Reusable components, composed of an RDBMS where information for each artefact are stored (dating, sources, description of each of the artefact's change) and of XML files (geometrical information for each version of an artefact. These components, developed in previous experiments [START_REF] Blaise | From artefact representation to information visualisation: genesis of informative modelling[END_REF][9], were here only updated with some recent archaeological evidence. • The visualisation itself -an SVG file [START_REF] Rathert | Knowledge visualisation using dynamic SVG Charts[END_REF] produced on the fly 4 for standard web browsers (Perl classes -architectural ontology -and modules developed in the abovementioned experiments), with user interactions monitored inside Javascript modules (zoom/Pan, show/hide, contour projection, queries, etc.). Finally, the visualisation also acts as an interface: specific menus are available either over contours or over events that allow the querying of various data sets (RDBMS, XML files, as well as static 3D models corresponding to each contour).

Evaluation, limitations, perspectives

The evaluation was carried out with a group of eleven students in mechanical engineering (no background in architecture or historic architecture, no knowledge about the case study). Willingly, the session lasted less than 30 minutes, including the time needed for us to present the visualisation's objective and graphic codes. The visualisation was projected on a white screen; with zoom levels fitted to the various questions.

Testers were asked 14 questions distributed into 5 groups (space reading, chronology, uncertainty, comparing/analysing individuals, correlating space and time).

In the first three groups of questions we tested the visualisation's readability: testers were asked to decode the visualisation and to write down "how long does transformation t lasts", "which in this group of artefacts does not fit", etc. Facing questions of the two last groups, testers had to do some reasoning with questions like "identify pattern p", "can modification on artefact λ be related to events in the neighbourhood", etc.

Results show an error rate of 3%, mostly on readability questions. With such a low error rate, our questionnaire may be considered as not demanding enough. But answers were easy to give provided the visualisation was easily and unambiguously understood: and this is precisely what we were eager to verify. The evaluation does indicate the visualisation itself is almost self-evident, and moreover that some domainspecific notions (levels of credibility or uncertainties in particular) were correctly grabbed by full beginners in the field. We make no other claim about the evaluation. What has to be said is that the visualisation is a result of the interpretation of more than 400 bibliographic sources. Accordingly an evaluation that would allow comparing "reasoning from sources" and "reasoning with the visualisation" is far from being easy to build. But it is clear that we do not consider this evaluation as a definitive one. A number of issues remain open, for which future works will be needed:

-Case-study bias. The spatial layout of artefacts considered is our case study is well suited to the visualisation (artefacts kept close to one another for instance). A more complex spatial configuration, with more artefacts, with artefacts scattered in a wider area, or with more complex overlapping, may result in less convincing visual results. Tests on other case studies should be considered.

-Some additional features would definitely improve the visualisation's usability: independent zooming on spatial/temporal data, multi-resolution and partial zooming on timelines for short events, alternative levels of details for contours, etc. Our objective was not to come out with a ready-to-use system, but to test a concept: it is possible that this concept may be better implemented with a different technical set.

-In this experiment, only contours of artefacts are offered for users to carry out spatial analysis. Spatial analysis should be broadened to other 2D/3D characteristics.

-What the visualisation really does is nothing more than placing time and space side by side, and having them interact. Whether a deeper integration would bring better results remains to be established.

Conclusion

We introduce a context+focus visualisation called concentric time aimed at summarising the journey through time of groups of artefacts. The visualisation can be used for research purposes, but was evaluated with as ulterior motive testing its usability for a wider public, with possible museology applications. It was designed in order to combine inside a unique information space, spatial features and chronology, with as constraints uncertainty assessment, interface capabilities, and simplicity. The concentric time visualisation has been applied and evaluated on a case study, thanks to which we checked its support in carrying out reasoning tasks about architectural changes. Its limitations are numerous: what were investigated here really are the possible benefits of the concept-no claim is made on a generic, ready to use, system. On the other hand, the concept appears as fairly generic: its usability beyond clarifying architectural changes could be investigated (for instance in geoinformation). Providing models and visual tools to handle dynamics of change remains today a hot research topic: we view our contribution as demonstrating that besides facing sometimes complex knowledge modelling challenges, researchers in the field also face the challenging complexity of simple visual thinking.

Fig. 2 .

 2 Fig. 2.Starting from the red circle (that stands for the current year -2011), concentric circles represent a move towards the past. A circle is drawn for each decade (a) or century (b), except during the XXth century (this for readability purposes, considering very few changes occur then in our case study). The farer a circle is from the center, the older it is.

Fig. 3 .

 3 Fig. 3. Artefacts are connected to radial timelines (a, a') that distributes events in time and recount changes through two graphic variables: -Four colours : yellow (periods of stability, no known change, f), blue (morphological change, including construction, k), red (recurrent change, d), grey (destruction, c, meaning here that no visible trace is left above the ground).-Shapes: differentiate long-lasting events (duration > year, represented by rectangles -I, e, k) from punctual events (duration < a year, represented by a circle -d, j). When the dating of a period is uncertain, coloured rectangles as well as circles are filled in white, and only outlined in colour (i, j). When the dating of a period is hypothetical, the outlined is dashed (e, g). When the very existence of an artefact is uncertain, the size of coloured rectangles is modified (halfwidth, g). Visible on the right figure, an arc marks an artefact's destruction (c). It is withdrawn from the representation when zooming in, in order to lower the overall graphic weight. Timelines for artefacts that remain up to now are continued up to the central circle (a'). Capitalizing on the graphics' radial structure, dotted grey lines connect the timeline to a position on the central circle; white lines extend short timelines backwards in time (b).

Fig. 4 .

 4 Fig. 4. All along the timeline, a new contour is drawn for each "blue" (morphological) transformation. The contour is drawn either on the timeline or above/below the timeline (a point we re-discuss in the implementation section).Depending on what really happened, the actual geometry of the new contour may be left unchanged from version to version (a -adding a new storey for instance does not change the artefact's contour). However, the visualisation acts also an interface, each contour is connected to queries that are specific to it -allowing users to retrieve info (basic data, 2D/3D content, bibliography, etc.) about what really happened (this point is re-discussed in the perspective section). The timelines radiate around the central circle, but the contour's orientation to the north is maintained whatever angle the timeline is drawn at.

Fig. 7 .

 7 Fig. 7. Left, context view underlines features of the collection: for instance the proportion of remaining artefacts (a), 3 out of 24, and the impact of mid XIXth century destructions (b).

Classic gothic cathedrals in Bourges or Chartres were built in some decades, starting from the end of the XIIth century. The construction of the gothic cathedral in Tours started at the same period, but lasted until the XVIth century. And it is precisely that difference in time that explains some of the differences in the shapes of these edifices, i.e. differences in space.

The information can be imprecise to the extent that although an edifice λ did exist between time t1 and t2 in a given "area", we cannot precisely position it inside this area. However, un-localised e may have had an influence on others around. For instance, in our case study, an edifice called "Smatruz" (an open market hosting traders temporarily), localised "probably somewhere in the south west corner", was destroyed in mid XVth century. Its function being still needed, it is apparently drifted to nearby "Sukiennice", causing changes on this latter edifice, i.e. no reasoning on "Sukiennice" without representing un-localised "Smatruz".

There are several good reasons to this choice: -A number of artefacts were built over another, older one -in other words at the same x,y position. Drawing contours one over the other would result in a good degree of visual pollution. Accordingly contours are represented in such a way as to minimise overlapping. -If we were to represent a given version of an artefact on loading of the system, we would have to privilege one version of an artefact over another -and there is no reason to do so. Furthermore, if we were to privilege version 1 for artefact λ and version 3 for artefact δ, we could end up with showing side by side edifices that never were there at the same moment.So that means we should rather select a given moment in history -well in that case we would just not see the whole set of contours since they never were all present at the same time. In short, showing approximate contours enhances readability, and avoids spatial or temporal inconsistency.

A puzzling problem was how to position contours along the timeline (across, above, below? -see Fig4.). No systematic solution appeared as satisfactory. For each contour we specify manually, inside a javascript "configuration" file, an x,y move from its "natural" position (which is the center of a blue rectangle / circle). The visualisation is dynamically written with "natural" positions, and an onload jscript event moves the contours according to the configuration file. Contours can also be freely re-positioned by users (mouse dragging).