Nicolas Chrétien
email: nicolas.chretien@femto-st.fr

Bruno Franc
email: bruno.francois@femto-st.fr

Gilles Martin
email: gilles.martin@femto-st.fr

Ballandras Sylvain

Time & frequency Femto-St

Jean-Michel Friedt
email: jmfriedt@femto-st.fr

Pratyush Varshney
email: varshney.pratyush@gmail.com

On the Need for Low Phase Noise Oscillators for Wireless Passive Sensor Probing

Keywords: surface acoustic wave, sensor, passive, wireless, phase noise, RADAR

Two oscillators are needed for passive wireless sensor readers: a radiofrequency local oscillator generating the carrier within the bandpass of the sensor, and a clock triggering analog to digital conversion of the signals returned by the sensor. We assess the influence on measurement resolution of these two oscillators, hinting at some design rules of the sensor based on the characteristics of the oscillators as driven by cost, size, or power consumption. We demonstrate that local oscillator phase noise is a significant parameter in assessing the resolution of passive acoustic sensor probed through a wireless link, although with different characteristic conditions whether the sensor is in a delay line configuration (short term -far from carrier -response) or a resonator (long term -close to carrier -response).

I. INTRODUCTION

Wireless passive sensors are either piezoelectric or dielectric transducers coupling with an incoming electromagnetic field following conditions dependent on the physical property under investigation. For instance, surface acoustic wave (SAW) delay lines convert an incoming electromagnetic pulse to a mechanical wave propagating on a piezoelectric substrate. Mirrors patterned on this substrate reflect a fraction of this wave back to the interdigitated transducer (IDT) connected to the antenna: the direct piezoelectric effect converts these acoustic pulses to electromagnetic signals detected by the receiver. Hence, a passive acoustic delay line reader operates following principles similar to RADAR, with a delayed echo not associated with reflections of the emitted signal over dielectric or conductive interfaces, but with delays associated with a measurement. Thus, all the well known RADAR techniques have been applied to passive wireless sensing, whether dielectric [START_REF] Allen | The use of groundpenetrating radar with a cooperative target[END_REF], [START_REF] Sridhar | A hydrogel-based passive wireless sensor using a flex-circuit inductive transducer[END_REF], [START_REF] Jatlaoui | Original identification technique of passive em sensors using loaded transmission delay lines[END_REF] or based on piezoelectric substrates [START_REF] Pohl | Radio signals for SAW ID tags and sensors in strong electromagnetic interference[END_REF], [START_REF] Reindl | SAW devices as wireless passive sensors[END_REF], [START_REF] Bulst | State of the art in wireless sensing with surface acoustic waves[END_REF], [START_REF] Stelzer | A low-cost interrogation unit and signal processing for a SAW identification-tag for a pressure sensor[END_REF], [START_REF] Bruckner | A hightemperature stable SAW identification tag for a pressure sensor and a low-cost interrogation unit[END_REF], [START_REF] Schuster | Performance evaluation of algorithms for SAW-based temperature measurement[END_REF]: wideband pulsed RADAR [START_REF] Wen | Detecting and evaluating the signals of wirelessly interrogational passive SAW resonator sensors[END_REF], [START_REF] Hamsch | An interrogation unit for passive wireless SAW sensors based on Fourier transform[END_REF], FMCW RADAR [START_REF] Stelzer | Wireless SAW sensors for surface and subsurface sensing applications[END_REF], [START_REF] Beckley | Noncontact torque senseor based on SAW resonators[END_REF], FSCW RADAR [START_REF] Scholl | Wireless passive SAW sensor systems for industrial and domestic applications[END_REF], [START_REF] Bruckner | A hightemperature stable SAW identification tag for a pressure sensor and a low-cost interrogation unit[END_REF], [START_REF] Schuster | Performance evaluation of algorithms for SAW-based temperature measurement[END_REF], and chirped RADAR [START_REF] Springer | A robust ultra-broadband wireless communication system using SAW chirped delay lines[END_REF], [START_REF] Puccio | SAW sensors using orthogonal frequency coding[END_REF].

Since it is well known that the local oscillator characteristics drives the detection capability of RADARs as will be discussed in the first introductory section of this presentation (Fig. 1), one might consider how local oscillator phase noise In case of a static target, a CW RADAR receiver noise detection limit is associated with the local oscillator frequency fluctuation between the emitted pulse at time t and the received pulse delayed by τ , the two-way transit duration.

affects passive wireless sensor resolution [START_REF] Kalinin | Influence of receiver noise properties on resolution of passive wireless resonant SAW sensors[END_REF]. We extend in the third section the discussion from the classical passive target to the short-term -wideband -delay line acoustic sensor configuration in which the phase noise characteristics far from the carrier defines the measurement resolution. Since wideband acoustic delay lines are only compatible with the allocated 2.45 GHz band, we consider in the fourth section the same approach applied to narrowband resonators, compliant with the narrowband 434 MHz radiofrequency band, and the phase noise characteristics now shifted to the region close to the carrier. Such considerations will bring us in the fifth section to consider a second oscillator usually found on such circuits: the clock defining the analog to digital conversion rate. Thus, the reader is led throughout this paper to consider the various regions of the phase noise spectrum as a limiting factor for acoustic sensing resolution depending on the transducer characteristic time constants.

II. PHASE NOISE INFLUENCE ON CW RADAR

This introductory section reminds the reader of basic concepts related to phase noise of oscillators and their effect on RADAR detection capability. We will focus on the continuous wave (CW) RADAR where the explanation is straight forward.

CW RADARs are used whenever a velocity information is considered without ranging capability: a radiofrequency (RF) wave is generated by an oscillator. This signal is on the one hand fed to an antenna (after being amplified by a power amplifier, PA) and on the other hand a fraction of the output of the oscillator is sent to one input of a mixer. The second input of this mixer is fed with the signal detected by either a second antenna in a bistatic configuration, or at the output of a circulator in a monostatic configuration, after amplification by a low noise amplifier (LNA). The output of the mixer m is the product of the frequency generated at time t by the oscillator but shifted by the Doppler frequency due to target motion δf , and the oscillator frequency delayed by a duration τ due to the electromagnetic wave propagation in air to and from the target:

m = cos (2π (f (t) + δf )) × cos (2π (f (t + τ ))) ∝ cos (2π (f (t) + δf ± f (t + τ ))) (1)
with only the difference term remaining after filtering the output by a low pass filter aimed at removing the signal at frequencies above f . Let us consider the case of slowly moving targets, where δf will be considered negligible; then m ≃ cos (2π (f (t) -f (t + τ ))). Ideally this term should vanish when the target is not moving and assuming the oscillator ideally stable, i.e., f (t) constant; the only beat frequency would be associated to δf . However, oscillators do exhibit phase noise, and thus f (t + τ ) and f (t) differ: the phase noise spectrum of an oscillator is defined as the Fourier transform of the autocorrelation function of the oscillator output frequency [START_REF] Rubiola | Phase Noise and Frequency Stability in Oscillators[END_REF]. The classical CW RADAR detection limit concludes that a moving target will only be detectable if its RADAR cross section is large enough so that the returned power (echo) is stronger than the power spectrum of the local oscillator: the phase variation is expressed in dBc/Hz, or a power with respect to the carrier power at an offset 1/τ from the carrier frequency. As a concluding remark, long range RADAR is interested in the behavior of the oscillator close to the carrier (since τ = 2d/c with d the distance to the target and c the velocity of an electromagnetic wave: d = 5 -50 km yields τ = 33 -333 µs in vacuum and thus the behavior of the oscillator at 3 to 33 kHz from the carrier is of interest). On the other hand, RADAR aimed at detecting moving walking people with targets in the sub-100 m range will only be affected by the phase noise above 1.5 MHz from the carrier.

III. APPLICATION TO SAW REFLECTIVE DELAY LINE

SENSING

One implementation of SAW delay line readers acts exactly as a CW reader: a carrier is chopped in pulses containing as many periods as there are electrodes in the sensor IDT (Fig. 2).

The reader on the one hand emits these pulses whose frequency is centered on the oscillator frequency output, and the returned signal from the sensor is centered on the same frequency, but shifted in time by a duration dependent on the physical property under investigation (which most significantly affects the acoustic wave velocity on the piezoelectric substrate). Thus, the mixer output exhibits a series of pulses whose rough delay is estimated through maximum returned power (threshold) or cross-correlation; but it is well known that only a phase measurement (with 2π uncertainty) provides the required high accuracy on the acoustic velocity and thus the measured physical quantity [START_REF] Kuypers | 2.45 GHz passive wireless temperature monitoring system featuring parallel sensor interrogation and resolution evaluation[END_REF], [START_REF] Reindl | Wireless passive SAW identification marks and sensors[END_REF]. Let us now add to the time delay from the acoustic wave propagation another contribution to the detected phase: the local oscillator intrinsic noise as characterized by its phase noise. The phase noise of a signal V (t) = (V 0 + ε(t)) sin (2πf 0 t + ∆ϕ(t)) is defined as [START_REF]Phase noise characterization of microwave oscillatorsphase detector method[END_REF] the phase fluctuations in a 1 Hz-wide bandwidth

S ∆ϕ = ∆ϕ 2 RM S
measurement bandwidth rad 2 /Hz and the classical representation of the noise spectrum is given by L(f ) = 1 2 S ∆ϕ (f ) = 10 × log 10

P SSB P S
dBc/Hz. Based on these informations, we will compute the phase noise fluctuations of the local oscillators during time intervals τ which are now given by the travel duration of the electromagnetic wave in the medium surrounding the sensor (negligible since readout ranges are in the tens of meters at most, or tens of nanoseconds) and the acoustic delay which is typically in the 1 to 5 µs range: the offset to the carrier of interest to acoustic sensing is in the 200 kHz to 1 MHz range.

This frequency range usually lies above the Leeson frequency f L = f LO /(2Q LO ) with f LO and Q LO the local oscillator resonator frequency and quality factor. This characteristic frequency defines a frequency offset from the carrier at which the resonator no longer acts as an energy tank and becomes transparent to the feedback amplifier noise. Above this frequency, the phase noise of the oscillator is constant and solely defined by the power injected in the resonator, the noise factor of the feedback amplifier and the operating temperature. We shall come back to such considerations in the design section.

Two practical applications will focus on a poor oscillator assumed to exhibit -130 dBc/Hz, and an excellent oscillator assumed to exhibit -170 dBc/Hz as the frequency offset of interest. Another numerical application using -90 dBc/Hz is justified by the fact that the returned signal noise floor is the maximum of either the initial oscillator noise floor raised by the power amplifier (PA) and low noise amplifier (LNA), or the LNA noise floor set by its thermal noise F LN A k B T /P R with F LN A ≃ 1.5 dB the noise factor of the reception amplifier, 10 log 10 (k B T ) = -174 dBm the product of the Boltzmann constant with the temperature T = 290 K, and P R the received power. From this consideration, the measurement resolution will first be constant as long as the LNA noise floor is lower than the LO noise floor, and drops once the returned power becomes so low that the LNA noise floor rises above the LO noise floor. The received power is related to the emitted power P E -limited to P E = +10 dBm by radiofrequency emission regulations in 434 and 2450 MHz ISM bands -through the free space propagation losses and the sensor insertion losses. Free space propagation losses F SP L = 4πdf c are associated with energy distribution on a sphere generated by the emitter, and in the case of a RADAR the link budget requires the use of F SP L 4 since the target itself acts as a point-like source generating a spherical wave. The SAW sensor insertion loss IL is a significant source of energy loss when probing SAW delay lines since a typical IL value is -35 dB. Thus, P R = P E × F SP L 4 × IL and switching to a logarithmic description, the noise floor on the return branch reaching the mixer is either the floor of the oscillator raised by the noise floor of PA and LNA, or the noise floor of the LNA amplifier F LN A,dB + 10 log 10 (k B T ) -10 log 10 (P E × F SP L 4 ) -IL). The lower the oscillator phase noise floor, the smaller the range at which the LNA noise floor becomes dominant, as shown in the numerical application of Table I.

In such cases, the phase variations due to the local oscillator are ∆ϕ RM S = √ 2 × 10 -(130..170)/10 rad/ √ Hz. Since we focus on measuring the phase within a 30 ns long pulse, the measurement bandwidth is 60 MHz and ∆ϕ RM S = √ 2 × 10 -(130..170)/10 × 60 × 10 6 rad whose numerical application yields to phase fluctuations from 0.2 o to 0.002 o (for -130 and -170 dBc/Hz cases respectively).

We must now relate these phase fluctuations with the phase variations due to a physical quantity variation: we focus on a temperature sensor. An acoustic sensor exhibits a phase rotation for every period, i.e., for a propagation length of one wavelength λ. The elastic wave propagates on the piezoelectric substrate at velocity v and the time-difference due to the two-way trip d from IDT to the mirror yields a phase shift of

∆ϕ = 2π × d/λ = 2π × d × f /v
The variation with temperature T of this phase difference is associated with the velocity variation, so that

∂∆ϕ ∆ϕ T = ∂v v T ⇔ ∂∆ϕ(T ) = 2π d × f v × ∂v v T
All quantities in this equation are known: for a LiNbO 3 substrate, we consider that v ≃ 3000 m/s, ∂v/v ≃ 60 ppm/K. Selecting d = 10 mm and f = 100 MHz (as used in [START_REF] Friedt | Surface acoustic wave devices as passive buried sensors[END_REF]), we conclude that 2π × 60 × 10 -6 × 10 -2 × 10 8 /3000 = 0.13 rad/K= 7.2 o /K.

By extending this analysis to various experimental parameters, we compare the local oscillator phase noise fluctuation implication on the measurement resolution in Table II.

We conclude that the local oscillator stability becomes a significant hindrance to high resolution temperature measurements, and reaching the mK resolution as was done with the Hewlett Packard HP2830A resonator-based probes is challenging.

Beyond the compliance with radiofrequency emission regulations, the use of ultra-wideband (UWB) interrogation strategies, e.g., Ground Penetrating RADAR based approaches [START_REF] Friedt | Surface acoustic wave devices as passive buried sensors[END_REF], yields the question of optimum operating frequencies. Indeed, we have seen that the time delay is a function of the acoustic velocity and propagation path length (defined respectively by selecting appropriate single-crystal piezoelectric substrate orientations, and design considerations in positioning the mirrors on the sensor surface), but also of the operating frequency:

∆ϕ = 2πd/λ = 2πdf /v = 2πf τ
where τ is the propagation duration of the pulse, i.e., ∂∆ϕ = 2πf ∂τ ⇔ ∂τ = 1/(2πf )∂∆ϕ, providing the relationship between phase noise and delay noise through the inverse of the frequency.

The remaning design issue lies in the selection of the echo pair used for computing acoustic propagation time delay and thus identifying the physical quantity under investigation. The first and last echos of a tag (start and stop bits) are usually considered for such purposes. However, the further away mirrors are, the longer the delay and thus the larger the local oscillator fluctuations, associated with phase noise rise. One should thus take care that the inverse of the propagation delay does not reach the Leeson frequency f L where the noise floor meets the rising phase noise slope: f L = f 0 /(2Q). Considering a (very favorable) Q = 20000 resonator used for generating a 2.45 GHz oscillator, f L = 60 kHz and the associated propagation delay is 16 µs, far above any practical limitation (such a delay would be associated with a 24 mm-long propagation path). However, for a more reasonable Q = 2000 [START_REF] Knauer | Design, fabrication, and application of GHz SAW devices[END_REF], the Leeson frequency reaches 600 kHz or a propagation delay of 1.6 µs. In this case, using echos returned by mirror at extreme positions of the delay line should be avoided (i.e., exhibiting propagation delays larger than 1.6 µs) and adjacent echos should yield results with higher resolutions.

IV. SENSOR DESIGN WITH RESPECT TO LOCAL OSCILLATOR CHARACTERISTICS

Consider two applications: a 434 MHz delay line (which would not comply with RF regulations) and a 2450 MHz delay line. Resonators in the former frequency range exhibit typical quality factors of 10000, and f L = 22 kHz, well below the 200-1000 kHz range we have been considering: a classical delay line design will be probed at best by the reader. However at 2450 MHz, since the product Q × f is constant for a given technology, the local oscillator Q drops to 1800, and the Leeson frequency rises to 680 kHz. The designer of a 2450 MHz delay line would be wise to avoid the phase noise rise below f L and limit the maximum time delay on the acoustic path to 1.5 µs. Since electromagnetic clutter fades within the first 700 ns (assuming 100 m range) and the typical pulse length is 40 ns spaced by at least 100 ns to account for manufacturing variability, limiting the delay to 1.5 µs still leaves enough space for 5 reflections, more than enough for multi-parameter-sensing (one reference pulse and 4 pulses for probing 4 different physical quantities, e.g., temperature [START_REF] Fachberger | Wireless SAW based high-temperature measurement systems[END_REF], pressure [START_REF] Buff | Universal pressure and temperature SAW sensor for wireless applications[END_REF] and two chemical compounds [START_REF] Dong | A multi-resolution passive SAW chemical sensor[END_REF].

V. APPLICATION TO SAW RESONATOR SENSING SAW resonator probing aims at identifying a characteristic frequency: in one embodiment of this approach, a frequency sweep network analyzer sequentially probes multiple frequencies in order to identify the frequency at which the sensor returns a maximum power. SAW resonators store energy during the electromagnetic signal emission phase, and release this energy (as an electromagnetic wave at the sensor resonance frequency f 0 ) during the listening stage: the time constant of each step is Q/(πf 0 ) with Q the sensor quality factor. The fastest approach to the best of our knowledge [START_REF] Friedt | Remote vibration measurement: a wireless passive surface acoustic wave resonator fast probing strategy[END_REF] for probing a resonance frequency of a resonator requires two signals at different frequencies, one above and one below f 0 (Fig. 5).

2Q/( f) π |S11| f Figure 5
. For a dual-mode resonator, required for a differential measurement, a minimum of 4 measurements each lasting 2Q/(πf ) seconds, with f the resonance frequency of one mode and Q its quality factor, is needed (red). A more classical approach of a frequency sweep network analyzer requires up to 128 measurements in the 434-MHz European ISM band (blue).

Hence, the minimum measurement duration is 2Q/(πf 0 ) for each probed frequency. For f 0 = 434 MHz and Q = 10000 in a dual resonator configuration, eight time constants (two resonators, and for earch two-measurement points, and for each one time constant for loading and unloading the resonator) yield 59 µs measurement duration, so that the oscillator stability at 17 kHz from the carrier is of interest. The phase noise S ∆ϕ and frequency noise (or stability) S ∆f at f from the carrier are related through

S ∆f = f 2 × S ∆ϕ = ∆f 2 RM S
BW So, with a measurement bandwidth BW of 2f , frequency fluctuations are given by

∆f RM S = BW × f 2 × 2L(f )
The result of this calculation is summarized in table III, assuming that the Q × f product is constant, as is usually considered for a given technology, with values representative of SAW resonators patterned on a quartz substrate. The temperature resolution -last column of Table III -is computed assuming a 60 ppm/K sensitivity. This last result scales as the temperature sensitivity of the substrate: a sensor allowing for a 170 K measurement range within the 1.7 MHz wide 434 MHz European ISM band -as sold by SENSeOR (Mougins, France) -only exhibits a 5.7 ppm/K sensitivity and hence the values in the last column are multiplied by 10.

Here again, for resonators acting as temperature sensors with 2.5 kHz/K temperature sensitivity (in order to fit a 170 K measurement range within the 1.7 MHz wide 434 MHz Industrial, Scientific and Medical (ISM) band, accounting for manufacturing variations), the 25 Hz frequency resolution (10 mK resolution) is only met if a reasonably stable local oscillator is used as reference, with a phase noise below -105 dBc/Hz. This result is consistent with the phase noise spectra provided in [START_REF] Friedt | A wireless interrogation system exploiting narrowband acoustic resonator for remote physical quantity measurement[END_REF], with a phase noise around -105 dBc/Hz in the 500-5000 Hz carrier offset range at the 434 MHz DDS output. The former range boundary is met when probing 128-samples in a frequency-sweep network analyzer approach: 128 points each requiring 2Q/(πf 0 ) requires a duration of 1.8 ms or an update rate of 530 Hz.

VI. ANALOG TO DIGITAL CONVERSION JITTER

We now change oscillator type to consider the analog to digital conversion (ADC) stage. The phase measurement requires two simultaneous measurements of the I and Q components of the returned signal after mixing with the local RF oscillator. The typical pulse duration is 30 ns so that the ADC bandwidth must be at least in the f s = 60 MHz range, or practically (3 points/period at least) 100 MHz.

Measuring a phase with 0.13 rad resolution over the full 2π range requires bits = 6 bit resolution. Since the jitter on the clock controlling the ADC yields a resolution loss (linear scale) of SN R = (2πf s σ t ), the jitter σ t must not exceed σ t ≤ 2 -bits /(2π × f s )

which is here equal to 42 ps [START_REF] Redmayne | Understanding the effect of clock jitter on high speed ADCs[END_REF]. However, increasing 10fold this resolution yields a 9 bit ADC resolution and a maximum jitter of 5 ps. On the other hand, let us estimate the jitter induced by an oscillator exhibiting a -130 dBc/Hz phase noise level in the 200 kHz-200 MHz range, representative of the influence of the clock controlling the ADC sampling at 100 MS/s for a maximum duration of 5 µs. The RMS jitter (in seconds) is given [START_REF] Kester | Converting oscillator phase noise to time jitter[END_REF] by

σ t = √ 2 × 10 -130/10 × 10 8 (2π × 10 8 )
which is equal to 7 ps, dropping the lower integration limit (200 kHz) by assuming that the constant phase noise level extends to the carrier. Thus, although even a very poor reference oscillator controlling the ADC meets the requirements of 9-bit resolution needed for high resolution temperature measurements, care should nevertheless be taken to reach sub-10 ps jitter. As an example, the Digital PLL generating the clock output of an iMX27 CPU as used on the APF27-board from Armadeus Systems (Mulhouse, France) for prototyping our experiments is specified at a maximum of 200 ps, hardly usable for the application described here [START_REF]MX27 and i.MX27L data sheet[END_REF].

VII. CONCLUSION While the debate on the advantages between delay line and resonator approaches is still ongoing, local oscillator characteristics brings some hint on which strategy might bring the most accurate result. From a local oscillator perspective, moving the frequency offset as far as possible from the carrier, i.e., allowing for as short a duration between various measurements of the sensor characteristics as possible, clearly hints at an advantage towards delay lines. However, this partial picture does not include the receiver noise level, especially the high bandwidth on the ADC sampling required to recover and digitize the fast delay line response: only an extremely stable (low jitter) clocking circuit for the receiver ADC will provide measurements with resolutions comparable to those of resonators. Furthermore, as opposed to FMCW or frequency sweep approaches which require tunable frequency sources (VCO, frac-PLL, DDS), a pulsed (UWB-like) delay line approach only requires a fixed frequency source generating a stable signal within the bandpass of the sensor, hence allowing for improved stability. Such results are most significantly the target of high quality factor piezoelectric resonator based oscillators aimed at reaching the targetted radiofrequency band.

Design rules concerning the oscillator characteristics are provided for delay lines: the maximum two-way trip duration should be lower than the inverse of the Leeson frequency, while only low noise floor enables high resolution measurements as explicitly stated with relationships between local oscillator phase noise densities and measured returned signal phase resolution. For resonator probed through a frequency sweep network analyzer approach, the tunable local oscillator source is clearly a limiting factor in the measured resonance frequency resolution.

Figure 1 .

 1 Figure1. In case of a static target, a CW RADAR receiver noise detection limit is associated with the local oscillator frequency fluctuation between the emitted pulse at time t and the received pulse delayed by τ , the two-way transit duration.

Figure 2 .

 2 Figure 2. Typical response from a SAW delay line, here from a temperature sensor with 8-bit coding sold by CTR Carinthian Tech Research (Villach, Austria), here excited by a 40-ns long pulse centered on 2.40 GHz. The pulse at 0 s is the excitation pulse, and the returned echos are located between 1 and 2.2 µs.

Figure 3 .

 3 Figure 3. Demodulation circuit for probing wireless passive SAW sensors. The mixer might be replaced by an I/Q demodulator in practical systems.

Figure 4 .

 4 Figure 4. Phase noise of a 2.45 GHz source generated by an Analog Devices ADF4360-0 Phase Locked Loop (poorly controlled), and Rohde & Schwartz SMA 100A tabletop frequency synthesizer set to 2450 and 434 MHz.


Nazim Fatès

Nikolaos Vlassopoulos

A Robust Scheme for Aggregating Quasiblind Robots in an Active Environment

Keywords: Decentralized Gathering, Agent Aggregation, Swarm Robotics, Reaction-diffusionchemotaxis

The question of how to aggregate autonomous agents with limited abilities in the absence of centralized control or coordination is known as the Decentralized Gathering Problem. We present a bio-inspired aggregation scheme that solves this problem and study a first application of this scheme to a small team of robots. The robots (Alice and Khepera III) obey simple rules and have only a rudimentary perception of their environment. The collective behavior is based on stigmergic principles and uses an active environment to relay the communications between robots. This results in an aggregation process that shows good properties of robustness and that can in principle be extended to swarms of robots.

INTRODUCTION

Let us consider a system consisting of a great number of autonomous robots randomly spread on a plane area, the environment. The robots need to gather to realize a given task, for instance to exchange information about the area or to collaborate to remove an obstacle. In the case where the agents are not aware of their position, orientation nor of the presence of other agents in their neighborhood and if furthermore there are no means of centralized coordination or communication, the question of their aggregation is known as the Decentralized Gathering Problem.

This problem has been studied before, using different assumptions about the agents and the environment. In an early work, [START_REF] Sugihara | Distributed Algorithms for Formation of Geometric Patterns with many Robots[END_REF] deal with the distributed coordination of a swarm of robots in order to form geometrical shapes. The work presented by [START_REF] Garnier | The Embodiment of Cockroach Aggregation Behavior in a Group of MicroRobots[END_REF][START_REF] Garnier | Self-organized aggregation triggers collective decision making in a group of cockroach like robots[END_REF] studies the emergent aggregation of robots that mimic the behavior of Blattela germanica, while [START_REF] Ando | Distributed memoryless point convergence algorithms for mobile robots with limited visibility[END_REF] and [START_REF] Flocchini | Gathering of Asynchronous Robots with Limited Visibility[END_REF] study the aggregation problem under the assumption of limited visibility. A theoretical study by [START_REF] Prencipe | Impossibility of Gathering by a set of Autonomous Mobile Robots[END_REF] explored the complexity of the problem of decentralized gathering with simplifying hypotheses.

One approach to solving the decentralized gathering problem consists of imitating the social amoeba species Dictyostelium discoideum [START_REF] Fatès | Solving the decentralised gathering problem with a reaction-diffusionchemotaxis scheme -Social amoebae as a source of inspiration[END_REF]. The global dynamics of this aggregation scheme were first studied in the scope of cellular automata (CA) modeling. However, this scheme is more generic as it can be considered as a high-level behavioral description. The key idea is the existence of an active environment that conveys simple messages among the agents, the virtual amoebae. The agents interact with the environment by either initiating the transmission of a message or by detecting the existence of messages in their local neighbourhood. These two ingredients form the basis for a stigmergic behavior.

The first advantage of this scheme is that both the agent description and the active environment are simple and straightforward to implement, at least in a simulation environment. The main idea behind the amoebae scheme is that the use of an active environment delegates to this environment parts of the complexity that would otherwise be implemented in the agents.

In previous simulations, the aggregation scheme was shown to be robust to different perturbations of the agents and the environment, e.g., to the presence of obstacles and noise in the movement of the agents. The first aim of this paper is to investigate how to transpose the aggregation scheme in the case of robotic agents. The second objective is to examine whether the robustness of the aggregation process is conserved with robots. We examine two different series of experiments that use physically different implementations of the same active environment. The first approach we will present is based on the ALICE micro-robots1 . These robots feature a compact design but have a restricted set of abilities: the robots we had at our disposal were equipped with only two sensors that could measure the light intensity of a source located above the robots. The second approach uses a modified version of the Khepera III robots, in conjunction with an "interactive table" (ROMEA) that has been specifically developed for conducting robotic experiments (Beaufort et al., 2011). Although these robots have superior abilities, we made a minimalistic use of their features: our goal is indeed to present a simple method for aggregating the robots, in particular we want the experiment to be set rapidly without having to tune and calibrate too many parameters.

This paper is organized as follows: Section 2 presents the aggregation scheme while Section 3 describes the implementation on the two different types of robots. Following, Section 4 presents the advantages and drawbacks that we observed during our experiments, as well as the differences we observed between the two types of robots and environment implementations. Finally, we discuss the results in Section 5.

THE VIRTUAL AMOEBAE AGGREGATION SCHEME

From Biology to Computer Science

The virtual amoebae aggregation scheme draws its inspiration from the behavior of the social amoebae species Dictyostelium discoideum. During its life-cycle, this species undergoes a transformation from a mono-cellular to a multi-cellular organism. This transformation occurs in different phases. Initially, and when the environment is rich in nutrients, the individuals are scattered. When the concentration of nutrients decreases, the amoebae start to aggregate so as to form a single body, which moves until it finds a suitable environment. Once such an environment is found, the organism transforms again, forming a fruiting body. The fruiting body releases spores in the environment, and the cycle starts again.

The aggregation phase of D. discoideum is based on chemotaxis, that is, the ability to move preferentially towards a source that emits chemical products. In nature, the production of the chemical agent (cyclic AMP) within each amoeba is strated when the nutrients in the environment are lacking. The cAMP is secreted internally by the amoebae, and then emitted, which causes neighboring individuals to get closer to the emitting cell by following the gradient of the concentration of cAMP. This behavior amplifies until it results in the formation of small clusters that eventually aggregate to form a single body.

Our amoebae aggregation scheme is based on a similar concept. At each discrete time step, an agent may emit a "chemical message" with a given probability. The message is then transmitted through the environment detected by other agents, which, in reaction, move towards the source of the message. Compared to the natural phenomenon, the amoebae aggregation scheme uses simplified hypotheses. In order to reduce the complexity of the agents, the agents do not follow the gradient of the transmitted messages, but rather simply detect the wavefronts (the reason for such a choice will be explained later). This allows for the agents to work in an "all-ornothing" fashion, which, in turn, simplifies the requirements for agent sensors, so that they need minimal calibration. Further, to simplify the design of the environment the "messages" do not use "classical" diffusion as their transport mechanism. Instead, the aggregation scheme is based on the existence of an active environment that allows the waves to travel at a constant speed and without any attenuation or modification. This allows us to simplify the agent-environment interaction, since an agent can initiate a "message" wave in the environment by making a local modification in the environment.

The Greenberg-Hastings Cellular Automaton

The dynamics of the environment we employ emulates the Greenberg-Hastings cellular automaton (GHCA). This model emulates reaction-diffusion processes and has the following advantages: (a) it is easy to implement and simulate in software, (b) the reaction-diffusion wavefronts travel over arbitrary large distances in the environment without any attenuation, (c) wavefronts skirt round obstacles. Formally, the discrete space is modeled with Λ= {1 L}× {1  L} as a square two-dimensional array of cells where L determines the environment size. Each cell c=  c x ,c y  of Λ is associated with a state in the set {0, ... , M}. We call 0 the neutral state, the states M -1,  ,1 the refractory states, and M the excited state. The state of a cell at time t is denoted by σ c t . The state of each cell is updated in discrete times according to a function of its current state and the states of the cells in its neighborhood. If we denote by E c t the fact that the neighborhood of a cell c at time t contains an excited state, then the state evolution of a cell, c is described by the following set of equations (Fig. 1):

σ c t+1 = { M, if σ c t =0∧ E c t σ c t -1, if σ c t ∈ M -1,,1 0, otherwise

}

In this work we consider three different types of neighborhoods: the 4-connected neighborhood, the 8-connected neighborhood, and a "circular neighborhood" that simulates the diffusion of isotropic waves. In all cases we assumed that the environment has free (or absorbing) boundary conditions; as a consequence, waves that reach the boundaries of the environment are absorbed.

In practice, this cellular automaton imitates the dynamics of reaction-diffusion processes: waves that start from excited cells expand outwards, so as to form wavefronts of excited cells followed by cells in refractory states (see Figure 2). When two wavefronts of excited cells collide, they annihilate.

Coupling of the Virtual Amoebae with the Environment

The basic model assumptions for the behavior of the virtual amoebae can be summarized as follows:

 Agents are autonomous and there is no direct communication among them.

 Each agent interacts only with its immediate environment (local interaction).  Agents are not aware of their position, orientation or relative positions of other agents in the environment (quasi-blindness).

With the above in mind, the behavior of each agent can be described as follows:

1. If a reaction-diffusion wavefront is detected in the neighborhood of an agent, and the agent is not on a refractory cell then the agent moves towards the wavefront.

2. If nothing is detected, then the agent may initiate a reaction-diffusion wave with a fixed probability λ , the firing rate.

Characteristic Results from the CA Implementation of the Model

Quantitative results regarding the aggregation characteristics have been presented in [START_REF] Vlassopoulos | How fast can Virtual Amoebae aggregate? Analysis for the Optimal Firing Rate in an Instance of the Reaction-Diffusion-Chemotaxis Aggregation Scheme[END_REF]. The most promising observations deal with the robustness of the aggregation scheme in the presence of obstacles in the environment as well as the the presence of perturbations in the movement of the agents. The ability of the virtual agents to aggregate in spite of various perturbations is mainly due to the properties of the GHCA. The reaction-diffusion wavefronts are capable of traveling around obstacles and regenerate when they travel through "holes" in obstacles. Figure 3 shows an example of aggregation in a rather constrained environment, consisting of rooms separated by doors. To study the behavior with respect to perturbations in the agent movement, the original model introduced a "thermal motion" factor, so that each agent would jump to a neighboring cell with some probability p a . Although this type of perturbation is quite different from what we would expect to observe in the application of the scheme to robot agents, still, it gives us a reason to believe that the model can be applied to robots. It is interesting to mention that the aggregation of the agents persisted even when the probability p a was quite high (up to 10%). Figure 4 shows an example for p a =2%. All these simulations results were rather encouraging ; the challenge is now to see if this robustness can somehow be "transposed" in the case of robotic agents. t=0 t=200 t=400 t=600

Figure 4: Example of the aggregation of the virtual amoebae in the presence of noise in the agent movement

FROM VIRTUAL AMOEBAE TO ROBOT AGENTS

The main differences between the original CA-based instance of the aggregation scheme and its robot counterpart concern space and time. In simulations, both the agents and the environment operate in discrete time but in the physical experiments the movements of the robots are continuous. Moreover, in the discrete model, the actual size of the agents is neglected and each agent is simply assigned to a cell location. In the physical experiment, the agents have a nonnegligible size and they cover many environment cells. As a consequence, the relative length of the environment cells with respect to the agent dimensions become important, since they affect how the agents perceive and interact with the environment. The discretization time, that is, the time unit which separates two consecutive CA steps is also a key parameter. The other practical problem to solve is how to make the agents aware of the environment and how to enable them to interact with it. This section describes the solutions we studied for solving these problems.

Chemotaxic Behavior

Our starting point for designing the robots behavior was that we needed to adapt to the abilities of the ALICE robots (see Figure 5 for a picture). These robots are equipped with two wheels; their speed and direction is controlled by setting the speed of each wheel. Two light sensors are located on the top of the robot. These sensors are connected to analog-to-digital converters and are subject to high noise ratio. With such robots, the main problem in applying the aggregation scheme is related to a problem of directional ambiguity: when detecting the waves with only two sensors, the robot can not distinguish between two incoming wavefronts with the same angle. This ambiguity is illustrated on Fig. 6: as the only available information is the different time of arrivals of the wavefronts, the agent can only determine if the wave is more or less parallel to the line that connects the two sensors. In other words, in the best case we are able to determine the absolute value of the angle of the incoming waves, but not its sign. The strategy we have chosen to overcome the directional ambiguity is to design a special movement for the robots. We programmed the robots to turn for an angle of 90 O towards the direction of the incident wave, that is, towards the first sensor that detects an incoming wave (Fig. 7). As a consequence, the robot will eventually turn towards the source of the waves, although this might require that more waves reach the robot.

To understand why we did not want to follow the gradient of the incident reaction-diffusion waves, consdier the following case. Imagine that a reaction-diffusion wave, whose intensity gradually varies from a maximum value to zero is first detected by one of the sensors and after a moment from the other sensor. Assume that the agent eventually turns so that its direction is perpendicular to the incoming wavefronts, that is, it is now well positioned to go towards the source. Since there are only two sensors and as they are placed perpendicular to the movement direction of the agent, they are now both parallel to the wavefronts. As a consequence, they become unable to accurately detect the gradient, since the light intensity varies in the same way on the two sensors, as it shown in Figure 8. This problem appears both when the robot faces towards or away from the reaction-diffusion waves. To avoid this second ambiguity, we programmed the robots to only react to the presence of wavefronts, that appear as a sudden change in the detected luminosity. When such a change is detected, the robot turns towards the sensor that detected the change. This rules allows a simple calibration and control logic for the sensors. More specifically, the decision on whether the robot should turn left or right is made only by calculating the difference of the measurements of the left and right sensor and comparing this difference to a threshold T L .

The effect of these rules is that the agents follow a ``zig-zag'' track towards the source of the waves (see Fig. 9).

Figure 9: "Zig-zag" approach for moving towards the source of reaction-diffusion waves

Adapting the Agent Behavior

One other problem with had to deal with was to ensure that the robots would not be affected by the reaction-diffusion waves that they initiated. This issue stems from the fact that the robot dimensions span several cells. An excitation that starts from a cell in the robot and extends outwards will eventually reach the sensors, and therefore the robot will "misunderstand" its own emission as an incident wave. Although we used a different approach to resolve this issue for the ALICE and KheperaIII robots, the most simple modification to the above model is to insert a refractory period after the initiation of a wave. During this refractory period the robot ignores all incoming waves.

The rules that govern each robot is shown in Fig. 10. Note that this figure also displays the agitation movements as well as the obstacle avoidance, which we will explain next.

The finite-state machine shown above covers the most complex cases that have been implemented in the KheperaIII robots. For the much simpler ALICE robots, where obstacle avoidance, agitatory moves and firing were not implemented on the robots, the state machine consisted of only 3 states, i.e. IDLE, MOVE LEFT and MOVE RIGHT.

Experimental Setup for the environment

As we have already seen, ALICE robots are equipped with two light sensors that are located on top of the robot. This arrangement was the basic guideline on the setup of the aggregation experiments, since the reaction-diffusion waves should be projected on the surface where the robots are located, from above. Following an idea by G. Theraulaz, we used a setup that was based on a camera and a projector: the images from the camera were processed in real-time in order to identify the positions of the robots in the environment and to map them to cell locations in the CA array. Knowing the agent locations and the state of the array, we simulated the evolution of the environment using FiatLux CA simulator2 and projected the resulting images back on the surface where the agents were located. Figure 11 shows a picture of one experiment, where we can see the agents and the projected waves.

Interaction between robots and environment

The main problem of this setup is that although the agents can "read" their environment and detect incoming wavefronts, it is not possible for them to "write" directly on the environment, that is,. to initiate a wave. In order to overcome this problem, we delegated the wave initiation part of the agent logic to the environment. However, since the fire decisions are not taken by the agents, it becomes impossible for them to know whether they should enter a refractory state and ignore their own emitted waves. The solution to this problem was again implemented in the environment, and more precisely by generating a "mask" that prevented from painting the first steps of the evolution of the reaction-diffusion process. In the case of the KheperaIII robots, setting up the experimental environment was greatly simplified, since the interactive ROMEA table and the existing programming interfaces provided all the tools required to track the position of the robots and simulate the reaction-diffusion process (see Fig. 12). The only drawback compared to the ALICE experiment was that we needed to work in total darkness because of a special tracking device that uses infrared leds that are situated below the KheperaIII robots.

Robot Parameters

The final set of the parameters that characterize the model varied significantly in the two implementations. In the simplest case of the ALICE robots, the only things that the robot needed to know were the speed and movement duration for each wheel, as well as the value of the threshold, T L . In the case of the KheperaIII robots, where the fire decisions were not merged with the environment logic, the parameter set included the fire probability, as well as the parameters for detecting other robots and the edges of the table. The experiment uses a "circular" neighborhood.

A quick tuning of the parameters

It is important to note that setting up the aggregation scheme for Alice or KheperaIII robots was done by tuning only a small number of parameters. More precisely, the main parameters that characterize the agents are the firing rate λ and the refractory time, the turning angle and movement time constants. Regarding the environment, the main parameters are the number of states M and the discretization time, that is, the time between two CA evolutions of the environment.

The tuning of these parameters were that the following balance was respected :

1. In order to speed up the gathering process, both the robots and the environment should evolve as quickly as possible. 2. In order not to mix up informations sent by each robot, the dynamics of the environment, that is, the speed of the reaction-diffusion waves, should stay higher than the robots speed. This in particular guarantees that no wavefronts will be detected more than once. 3. The environment should not evolve too rapidly in order for the robots to have a sufficient time to detect the waves.

To tune the turning angle, we implemented a simplified model in software to emulate as closely as possible the ALICE robots, with configurable sensor and wheel positions as well as turning radius. Using this model we ran several simulations with various turning angles and also trying out cases where the angle was subject to noise. In all cases, the robots managed to converge, with the fastest average convergence time obtained for a fixed turning angle of 90 O .

EXPERIMENTS AND OBSERVATIONS

Readers are encouraged to consult the AMYBIA project web-page 3 , which displays videos of our ALICE and KheperaIII experiments. Some experiments are also accessible through an applet that replays the recorded robot positions and orientation.

Clear environment with Alice and KheperaIII Robots

Our first set of experiments was carried out with ALICE robots in an environment without obstacles. In such conditions, we could observe that the aggregation always succeeded, despite numerous errors in the robots' movements. For instance, although the robots that we used for our first experiments (ALICE) were calibrated so as to turn by 90 O when detecting a wavefront, their real turning angle would vary from 0 to sometimes much over 90 O . However, this did not obstruct the aggregation process and in all experiments we performed, the robots managed to gather in a tight area. Figure 13 presents an example of an aggregation process with four robots. As it can be seen, the four trajectories are convergent but are also very noisy. From an observer point of view, the aggregation also looks slow : with our settings, the time needed to make the robots come in a close area typically varied from 1 to 5 minutes. The key point is that these drawbacks should not be interpreted as a default of our system: it is precisely the use of randomness in the triggering of the waves and in the movements of the robots that guarantees the system's properties of robustness. As observed from the simulations we carried out, it should also be understood that it is this "clever" use of randomness that will allow the system to scale and be applied to a large team or to a swarm of robots. Indeed, all happened as if the system "wants" to minimize the distance between robots and each time step. From an observer point of view, one sees that the "enveloppe" that contains all the robots progressively shrinks until robots get close to another.

The only case were the aggregation process was blocked was when spurious sources of emissions of waves existed. In our case, such errors resulted from "false" detections but one can also imagine that they come from robots that will always be emitting. We observed that these errors caused no problem when they were were rare, but in the case where a "ghost source" continuously triggered waves, the aggregation process was significantly perturbed. Indeed, the robots could not separate the consecutive waves which brought them to turn randomly.

Obstacle Avoidance with KheperaIII Robots

Our next step was to consider environments with physical obstacles that we added to the environment. This experiments were made to study how the gathering process is affected by a change which, according to our assumptions, the agents are not aware of. From the perspective of the CA environment that propagates the reaction-diffusion waves, permanent obstacles were represented as cells with a fixed state.

From the robot's perspective, we used three strategies to cope with such obstacles. The common denominator of these strategies was that whenever an obstacle was present in the path of a robot, the robot would immediately stop moving, whatever its current state.

1) The first and simplest method was to add random moves to the agents. Such moves would occur when the agent was idle, that is, not performing a chemotaxis move or being in a refractory state. They consisted of a turn towards a random direction followed by moving for a fixed distance towards this direction. The use of this strategy increased the aggregation time, but performed good in the presence of both permanent and non-permanent obstacles.

2) The second strategy was to make a robot move backwards on a given distance when it meets an obstacle. We observed that this method has a good performance in the case where there are no permanent obstacles but that it created deadlocks in some configurations with permanent obstacles. Such deadlocks occur when the size of the obstacle is greater than the "backtracking" distance of the robot. Figure 15 illustrates a deadlock scenario. First, the robot detects a wave and turns towards the obstacle; it then backtracks and now faces the obstacle. When the next wave is detected, it turns left and approaches the obstacle. When the robot detects the obstacle, it backtracks and this time faces it left. The whole process repeats again and again, introducing a deadlock. The third strategy we used was a slight improvement on the backtracking method. Instead of "blindly" moving backwards, and since it was easy to determine which side of the robot is not blocked by the obstacle, the robot would move towards the "most free" area. Surprisingly, this method performed even worse than the simple backtracking and increased the deadlocks instead of resolving them. The reason for such a behavior was that the robots were always attracted in the same area and had less possibilities to escape the deadlock by pure chance.

3) The most promising strategy we employed was a variation on the first method which consisted of a combination of two random moves. First, the agitatory movement was preserved. Second, whenever an obstacle was blocking a robot and the robot wanted to move towards that direction it would perform a sequence of several agitation movements, increasing the chances to move towards a free area. This strategy resolved many deadlocks in the process and gave the best results in terms of a trade-off between minimizing the aggregation tie and keeping the individual behavior of the robots simple (see Fig. 14 for an example of trajectory).

There are of course many other possible methods for avoiding deadlocks. For instance, the robots follow the gradient of the reaction-diffusion waves instead of just detecting the wavefronts. However, as already mentioned, this would imply a more complex and possibly less robust behavior of the system.

In short, the two basic features of the amoebae aggregation scheme is that it works using only local operations and stigmergy, as well as the fact that it delegates most of the communication overhead to the environment. The operations are local both for the agents and for the CA-based environment. This significantly simplifies the agent design and agent logic, but also requires the existence of active cells to relay the messages.

CONCLUSIONS AND FUTURE WORK

We presented our first results on the application of a decentralized aggregation scheme on two different types of robots. Our method to achieve the aggregation used a stigmergic principle and employed only local interactions of the robots with their environment. The experiments carried out allowed us to observe that the aggregation scheme is simple to implement and resists various perturbations on the robots. The main advantage of our method is to use only minimal assumptions on the robots abilities and on their sensors (number, precision, fidelity, etc.). The main difficulty for implementing the aggregation scheme is the design of an active environment to relay the messages without attenuation. In our experiments, the environment was simulated by an external (centralized) computation and projected on the robots. We also considered the case where this environment contained permanent obstacles and observed that, contrarily to the simplified discrete model, we need to add extra feature to the robots in order to solve some "deadlock" situations.

So far, our experiments involved no more than a few robots. Nevertheless, from the CA simulations we have observed that the model maintains good aggregation properties even for high number of agents (up to a few hundreds). These aggregation properties are maintained through the change of scale mainly because of the bio-inspired nature of the scheme and its "intelligent" use of randomness (no need for electing a "leader" is involved). In the case were a swarm of robots is used, the challenge would be to keep a "reasonable" aggregation time. This goal could be achieved by minimizing the number of colliding waves that annihilate. Indeed, not all the emitted waves reach all agents and therefore this somehow slows down the aggregation. From a previous research on the aggregation scheme [START_REF] Vlassopoulos | How fast can Virtual Amoebae aggregate? Analysis for the Optimal Firing Rate in an Instance of the Reaction-Diffusion-Chemotaxis Aggregation Scheme[END_REF], we know that there is an optimal value of the firing rate. It is thus an open question to know how this optimal value depends on the number of robots and on other parameters.

As a part of our ongoing research, we are developing a model that simulates the KheperaIII -ROMEA interaction. The model used allows us to test and verify our observations and try new algorithms for obstacle avoidance and deadlock prevention. In order to gain some realism, we intend to develop intermediate layers to run simulations on larger or more complex environments. Our goal is to investigate whether the good scaling properties observed with the discrete CA simulation model are kept with these more realistic models.

Finally, we would like to study how the system behavior varies as we change the assumptions made on the agents abilities. Indeed, increasing the "level of cooperation" of the agents with the environment could provide a more efficient gathering process, at least in the case where obstacles are present. This raises the question of how to set the trade-off between simplifications made in the agents and the robustness and efficiency of the aggregation process.

Figure 1 :

 1 Figure 1: Cell rules as a state transition diagram

Figure 2 :

 2 Figure 2: Reaction-diffusion process with a lattice size L=25, M=4 and5 initially excited cells. Excited cells are in bright red, refractory cells are in a red and grey mixture, neutral cells are in white.

Figure 3 :

 3 Figure 3: Agent aggregation in the presence of obstacles, for times t=0, 20, 100 and 700. Agents are shown in green, obstacles in blue, neutral cells in white and excited cells with red.

Figure 5 :

 5 Figure 5: ALICE robot

Figure 6 :Figure 7 :

 67 Figure 6: An example of the directional ambiguity, where the robot cannot distinguish between the two incident wavefronts

Figure 8 :

 8 Figure 8: Example of the gradient ambiguity: the sensor will see the same temporal variations for a wave moving leftwards and a wave moving to the right

Figure 10 :

 10 Figure 10: The complete state machine of the robot agents

Figure 11 :

 11 Figure 11: Snapshot of an experiment with ALICE robots, as taken from the camera that detects the robot positions.

  This obstacle detection feature constitutes the major difference between the ALICE and KheperaIII experiments. Indeed, although the absence of obstacle avoidance in the first set of experiments resulted in the ALICE robots blocking and pushing each other, this behavior would only delay the aggregation. The robots would eventually unblock and succed to aggregate in a tight area. On the other hand, with the KheperaIII robots, using obstacle avoidance is mandatory since the robots should neither push each other, nor hit on the walls of the ROMEA table.

Figure 12 :

 12 Figure 12: Snapshot of an experiment using the Khepera robots on the interactive ROMEA table.The experiment uses a "circular" neighborhood.

Figure 13 :

 13 Figure 13: Recorded trajectories for the aggregation of 4 robots.

Figure 14 :

 14 Figure 14: Recorded trajectories for the obstacle avoidance algorithm that uses multiple agitation moves (strategy 3)

Figure 15 :

 15 Figure 15: Example of deadlocks introduced by the "backtracking" obstacle avoidance method

Designed at EPFL, Lausanne, Switzerland, 1998, see: http://en.wikipedia.org/wiki/Alice_mobile_robot

http://fiatlux.loria.fr

http://amybia.loria.fr/


Arnaud Badel

Blandine Rozier

Kohki Takahashi

Satoshi Awaji

Simulation of Local Dissipation Phenomena in the REBCO Insert of the 25-T CSM Magnet: Understanding and Preventing Destructive Thermal Runaway

Keywords: HTS magnets, quench propagation, very high field and nmr magnets. I

We recently developed an electrothermal model to study the occurrence and evolution of dissipative zones in a Rare-Earth BaCuO (REBCO) high-temperature superconductor (HTS) coils, considering the local critical current (I c ) inhomogeneity with a centimeter-scale resolution. In this study, we present an application of this model, implemented in Comsol Multiphysics, to the REBCO insert of the 25-T CSM magnet developed in the High Field Laboratory for Superconducting Material of Tohoku University. The destructive thermal runaway event or "quench" observed on that coil will be introduced, followed by the modeling approach that was used to model it. We could simulate the coil behavior with a single 2-cm-long defect and verify that a more sensitive detection system in the millivolt range could be used to safely protect such magnets even from such a local defect. In the case of a more usual inhomogeneity distribution of I c , the temperature rise is more homogenous and safe protection can be reached even at higher detection levels.

R ARE-EARTH BaCuO (REBCO) High Temperature Su- perconductors (HTS) are attractive for very high field inserts due to their large engineering current densities under high magnetic fields when operating close to 4.2 K, combined with their good mechanical strength. That is the reason why it was the first choice for designing the HTS insert of the 25 T Cryogen-free Superconducting Magnet (CSM) developed at the High Field Laboratory for Superconducting Material [START_REF] Awaji | New 25 T cryogen-free superconducting magnet project at Tohoku university[END_REF]. Specific developments were conducted on REBCO tapes produced by Fujikura for that project. They could typically carry around 400 A/cm-width under 20 T transverse field at 20 K and sustain more than 780 MPa without degradation [START_REF] Fujita | Characteristics of REBCO coated conductors for 25 T cryogen-free superconducting magnet[END_REF]. The initial design A. Badel was with CNRS and Université Grenoble Alpes, Grenoble 38042, France. He is now with Tohoku University, Sendai 980-8577, Japan (e-mail:, arnaud.badel@imr.tohoku.ac.jp).

B. Rozier is with the Université Grenoble Alpes, Grenoble 38042, France (e-mail:,blandine.rozier@g2elab.grenoble.inp.fr).

K. Takahashi and S. Awaji are with Tohoku University, Sendai 980-8577, Japan (e-mail:,kohki@imr.tohoku.ac.jp; awaji@imr.tohoku.ac.jp).

TABLE I 25 T CSM REBCO INSERT CHARACTERISTICS

consisted in a 14 T Low Temperature Superconducting (LTS) outsert and an 11.5 T REBCO insert cooled down independently, respectively to around 5 K and 7 K, using a conduction cooling system detailed in [START_REF] Takahashi | Design and test results of a cryogenic cooling system for a 25-T cryogen-free superconducting magnet[END_REF]. The characteristics of the REBCO insert are summarized in Table I.

Commercial REBCO tapes display performance inhomogeneity along their lengths, with mm-scale variations of the local critical current measured continuously at 77 K in self-field conditions. Combined with their very good thermal stability, it leads to a high risk of damaging hot spots when operating close to the estimated critical current, as in that case the critical current may be overstepped locally on the weaker zones, leading to local damaging hotspots.

In the case of the 25 T CSM REBCO insert, the choice was made to use a very large (more than 50%) current margin [START_REF] Miyazaki | Design of a REBCO insert coil for a cryogen-free 25-T superconducting magnet[END_REF] in order to limit the risk of such local critical current overstepping. In principle, the magnet could operate as high as 20 K. Moreover, Fig. 1. Top, evolution of detection bridge voltage (Balance V) and coil current up to the thermal runaway event (from [START_REF] Awaji | First performance test of a 25 T cryogen-free superconducting magnet[END_REF]). Bottom, current evolution during discharge (interpolated from experiment data).

all the REBCO tape lengths used in the magnet were controlled against unusual local defects using continuous measurement under 77 K self-field [START_REF]TAPESTAR-Family: Quality Management for Coated Conductors[END_REF]. They typically showed maximum local variation of ±7.5 % around the Ic 77 Ks.f mean value. Then, each of the pancakes were tested individually at 77 K to ascertain that the conductor was not damaged during the process [START_REF] Awaji | 10 T generation by an epoxy impregnated GdBCO insert coil for the 25 T-cryogen-free superconducting magnet[END_REF]. In operation, both the outsert and the insert were protected by classical bridge detection system triggering discharges in dump resistors. The bridge threshold voltage in the case of the REBCO insert was set to 1-2 V, with a delay of 0.5 s. That protection system was mostly meant to protect the REBCO insert in case of quench induced by the quench of the outsert. The REBCO insert was tested successfully in standalone, but a destructive thermal runaway occurred during the first ramp up under background field, at 124.6 A (94% of the nominal value) [START_REF] Awaji | First performance test of a 25 T cryogen-free superconducting magnet[END_REF]. The field reached at that moment 24 T. Later on it was replaced by an 11 T BSCCO insert [START_REF] Hanai | Development of an 11 T BSCCO insert coil for a 25 T cryogen-free superconducting magnet[END_REF], with which the 25 T CSM is now operational [START_REF] Awaji | First performance test of a 25 T cryogen-free superconducting magnet[END_REF].

The evolution of voltage and current around the quench event is shown Fig. 1. The objective of this work is to use numerical modelling in order to study the possible origins and the dynamic of that event, and to evaluate how this behaviour would have been improved using a more sensitive detection system. The results of the simulation study will be presented and discussed in part III, after the modeling approach is introduced in part II.

II. MODELLING APPROACH

REBCO coated conductor are thin tapes with multilayered structure. In the case of Fujikura's tape for the 25 T CSM insert, the tape is 5 mm wide, with a 2.5 µm REBCO layer grown on a 75 µm-thick Hastelloy C276 substrate through a number of thin buffer layers of about 0.5 µm in total. It is surrounded by 2 µm of Silver and 20 µm of Copper (Fig. 2). Both the thermal and electrical properties of REBCO tapes are strongly non-linear. Together with the high aspect ratio of the layers, it would make the solving of the local transient coupled thermal-electric behaviour everywhere in the magnet very time-consuming, if not impossible in the general case. For this reason, we use homogenization as well domain study reduction in our model as will be explained in the next paragraphs. The model is implemented in Comsol Multiphysics as a transient heat transfer problem. The electric model is contained in external functions that are called for each blocks before each time step to evaluate the local heat generation. We will introduce briefly here the model hypothesis approach, the modelling implementation is presented in more details in [START_REF] Badel | Modelling of 'Quench' or the occurrence and propagation of dissipative zones in REBCO HTS coils[END_REF].

A. Electrical Model

From the electrical point of view, we homogenize the tape at the cross-section level but discretize its length. This is coherent with the available characterization data for the tape critical current that is given at a given position on the length but for the whole tape. Fig. 3 show an example of typical in-length characterization data for Fujikura's tape used in the REBCO insert.

As in our case the magnet is wound with single tape conductor, this results in an equivalent circuit consisting in a string of seriesconnected non-linear variable resistances.

For each elements of the 1D in-length discretization of the electrical problem (called "blocks" in the following) the local critical current must be evaluated in order to compute the equivalent resistance. This local critical current depends from the local tape performances and the operation conditions in terms of field (amplitude and orientation) and temperature. We first evaluated a "nominal" distribution of Ic based on short sample characterization, then modulated it by a "tape performance" coefficient varying around unity in a range consistent with the spread observed in continous Ic 77 K s.f measurement. The exact variation of Ic in a coil cannot be obtained (for existing coils) or foreseen (for future ones). Our aim by modulating the "expected" Ic with distributed defects is simply to give to the model an input Ic distribution that is plausible, and study its influence on the overall behaviour. Fig. 4 display an exemple of Fujikura carcaterization data [START_REF] Fujita | In-field critical current property of IBAD/PLD coated conductors[END_REF] used as reference values.

Once the critical current for each block is obtained, the equivalent resistance per block is evaluated taking into account the current sharing with the other layers, using iterative methods.

B. Thermal Model Definition

Considering the pancake strcuture of the winding, we can limit the domain to 2D and study one pancake independently from the others. Knowing from preliminary studies where the weakest Ic is most likely to be found, the thermal model geometry can be defined.

For the numerical modelling of the transient thermal problem, we homogenize the tape thermal properties. The layered structure leads to a strong anisotropic thermal conductivity that must be considered. The co-wound insulation is a 55 µm meterthick teflon-coated polyimide tape. It is a distinct layer between each turn in the model, in order to keep the possibility to vary its thickness and to apply thermally resistive interfaces between the turns.

Finally, the 25 T CSM insert pancakes were actually impreganted but the impregantion material was mostly kept out of the winding thanks to the winding tension. We neglect it in our model.

III. APPLICATION TO THE 25 T CSM REBCO INSERT

A. Magnetostatic Study and Domain Reduction

2D axi-symmetric magnetostatic studies are conducted to get the field contributions of both the outsert and the insert for a given combination of insert and outsert current. It is used to to evaluate the "expected" local Jc on every turn of every pancakes (Fig. 5), which are then averaged on the turn width to obtain the local Ic. In the Fujikura tapes used in the REBCO insert, the variation of Ic at 77 K s.f. was typically below ±7.5 %. A study of all the turns in which the expected Ic is at most 7.5 % higher than the lowest expected Ic (i.e., lower than 157 A) should include the area in which the risk of dissipation occurrence is maximum.

In our study, it represents most of the top pancake (from turn 125 to 400) and does not concern the others. However, the total tape length to be considered is still more than 190 m, which would generate a model too heavy to handle. In that study we limited ourselves to 100 turns, centered arround the weakest one, that is turn 270. In the eletrical model the conductor length is discretized in about 3000 elements. For each of these elements, which are approximately 2 cm long, an Ic value was thus set based on the turns expected Ic modulated by a random coefficent between 0.925 and 1.075.

The initial operating temperature is 6.5 K, obtained from measurement on the real insert. To represent the conduction cooling in the 2D model, an out-of-plane heat exchange of 2 W/m².K is considered on both pancake flat surfaces from a cold plate at 6 K. This value is derived from the available cooling power of the system at that temperature.

B. Simulation Procedure

In the model study the background field is considered already established at nominal value (14 T in the center). In order to shorten the simulation time, the insert current is rised at a constant ramp rate of 0.1 A/s, while in reality it was closer to 0.017 A/s. Simulations were conducted on few cases at the two speeds to verify that such ramp rate change had no significant impact. The total voltage across the studied domain is evaluated, and dicharge of the current is initiated after it reaches a given detection threshold. The current discharge profile after the quench protection trigger is known from experiment (see Fig. 1), so the simulation was made with driven current, the preregistered current profile being imposed during the discharge. This way were able to circumvent the limitations of the implementation of our model in Comsol that does not allow easily the coupling of electrical circuits. We are currently developping a more flexible implementation of the model [START_REF] Badel | Modelling of 'Quench' or the occurrence and propagation of dissipative zones in REBCO HTS coils[END_REF] that includes such coupling, but for now it cannot manage a geometry with the dimensions considered here.

C. Results and Discussion

The 25 T CSM REBCO insert was damaged during operation at current lower than the nominal, which was already defined with a large margin. The most likely reason was the existance of a sharp local defect. We could reproduce a thermal runway of the magnet at the observed current (around 125 A) with a single 2 cm-long defect having only 18.8 % of the expected Ic. The resulting voltage and hotspot temperature evolution are shown Fig. 6 for various detection threshold.

We first used the protection setting used in the experiment, 1 V. In that case we could observe that thermal runaway occurs but does not destroy the coil immediately. However, the model evaluate the hotpost temperature at 600°C (destruction) 2 s only after the discharge is trigered, instead of more than 6 s in experiment. This may be explained by the fact that our model represents only a quarter of a pancake volume as only 100 turns are modelled out of 435. In case of large heat generation, this limits the spread of the heat diffusion and leads to an overestimation of the peek temperature.

Then we studied the possibilities offered by early detection in such experimental case, varying the threshold value (Fig. 6). The thermal runaway can be prevented even considering the rather slow discharge of the insert and its coupling with the outsert, provided that the detection threshold is lower than 20 mV. This value was obtained with no protection delay. A trigger around 12-18 mV would make it possible to wait more than a second (10 s on the graph corresponds to a 1 A variation before the discharge) to confirm the detection event. In order to minimize fatigue due to local temperature rise, the trigger should be set between 5 and 10 mV so that the temperature is kept below 50 K.

However, these studies are made in the worst possible scenario, that of a single sharp defect. For comparison purpose we then studied the same insert, but with Ic properties homogeneously reduced to 23 % of their expected value, so that the thermal runaway occurs at the same current. The dissipation is more homogenous (Fig. 7). A 20 mV threshold value leads to a maximum temperature of only 17 K.

IV. CONCLUSION

We recently developped a model to study the occurrence and evolution of dissipative zones in REBCO HTS magnets and the resulting thermal runway phenomenon, taking into account the local inhomogeneity of REBCO tapes performances. This model was used to investigate the destructive thermal runway that occurred in the REBCO insert of the 25 T CSM magnet developed at the HFLSM of Tohoku University. Although imperfectly, we could reproduce such event and demonstrate that a more sensitive detection system, with threshold value in the mV range, could have been used to protect such insert effectively from sharp local defects, while detection system in the tens of mV range would be sufficient in case of a more typical distribution of inhomogeneity.

Fig. 2 .

 2 Fig. 2. Cross section of Fujikura tape used in the REBCO insert under study.

Fig. 3 .

 3 Fig.3. Typical Ic distribution (lower curve) of 5 mm-wide Fujikura tape used in the project, at 77 K self field (taken from[START_REF] Fujita | Characteristics of REBCO coated conductors for 25 T cryogen-free superconducting magnet[END_REF]).

Fig. 4 .

 4 Fig. 4. Variation of Ic per Unit width with Field in Transverse Direction at Various Temperature (from [8]).

Fig. 5 .

 5 Fig. 5. Top, Jc distribution in the REBCO insert based on magnetostatic study of the 25 T CSM in nominal operation. Bottom, Turn by turn variation of Ic in the top pancakes of the insert under nominal operation. lowest curve is top pancake, second one is 2 nd from the top and so on.

Fig. 6 .

 6 Fig. 6. Dissipating voltage (top) and hotspot temperature (bottom) evolution for various protection discharge thresholds. Origin of time was set arbitrarily.

Fig. 7 .

 7 Fig. 7. Temperature distribution at the beginning of the protection discharge for a 20 mV trigger and ± 7.5% distributed conductor inhomogeneity.