Jan Dentler
email: jan.dentler@uni.lu

Martin Rosalie
email: martin.rosalie@uni.lu

Grégoire Danoy
email: gregoire.danoy@uni.lu

Pascal Bouvry
email: pascal.bourvy@uni.lu

Somasundar Kannan
email: somasundar.kannan@uni.lu

Miguel A Olivares-Mendez
email: miguel.olivaresmendez@uni.lu

Holger Voos
email: holger.voos@uni.lu

Collision Avoidance Effects on the Mobility of a UAV Swarm Using Chaotic Ant Colony with Model Predictive Control

come

Introduction

Over the last decade, the usage of UAVs has been strictly increasing in areas like the construction industry, agriculture, defense and security. The advantages of UAVs are their remote sensing capabilities as well as their mobility and related extended operational space, which are ideally suited for area exploration tasks. For instance, the response time in search and rescue scenarios within large areas can be reduced significantly by using multiple autonomously operating UAVs. However, the coverage of the surveyed area is dependent on the coordination of the UAVs. A comprehensive overview on mechanisms to coordinate and control such swarms are given in [START_REF] Gazi | Swarm Coordination and Control Problems[END_REF] and [START_REF] Bertrand | MPC Strategies for Cooperative Guidance of Autonomous Vehicles[END_REF]. The detailed state of the art in multi-UAV area exploration from a perspective of optimization and artificial intelligence is summarized in Fig. 1 Use case: Implementation of area exploration scenario. CACOC gives UAV waypoints to the MPC which controls UAV models in order to simulate a realistic closed-loop system. [START_REF] Nigam | The multiple unmanned air vehicle persistent surveillance problem: A review[END_REF]. A typical approach which is covered in all of these contributions is a cooperative area exploration control with potential fields over a receding horizon. The nature inspired interpretation of this method has been first introduced in [START_REF] Gaudiano | Swarm intelligence: a new C2 paradigm with an application to control of swarms of UAVs[END_REF] as Ant Colony Optimization with pheromones to coordinate a swarm of UAVs. This deterministic path planning is differing from Dynamic Data Driven Application Systems (DDDAS, see [START_REF] Mccune | Swarm control of UAVs for cooperative hunting with DDDAS[END_REF] for an example applied on a swarm of UAVs), where the UAVs coordinate according to the data they collected.

Kuiper and Nadjm-Tehrani [START_REF] Kuiper | Mobility models for UAV group reconnaissance applications[END_REF] define an area exploration scenario with ten UAVs using a pheromone based approach for path planning. Their study compares a random mobility model to a pheromone-based mobility model using metrics for coverage. Some recent studies propose to replace the random part of such metaheuristics by chaotic dynamics to increase their performance, e.g. for Particle Swarm [START_REF] Gandomi | Chaos-enhanced accelerated particle swarm optimization[END_REF][START_REF] Pluhacek | Chaos particle swarm optimization with eensemble of chaotic systems[END_REF] or Firefly Optimization [START_REF] Gandomi | Firefly algorithm with chaos[END_REF]. Thus, previous work [START_REF] Rosalie | From random process to chaotic behavior in swarms of UAVs[END_REF] focused on the usage of chaotic dynamics to create a Chaotic Ant Colony Optimization to Coverage (CACOC) mobility model in order to increase the area coverage performance of a standard Ant Colony algorithm. However, the high-level simulations used to evaluate CACOC do not take into account real quadrotor dynamics. For this reason in [START_REF] Rosalie | Area exploration with a swarm of UAVs combining deterministic chaotic ant colony mobility with position MPC[END_REF] CACOC has been evaluated in an area exploration scenario using ten realistic UAV models, which have been introduced in the authors previous work [START_REF] Dentler | A real-time model predictive position control with collision avoidance for commercial low-cost quadrotors[END_REF]. However, in [START_REF] Rosalie | Area exploration with a swarm of UAVs combining deterministic chaotic ant colony mobility with position MPC[END_REF] the UAVs are operating on different altitudes in order to avoid collisions. In reality this is not always feasible due to airspace and physical limitations and furthermore leads to inhomogeneous sensor perception area sizes and ground resolutions.

The main contribution of this work is addressing this limitation by extending the CA-COC mobility model by a CA (CACOC+CA) mechanism. In addition, the paper analyses the performance of the proposed approach in a realistic use case exploration scenario with three UAVs, as illustrated in Fig. 1. The realistic UAV behavior is modelled by using an identified model of real quadrotor dynamics in combination with a MPC position control, as introduced in the authors previous publication [START_REF] Dentler | A real-time model predictive position control with collision avoidance for commercial low-cost quadrotors[END_REF]. Furthermore, the in [START_REF] Dentler | A real-time model predictive position control with collision avoidance for commercial low-cost quadrotors[END_REF] presented CA technique for static obstacles is applied within this work to avoid collision between UAVs (MPC+CA). The result are realistic closed-loop UAV dynamics which reflect the position tracking with MPC of waypoints that are generated with the CACOC mobility model. Based on this realistic UAV behavior, this work is evaluating the efficiency of the proposed CA-COC+CA with MPC+CA approach by means of coverage and CA metrics. For this purpose, new metrics are introduced to evaluate the CA capabilities of the utilized models in longterm simulations. These simulation serve as preliminary analysis steps before the deploy-ment of real UAVs is possible. Additionally, this paper provides more details about chaotic dynamics and the CACOC mobility model. This article is organized as follows. In section 2 the description of the CACOC UAV mobility model along with a CA mechanism is provided. This mobility model includes chaotic dynamics to obtain unpredictable and structured trajectories. Furthermore, the related state of the art and implementation details are given. Section 3 describes the MPC of UAVs including a CA mechanism as well as the underlying dynamic model of a UAV. The experimental scenario combining CACOC and MPC is given in section 4. This includes an evaluation of the proposed methods and a description of the metrics utilized for the statistical analysis. Section 5 is concluding the presented work and outlines future research directions.

UAV mobility models

This section first provides a brief description of the CACOC mobility model used to generate the target waypoints for the UAVs which lead to an optimal area exploration. This mobility model has been introduced in [START_REF] Rosalie | From random process to chaotic behavior in swarms of UAVs[END_REF][START_REF] Rosalie | Area exploration with a swarm of UAVs combining deterministic chaotic ant colony mobility with position MPC[END_REF]. Its algorithm and main properties are detailed within this section. Subsequently, we associate a theoretical CA method to the CACOC mobility model to enhance the capability of the swarm in order to allow the UAVs to fly at the same altitude.

Chaotic Ant Colony Optimization to Coverage (CACOC)

Typical applications of UAV swarms are area surveillance, search and rescue, etc. The advantage of swarms in such scenarios is the distributed sensing capability, which is further enhanced by the movements of each individual UAV. In order to fully exploit these advantages, the UAV movement has to be coordinated. A traditional method to optimize the area coverage is pheromone based Ant Colony Optimization (ACO). In order to improve the performance of this approach, previous work [START_REF] Rosalie | From random process to chaotic behavior in swarms of UAVs[END_REF] has combined ACO with the chaotic behavior of a dynamical system resulting in the CACOC mobility model. More precisely, CACOC replaces the random parts of the ant colony optimization algorithm [START_REF] Dorigo | Optimization, learning and natural algorithms[END_REF] by a chaotic behavior in such a way that the exploration capabilities of the UAV swarm is improved. The paper presented a high level simulation process in which CACOC has been used to describe the mobility of a swarm of ten UAVs. Its performance was compared to existing models using random, chaos and classical ant colony optimization algorithms [START_REF] Rosalie | From random process to chaotic behavior in swarms of UAVs[END_REF]. The results demonstrate empirically that CACOC permits an increase in performance for several aspects of coverage (overall coverage, recent coverage and fairness of the coverage). The CACOC performance on the coverage problem was empirically demonstrated using metrics as described in Sec. [START_REF] Danoy | Connectivity stability in autonomous multi-level UAV swarms for wide area monitoring[END_REF].

This section provides a brief description of this mobility model and its extension regarding CA. For more details the interested reader should refer to [START_REF] Rosalie | From random process to chaotic behavior in swarms of UAVs[END_REF]. We will first present a mobility model including chaotic dynamics (CROMM) before introducing CACOC that combines both chaotic dynamics with an Ant Colony exploration algorithm. One approach to generate chaotic dynamics is to use the Rössler system [START_REF] Rössler | An equation for continuous chaos[END_REF]. It is a well-known dynamical system that can exhibit various chaotic mechanisms [START_REF] Rosalie | Templates and subtemplates of Rössler attractors from a bifurcation diagram[END_REF].

The UAVs are considered to have a constant speed and choose a movement direction at each discrete time step : A for ahead, R for 45 • right and L for 45 • left. Then, the next choice is given by the first return map (Fig. 2). This map underlines the dynamical signature of the Rössler system giving ρ n+1 versus ρ n . Thus, the next action depends on the previous one:

if ρ n < 1/3 then direction is right (R); if 1/3 ρ n < 2/3 then direction is left (L); else the direction is ahead (A). Algorithm 1 CROMM mobility model This basic mobility model is named CROMM and its pseudo-code is detailed in Alg. 1. In that case, the good exploration performance of the UAVs is due to the periodic orbits of the system that lead to patterns (Fig. 2). The periodic orbits with low periods are considered as skeleton of the chaotic dynamics and are often visited during the simulation process using a Table 1 Pheromone action table for ACO UAV mobility model. le f t is the amount of pheromones sensed on the left of the UAV, ahead is the amount of pheromones sensed in front of the UAV and right is the amount of pheromones sensed on the right of the UAV; total = le f t + ahead + right.

Probability of action with pheromones Left

Ahead Right

p L = total-le f t 2×total p A = total-ahead 2×total p R = total-right 2×total
Runge-Kutta (4 th order) algorithm. As illustrated Fig. 2, the period 1 orbits leads to symbols AAAAA. . . (straight line); the period 2 leads to the pattern ARARA. . . (large circle) and the period 4 leads to the serpentine pattern RALARALA. . . For CACOC, we used an ant colony algorithm combined with chaotic dynamics. If there is no virtual pheromone to guide the UAV (pheromones are deposited by each UAV to indicate areas they already visited), CROMM is used. If there are pheromones, they have repulsive properties and the next choice of direction depends on the total amount of pheromones sensed around the UAV. Such a repulsive effect prevents the UAVs to reach recently visited areas since pheromones evaporate in time. We also used the first return map values (Fig. 2) to choose the next direction with the pheromones' perception instead of a random number. Consider that p L , p A and p R are inversely proportional to the total amount of pheromones sensed respectively to the left, ahead and right of the UAV and that p R + p L + p A = 1 (Tab. 1). Thus, with ρ n taken from the first return map (Fig. 2), the next direction is chosen according to these rules:

-if ρ n < p R next direction is right; -if p R ρ n < p R + p L next direction is left; -else the direction is ahead. The pseudo-code of CACOC is detailed in Alg. 2.
Algorithm 2 CACOC mobility model

CACOC and Collision Avoidance (CACOC+CA)

One of the drawbacks of the previously introduced CACOC mobility model is that it does not prevent UAVs from colliding. This can be avoided by flying the UAVs at different altitudes, which however might not be possible due to limited airspace and also leads to differing sensor perception areas and ground resolutions. For these reasons, it is preferable that UAVs are able to fly at the same altitude which requires CA mechanisms. Therefore, this work aims to include CA in the CACOC mobility model (CACOC+CA).

State of the Art

In the literature, several surveys detail the main methods used to implement CA for UAVs [START_REF] Albaker | A survey of collision avoidance approaches for unmanned aerial vehicles[END_REF][START_REF] Albaker | A conceptual framework and a review of conflict sensing, detection, awareness and escape maneuvering methods for UAVs[END_REF][START_REF] Mahjri | A review on collision avoidance systems for unmanned aerial vehicles[END_REF]. One of this article provides a classification of the CA system (see Fig.

3 [START_REF] Albaker | A survey of collision avoidance approaches for unmanned aerial vehicles[END_REF]). In [START_REF] Geng | Mission planning of autonomous UAVs for urban surveillance with evolutionary algorithms[END_REF], authors used evolutionary algorithms to find an optimal solution for path planing including CA constraints. However, this family of methods is not suitable for our problem because we are aiming to obtain unpredictable trajectories (without repetition of the same path). The same issue comes with methods proposing CA for static objects using a collision cone philosophy to predict the collision (see [START_REF] Tripathi | Reactive collision avoidance of UAVs with simple pin-hole camera based passive stereovision sensing[END_REF] for instance). Most of the literature uses a safety cone to detect the potential collision and initiate evasive maneuvers (see [START_REF] Leonard | Distributed reactive collision avoidance for a swarm of quadrotors[END_REF] for instance).

In 2012, Lalish & Morgansen [START_REF] Lalish | Distributed reactive collision avoidance[END_REF] introduced the Distributed Reactive Collision Avoidance (DRCA) algorithm which is a two-step process. It consists of an optimization-based deconfliction maneuver, followed by the longer-term deconfliction maintenance phase. Its reactive and force-field type approach has been recently extended for swarms of quadrotor helicopters [START_REF] Leonard | Distributed reactive collision avoidance for a swarm of quadrotors[END_REF]. Authors of [START_REF] Han | Proportional navigation-based collision avoidance for UAVs[END_REF] use the proportional navigation (PN) guidance law to build a CA system. This guidance law was originally designed for missile engagement scenarios and require persistent trajectories.

All of these algorithms handle various situations including heterogeneity and priority between UAVs in order to permit the UAVs to reach their destination close to the nominal trajectory. However, the CACOC mobility model is typically using low waypoint and direction update intervals. Accordingly, the chosen direction of a UAV is not persistent wherefore the cone evasion method or proportional navigation cannot be used efficiently. Furthermore, the goal of the cooperating UAVs is maximum area coverage which is not considered in the previously mentioned techniques. For example for several UAVs in close proximity, the goal is to spread them out by using a repelling mechanism. To the best of our knowledge, no current algorithm is available to suffice chaotic fast direction changes combined with the coverage problem capable of avoiding collisions of UAVs.

Collision Avoidance for CACOC

As previously mentioned, the CACOC mobility model is based on CROMM (Alg. 1). The latter one uses the properties of periodic orbits to obtain exploratory patterns: straight line, right turn and serpentine movements (Fig. 2). In order to avoid collision with CACOC, a collision evasion maneuver is chosen complementarily to these patterns: a large turn to the left. The evasion maneuver is initiated when one or more UAVs appear in a UAV neighborhood. This left turn is combined with the average direction of the neighbors (UAVs in the neighborhood) and lasts until the other UAVs leave this neighborhood. We use the velocities:

previous direction (p) velocity to the left (l) average velocity of the neighbors (0.8 × n) computed from p of the other UAVs ⇒ net velocity (d = l + 0.8 × n)

As illustrated by Fig. 3, the net velocity gives the direction of the UAVs. Then this direction vector is normalized to obtain the CA vector with the same norm as CACOC vectors:

d s = s l + 0.8 × n ||l + 0.8 × n|| (1)
with s the UAVs speed as a parameter of our model. For the sake of clarity, the mobility model is referred to as CACOC+CA in the remainder of this article. The algorithm of CA-COC+CA is given by Alg. 3.

Previous direction (p)

Velocity to the left (l)

Average velocity of the neighbors (0.8 × n)

Net velocity (d = l + 0.8 × n) UAV 2 UAV 1 UAV 2 UAV 3 UAV 1
Fig. 3 Velocities imposed to prevent collision: blue arrow represent the velocity to the left (l), the black arrow represents the average direction resulting from the other UAVs (n) and the red arrow is the direction given (net velocity d = l + 0.8 × n). This schema only represents the direction of the UAV to illustrate how CA is handled. The norm is not proportional to the movement executed by the UAVs (see Fig. 4 for the trajectories of the UAVs) and Eq.( 1) for the normalization. move according to the current state Fig. 4 shows two consecutive maneuvers of CA during a CACOC+CA simulation. According to the classification of the CA system (see Fig. 3 [1]), Fig. 4 is showing the distinctive evasive left turns. As a remark, large turns to the left can occur in CACOC (see Fig. 2). Hence, the large turn to the left permits to include a mechanism of CA without introducing a new complex pattern. In order to evaluate the proposed CACOC+CA approach statistically, its performance is evaluated in long-time simulations. As the goal is to deploy the UAVs at the same altitude, the CACOC+CA mobility model is expected to be slightly less efficient than a CACOC mobility model, due to the inhibiting character of CA. In the long-time simulations the computed CACOC waypoints are serving as target points for MPC controlled UAVs with CA inclusion which is subject of section 3. The consequence is a realistic evaluation scenario which allows to evaluate the performance of CACOC considering realistic UAV dynamics, controllers and CA.

UAV control

In order to describe the dynamic behavior of a real quadrotor system, previous work [START_REF] Dentler | A real-time model predictive position control with collision avoidance for commercial low-cost quadrotors[END_REF] presented an identified state space model of an AR.Drone 2.0 quadrotor. This state space description has been used as prediction model within a model predictive controller which has been applied to a real AR.Drone 2.0 system. The resulting control performance validated the quality of the identified model. For this reason, the same model is utilized for the evaluation of the presented CACOC+CA approach. The utilized UAV model is illustrated in Fig. 5. The state of the i-th UAV can be described by vector

x i = x i , y i , z i ,Ψ i , v f ,i , v s,i . (2)
For means of visualization, the time dependency of states and controls (e.g., x i (t)) is not explicitly shown here. The system state (2) consists of the quadrotor's x i , y i , z i position in Cartesian world coordinates (W ), its yaw angle Ψ i and the forward v f ,i and side ward velocities v s,i in the vehicle frame (V ). The UAV inputs are concatenated in its control vector

u i = u f ,i , u s,i , u z,i , u Ψ ,i , (3)
containing forward u f ,i , side ward u s,i , upward u z,i and yaw u Ψ ,i velocity. The identified system model of [START_REF] Dentler | A real-time model predictive position control with collision avoidance for commercial low-cost quadrotors[END_REF] is given as

ẋi = f i (x i , u i ) =         v f ,i cos (Ψ i ) -v s,i sin (Ψ i ) v f ,i sin (Ψ i ) + v s,i cos (Ψ i ) u z,i 1.6u Ψ ,i -0.5092v f ,i + 1.458u f ,i -0.5092v s,i + 1.458u s,i         . ( 4
)
This model is based on a second order model for the x and y-channel and single integrator systems for z and Ψ . This assumption of single integrator systems in z and Ψ implies that the target velocity in this channel is reached instantaneously. This is a suboptimal estimation, but does not influence the evaluation of CACOC as the UAVs are flying on the same altitude and do not change their orientation. As further remark, the AR.Drone 2.0 quadrotor behavior shows nonlinearities in the z-axis, due to unexpected behavior in the inherent sensing solution combining cameras and ultrasonic-sensors. However, except the nonlinear mapping of the forward and sideward velocities to ẋ and ẏ no further nonlinearities are considered by the model. Despite, the model is sufficiently precise to serve as prediction model for a real AR.Drone 2.0. The scalar parameters in (4) have been identified from a real AR.Drone 2.0 system by means of a motion capture system. As commercial drones are typically operated by velocity commands, a position controller is required to reach a desired position. [START_REF] Dentler | A real-time model predictive position control with collision avoidance for commercial low-cost quadrotors[END_REF] is using a fast MPC for this purpose which is based on the computation of optimal controls by solving an optimal control problem (OCP). In [START_REF] Rosalie | Area exploration with a swarm of UAVs combining deterministic chaotic ant colony mobility with position MPC[END_REF], this approach has been applied to ten UAVs using OCP min

u i J i = t f t 0 (x * i -x i ) Q (x * i -x i ) + u i Ru i dτ (5)
s.t. ẋi = f i (x i , u i ) (6)
c ≤ u f (t) 2 -1 u s (t) 2 -1 u z (t) 2 -1 u 2 Ψ -1 (7)
x i (0) = 50, 0, i, 0, 0, 0

with Q = D{ 1.5, 1.5, 1.6, 0.1, 0, 0 },

(9) R = D{ 1, 1, 1, 1 }. ( 10
)
The idea of MPC is to predict the UAVs behavior within a horizon τ = [t 0 ,t f ] by means of the system dynamics [START_REF] Dentler | Implementation and validation of an event-based real-time nonlinear model predictive control framework with ROS interface for single and multi-robot systems[END_REF] which are equal to the model (4). Based on this prediction, the control inputs u i are determined in order to minimize cost functional J i while satisfying constraints [START_REF] Dentler | A real-time model predictive position control with collision avoidance for commercial low-cost quadrotors[END_REF] and initial conditions [START_REF] Dorigo | Optimization, learning and natural algorithms[END_REF]. The position tracking in J i is achieved by a quadratic penalty (9) of the state error, while energy optimality is imposed by a quadratic penalty (10) of the controls.

The advantage of MPC is its ability to directly consider constraints in the OCP (e.g The AR.Drone 2.0 input limitations -1 ≤ u ≤ 1 in ( 7)) and a scenario independent control policy. In [START_REF] Dentler | A real-time model predictive position control with collision avoidance for commercial low-cost quadrotors[END_REF] this advantage has been used to avoid the collision of a UAV with a static obstacle by imposing the constraint as additional cost term

l CA (x i , x j , d min ) = a 1 + e -b d 2 min -(x i -x j ) (xi-xj) (11)
in J i . The idea is to impose costs of height a, if the Euclidean distance between the drone positions x i and x j is smaller than distance d min . For larger drone distances the costs vanish to 0. This switching behavior is approximated by a sigmoid with the steepness parameter b. The result is a repulsive cost-gradient around the violation distance. The disadvantage of this constraint handling is its cost trade-off with the trajectory tracking and energy optimality. Unfortunately, this cost gradient cannot be chosen arbitrarily steep, as it worsens the OCP condition. Therefore, the cost trade-off is also a trade-off between repulsion from the violation and solution feasibility, which does provide CA only for limited tracking errors and system velocities. Within this work the CA distance is set chosen to d min = 1.5m.

As l CA is considering the positions of two UAVs, the cost functions cannot be treated separately for each UAV as in ( 5)- [START_REF] Gandomi | Chaos-enhanced accelerated particle swarm optimization[END_REF]. As shown in previous work [START_REF] Dentler | A modularization approach for nonlinear model predictive control of distributed fast systems[END_REF] for distributed systems and in [START_REF] Dentler | Implementation and validation of an event-based real-time nonlinear model predictive control framework with ROS interface for single and multi-robot systems[END_REF] for UAVs, the CA constraints [START_REF] Gaudiano | Swarm intelligence: a new C2 paradigm with an application to control of swarms of UAVs[END_REF] are coupling the individual OCPs. Combining UAV dynamics, MPC for coupled multi agent systems and CA cost functions for three UAVs (i = 0, 1, 2) results in OCP min

u 0 ,u 1 ,u 2 J = t f t 0 2 ∑ i=0 (x * i -x i ) Q (x * i -x i ) + u i Ru i +l CA (x 0 , x 1 , 1.5) + l CA (x 0 , x 2 , 1.5) + l CA (x 1 , x 2 , 1.5) dτ (12) with l CA (x i , x j , d min ) = 1.0 1 + e -3.0 d 2 min -(x i -x j ) (xi-xj) . ( 13
) s.t. ẋi = f i (x i , u i ) , i = 0, 1, 2 (14)
c i ≤ u f , i (t) 2 -1 u s , i (t) 2 -1 u z , i (t) 2 -1 u 2 Ψ ,i -1 , i = 0, 1, 2 (15)
x 0 (0) = 0, 0, 1, 0, 0, 0 (16) x 1 (0) = 0, 15, 1, 0, 0, 0 (17) x 2 (0) = 15, 0, 1, 0, 0, 0 (18) Q = D{ 1.5, 1.5, 1.6, 0.1, 0, 0 },

(19) R = D{ 1, 1, 1, 1 }. ( 20
)
For means of visibility of the area coverage, the desired trajectory waypoints in the given scenario are limited to changes in the xy-plane, while the initial z i = 1 and ψ i = 0 are tracked. The resulting closed loop behavior is referred to as MPC+CA within this work.

Analog to previous work [START_REF] Dentler | A real-time model predictive position control with collision avoidance for commercial low-cost quadrotors[END_REF], a Condensed Multiple Shooting Continuation Generalized Minimal Residual (CMSCGMRES) [START_REF] Shimizu | A real-time algorithm for nonlinear receding horizon control using multiple shooting and continuation/krylov method[END_REF] approach is used to solve the MPC with OCP ( 12)- [START_REF] Nigam | The multiple unmanned air vehicle persistent surveillance problem: A review[END_REF]. The CMSCGMRES is parameterized with a maximum number of iterations k max = 6, horizon length T MPC = 1s, control update interval of ∆t MPC = 0.1s, n hor = 10 steps within the horizon, a forward difference step of h = 0.001s, solution tolerance of ε = 10 -4 , continuation factor ζ = 10 and the horizon expansion factor α = 2. More details about the parameters are given in [START_REF] Ohtsuka | A continuation/gmres method for fast computation of nonlinear receding horizon control[END_REF].

The resulting position tracking and CA is shown in Fig. 6. MPC+CA with inherent UAV model serves as realistic UAV simulator for the CACOC+CA evaluation which is subject of the next section 4.

Experiments

This section introduces the simulation framework and experimental setup, including all parameters. Subsequently, all metrics for the analysis of the experiment are detailed. Finally, the statistical results of the analysis are given. For the experimental evaluation of CACOC, respectively CACOC+CA, MPC+CA and three AR.Drone 2.0 UAVs V-REP [START_REF] Rohmer | V-rep: a versatile and scalable robot simulation framework[END_REF] models are implemented according to Fig. 7 using ROS [START_REF] Quigley | ROS: an open-source robot operating system[END_REF] as communication interface. The CACOC waypoints are updated every second ∆t CACOC = 1s, wherefore one CACOC step equals to ∆t sim = 1s of simulation. The UAV positions in V-REP are updated every ∆t V REP = 0.1s using the dynamics (4). As the controller is synchronized with V-REP, the control update interval is equally ∆t MPC = 0.1s. The simulation environment is illustrated in Fig. 8 1 . The simulation is performed using the parameters shown in Tab. 2.

Metrics

For the evaluation of the experiment the same three area coverage metrics as introduced in [START_REF] Rosalie | Area exploration with a swarm of UAVs combining deterministic chaotic ant colony mobility with position MPC[END_REF] are used. A short description is provided hereinafter. Fig. 9 underlines how the metrics are calculated for one simulation. The three UAVs start in the bottom left corner with the following coordinate on a grid: (0, 0), (0, 15) and (15, 0).

Coverage

The coverage metric is the percentage of cells of the total area visited during the whole simulation. The environment is a 30 m by 30 m square area. The coverage value varies during the whole simulation. To have a representative value of the coverage, we choose to compare the coverage value after 4000 steps for each model. This indicates the efficiency of the models to visit the total area. On the other hand, we want to evaluate the first steps of each model to compare their initial behavior. For this purpose the first 500 steps are used to extracted the slope of a linear regression a × x (Fig. 9).

Fairness of the coverage

The fairness measures, if all cells are regularly and equally scanned. This is computed as the standard deviation (STD) of their respective number of scans [START_REF] Schleich | UAV fleet area coverage with network connectivity constraint[END_REF]. To evaluate the fairness during the whole simulation, a linear regression a × x + b is performed using the last 3500 steps (Fig. 9). This measure is complementary to the coverage initial slope that only evaluates the initial UAV trajectories. Here, the slope value is considered as a measure to evaluate the fairness of the models without this initialization part.

Recent coverage ratio

This metric introduced in [START_REF] Danoy | Connectivity stability in autonomous multi-level UAV swarms for wide area monitoring[END_REF] represents the percentage of coverage during the last 100 iterations. These 100 steps correspond to the pheromones' evaporation time. For the computation of the recent coverage mean value (Fig. 9), the 100 first iterations are excluded from the simulation.

Distance between each couple of UAVs

The metrics used to quantify the efficiency of the CA mechanism are dedicated to run where only one collision is scheduled on purpose [START_REF] Leonard | Distributed reactive collision avoidance for a swarm of quadrotors[END_REF]. However, here the goal is to cover a wide area where multiple collisions can occur. As a consequence, new metrics are necessary to Fig. 10 For a given simulation, using CACOC+CA with MPC+CA, the distance between each couple of UAV is shown in meters. An estimation of the distribution is given with the mean value as well as the minimum.

(a) CACOC with MPC + CA (b) CACOC + CA with MPC + CA Fig. 11 For a given simulation, using CACOC+CA with MPC+CA, the distance between the desired position and the controlled one is shown in meters. An estimation of the distribution is given with the mean value as well as the minimum.

quantify the CA mechanism efficiency on a long-time simulation. An intuitive way of evaluating the CA mechanism is to measure the distance between each couple of UAVs during the whole simulation. This includes the mean value and the minimal value of this metric to make a comparison between CACOC+CA with MPC+CA and CACOC with MPC+CA. For a given simulation, Fig. 10 represents the estimated distribution of the measured values.

Distance between the desired and the current position

In order to quantify few CA occurrences in a simulation, the authors of [START_REF] Leonard | Distributed reactive collision avoidance for a swarm of quadrotors[END_REF] use the excess separation deviation from the path and the graph status (indicating if the UAV is doing CA). We introduce a new metric by computing the distance between the waypoint of a UAV and the position given by the MPC before receiving the new waypoint. During the whole simulation it indicates how far the UAV is from the current waypoint. We collect the minimum value and the mean value of this metric to quantify the CA.

For the CACOC+CA with MPC+CA example simulation, Fig. 11(b) shows that only few of the UAV2 values are far from the median value. This indicates that despite the CA included in CACOC, the CA of the MPC is used. For the CACOC with MPC+CA example simulation in Fig. 11(a), only the CA of the MPC is used resulting in several points far from the median value for all UAVs. We collect both mean value and minimum values to have an average value of the distance combined with the shorter distance reached during the simulation. As a remark, the static offset in the tracking error is expected, due to the constantly changing target position and missing integral part in the MPC control policy.

Experimental results

The previously described metrics are obtained with the same program to avoid bias. First, they are computed from the waypoints and secondly using the synchronized positions of the UAVs from the V-Rep simulation using ROS. As CACOC is a chaotic deterministic mobility model, only the initial conditions of the dynamical system are varied to create a wide range of simulations with the same mobility model without changing its parametrization as it has already been done previously [START_REF] Rosalie | From random process to chaotic behavior in swarms of UAVs[END_REF].

Comparison of the coverage metrics

Table 3 details the coverage metrics statistics to compare CACOC with MPC+CA and CA-COC+CA with MPC+CA. Data for coverage after 4000 steps follow the same distribution with 95% confidence according to a Wilcoxon test (p-value = 0.1357). This permits to confirm that for both, CACOC with MPC+CA and CACOC+CA with MPC+CA, we fulfill our objective by covering almost the whole area: for both the median value is 1 meaning that 100% of the area has been covered in more than half of the simulations. Fig. 12 shows that the slope of coverage values of CACOC+CA with MPC+CA values are lower than the values of CACOC with MPC+CA. Data for slope of coverage follow a different distribution with 95% confidence according to a Wilcoxon test (p-value = 0.0002189). This small difference is due to the introduction of the CA mechanism that prevents a very good spreading of the UAVs at the beginning of the simulation compared to the CACOC mobility model. However, this metric mainly quantifies the early steps of the algorithm. Hence, the performance reduction is not significant as the starting points of the UAVs are not in close proximity. We remind that the lower the fairness, the better the performance. Data for fairness follow different distributions with 95% confidence according to a Wilcoxon test (p-value = 0.04146). The p-value is close to the threshold value to accept the hypothesis, meaning that the results are similar. Fig. 13 details how CACOC+CA with MPC+CA is slightly less efficient than CACOC with MPC+CA. Although the two resulting distributions are not the same, the difference between the mean values is only 3%. As this metric evaluates the global performance of the algorithm, we can conclude that the usage of a CA mechanism does not significantly impact the coverage of the area.

Data for recent coverage follow different distribution with 95% confidence according to a Wilcoxon test (p-value < 2.2e -16) and Fig. 14 details the recent coverage values. In our case, for 3 UAVs over a 30 × 30 area with pheromones' duration of 100 steps, the optimal value is 3 × 100/(30 × 30) = 1 3 . The mean value of the recent coverage of CACOC+CA with MPC+CA values is 0.22% of the optimal value and accordingly lower than the average value of CACOC with MPC+CA. This marginal difference implies that even with a CA mechanism, the mobility model is efficient. Considering the overall coverage performance of the CACOC+CA with MPC+CA compared to CACOC with MPC+CA, we can say that even if the values are not following the same distribution, they are close to the value without CA. This small offset is reasonable comparing to the added security by CA. In the next section, we will compare the CA metrics to see the difference between the two methods.

Comparison of the collision avoidance metrics

The Figs. 15-18 illustrate the value distribution for each metric dedicated to the CA. We thus perform a pairwise non-parametric test (using the Wilcoxon signed-rank test): where we observe that the data follows different distributions with 95% confidence. Tab. 4 synthesizes the statistics for these CA metrics. For the average distance between the actual position of a UAV and its waypoint (Fig. 15), the average values are similar, but the standard deviation (STD) is almost doubled for the CACOC+CA with MPC+CA. This implies that CACOC+CA with MPC+CA with the lower mean performs slightly better, but the variation of its tracking performance is higher. This variation is explainable with a border lock phenomenon occurring when UAVs come close to border and CA constraints. As the UAVs should stay in the the given search area they are forced back into it. As a result, the UAV trajectories start to oscillate and stay consequently close together for a while. We remind that the CA of the MPC is rarely used when there is already a CA in the mobility model (see Fig. 11a where a significant amount of outliers is visible). Considering the minimum distance between the UAV and its waypoint (Fig. 16) the distributions are not the same but the range of values is similar. Taking into account the mean values (0.06 for CACOC with MPC+CA and 0.11 for CACOC+CA with MPC+CA in Tab. 4) we can consider that in average, the second method does not permit two UAVs to get very close to waypoints. We explain this shift because of the CA mechanism in CACOC with MPC + CA that can cause the UAVs to avoid each other and reach positions far or close to the given waypoints. This happens less frequently in CACOC+CA with MPC+CA because the second CA mechanism (MPC+CA) is typically not used. The distribution of the average distance between each couple of UAVs (Fig. 17) is totally different and shows outliers. These outliers are also a result of the border lock phenomenon, inducing oscillations around the area border constraints. They can furthermore be found in the minimum value plot of CACOC+CA with MPC+CA in Fig. 17. However, without this specific behavior, the average distance is higher by 15%. This implies that the UAVs are more often far from each other than width CACOC with MPC+CA. Finally, considering the minimum distance between each couple of UAVs (Fig. 18), the average value of CACOC+CA with MPC is twice as high as for CACOC with MPC+CA for CACOC+CA and reaches a mean value of more than 1.5m. Thus, CACOC+CA complements the CA implemented in MPC to reach the desired minimum distance of d min = 1.5m. The higher variance can be explained by the border lock phenomenon and is subject to further research.

Conclusion

To support the real deployment of UAVs for coverage purposes, we are aiming to design robust and resilient systems. The paper is extending the Chaotic Ant Colony Optimization to Coverage mobility model by a CA mechanism (CACOC+CA). Its coverage and CA performance is compared to the mobility model without CA (CACOC) in physical simulations within a dedicated robotic environment. The physical simulation includes a dynamic UAV model controlled with MPC including a CA mechanism (MPC+CA). The empirical evaluation of both "CACOC with MPC+CA" and "CACOC+CA with MPC+CA" was conducted on an area coverage scenario with 3 UAVs. In this context, three area coverage metrics and four CA metrics have been used to evaluate and compare the algorithm's performance. The obtained results show that CACOC+CA is able to accomplish covering tasks efficiently while complementing the potential field CA mechanism of the MPC controller, reducing the average tracking error of UAVs and increasing the mean minimal distance between UAVs close to the chosen MPC CA distance d = 1.5m. The variance of the minimal distance is assumed to be a result of the border-lock phenomenon and subject to further research. In terms of coverage, the high similarity of the results between CACOC with MPC+CA and CACOC+CA with MPC+CA demonstrates that it performs well in a real environment under consideration of a fine-tuning of the MPC approach. Thus, the implementation of CA with CACOC is advantageous to deal with collisions of UAVs, while the CA of MPC can be used to deal with unknown obstacles (eg. birds).

In order to transfer the proposed approach to real UAVs in future, the CA bounds have to be adjusted according to the localization precision. For scenarios with more UAVs, taking into consideration a CA constraint for every UAV pair does not scale well. This issue can be addressed by using a distributed MPC controller which only considers CA constraints to direct neighbors. Further work on the UAV control part will treat the reduction of the position tracking error, for instance, with the introduction of a disturbance model, the adaptation of the MPC control law, or a target position controller as already presented in [START_REF] Dentler | A real-time model predictive position control with collision avoidance for commercial low-cost quadrotors[END_REF]. Furthermore, more research work will be dedicated to address the border lock phenomenon. The presented two-sided simulator will also permit to work on the CACOC parameters (pheromone evaporation time, covered area by step, etc.) to increase the performance of the mobility model by comparing the theoretical and realistic results. Finally, we also planed to develop a timeefficient software to handle the MPC+CA algorithm for a swarm of 10 and 20 UAVs.

Fig. 2

 2 Fig. 2 Details of the CACOC mobility model with the first return maps from the Rössler system and the pattern it generates when there is no pheromone to guide the UAVs. The next choice depends on the previous one with: A for ahead, L for left and R for right. The period 1 orbits leads to symbols AAAAA. . . (strong violet straight line). The period 2 leads to the pattern ARARA. . . (lime green large circle). The period 4 highlighted in gray leads to the serpentine pattern RALARALA. . . (lime green and strong violet lines).

Fig. 4

 4 Fig. 4 Trajectories of the three UAVs doing maneuver to avoid collision.

Fig. 5

 5 Fig. 5 AR.Drone 2.0 Quadrotor with utilized coordinate conventions.

Fig. 6

 6 Fig. 6 MPC+CA: MPC position tracking for UAVs (AR.Drone 2.0) with CA. Target positions marked green, UAV positions red. The sequence is showing two intersecting UAV targets and the resulting evasion maneuver imposed by the MPC+CA.

Fig. 7

 7 Fig. 7 Implementation setup with ROS as communication interface between CACOC, MPC and V-REP

Fig. 8

 8 Fig. 8 Three MPC controlled UAVs in the V-Rep simulation environment following waypoints provided by CACOC.

Fig. 9

 9 Fig.9Metrics measurements for one simulation. Two values are extracted for the coverage: the coverage after 4000 steps and the slope of coverage for the 500 first steps. For the fairness the slope of fairness is computed after the 500 first steps. The average value of the recent coverage is calculated after the 100 first steps.

Fig. 12

 12 Fig.12Slope of the coverage (see Fig.9for details about the linear regression).

Fig. 13

 13 Fig.13Fairness slope of the coverage (see Fig.9for details about the linear regression).

Fig. 14

 14 Fig.[START_REF] Han | Proportional navigation-based collision avoidance for UAVs[END_REF] Recent coverage (see Fig.9for details about the average value).

Fig. 15

 15 Fig.15Average distance in meters between the UAV and its waypoint.

Fig. 16

 16 Fig.16Minimum distance in meters between the UAV and its waypoint.

Fig. 17

 17 Fig. 17Average distance in meters between each couple of UAVs.

  Fig. 17Average distance in meters between each couple of UAVs.

Fig. 18

 18 Fig. 18 Minimal distance in meters between each couple of UAVs.

Table 2

 2 Main experimental parameters.

	Parameter Name	Parameter Value
	Simulation area
	Geographical Area	30 m × 30 m
	Number of cells	30 × 30
	Simulation Setup
	Update interval of CACOC waypoints ∆t Experiments
	Mobility models	CACOC with MPC+CA, CACOC+CA with MPC+CA
	Number of UAVs	3
	Simulation steps	4000
	Independent runs	30

CACOC = 1s Update interval of UAV dynamics in V-Rep ∆tVREP = 0.1s Control update interval ∆t MPC = 0.1s Time duration Tsim = 4000s for 4000 steps UAV Waypoint UAVs speed 1 m/s Possible UAV actions ahead, 45 • left, 45 • right Initial UAVs position middle of the bottom of the map

Table 3

 3 Coverage metrics statistics. For each metric, the first line details CACOC with MPC+CA and the second line details CACOC+CA with MPC+CA and the best mean values are highlighted.

	Metrics	CA mean	min	median	max	STD
	Coverage after 4000 steps	0.99967 0.99944	0.99778 0.99778	1 1	1 1	0.00066 0.00070
	Slope of coverage	0.00188 0.00183	0.00179 0.00173	0.00187 0.00183	0.00196 0.00195	3.8377e-05 4.2101e-05
	Slope of the Fairness of the coverage	8.0783e-4 7.2143e-4 8.0077e-4 9.1150e-4 8.3092e-4 7.3900e-4 8.2478e-4 9.2471e-4	4.7068e-05 4.4969e-05
	Recent coverage	0.26730 0.26051	0.26606 0.25781	0.26727 0.26046	0.26894 0.26315	9.0553e-4 1.3482e-3

Table 4

 4 CA metrics statistics. For each metric, the first line details CACOC with MPC+CA and the second line details CACOC+CA with MPC+CA and the best mean values are highlighted.

	Metrics	CA mean	min	median	max	STD
	Average distance between UAV and waypoint [m]	0.50657 0.50572	0.45554 0.41288	0.49252 0.47357	0.71484 0.90022	0.04770 0.09822
	Minimal distance between UAV and waypoint [m]	0.05954 0.10562	0.00334 0.01600	0.05376 0.11084	0.15001 0.18513	0.02980 0.04506
	Average distance between couple of UAVs [m]	15.327 17.192	14.112 10.636	15.339 17.462	16.395 18.567	0.49876 1.2416
	Minimal distance between couple of UAVs [m]	0.82445 1.5546	0.26817 0.06670	0.79422 1.6684	1.3867 2.8743	0.24397 0.67708

A video illustrating the framework simulator is available at https://martinrosalie.gforge. uni.lu/_downloads/jint.mp4

Acknowledgements This work is supported by FNR "Fonds National de la Recherche" (Luxembourg) through AFR "Aides à la Formation-Recherche" Ph.D. grant scheme No. 9312118. The experiments presented in this paper were carried out using the HPC facilities of the University of Luxembourg [32] (see http://hpc.uni.lu).


Hiroki Takahashi

Sekai Iwama

Simon Clevers

Stéphane Veesler

Gérard Coquerel

Hirohito Tsue

Rui Tamura

In situ Observation of Polymorphic Transition during Crystallization of Organic Compounds Showing Preferential Enrichment by Means of Temperature-Controlled Video-Microscopy and Time-Resolved X-Ray Powder Diffraction

A solvent-assisted solid-to-solid polymorphic transition, which was assumed to be one of the key processes for the occurrence of preferential enrichment (PE), has been proved to occur by means of temperature-controlled video-microscopy (TCVM) and time-resolved in situ X-ray powder diffraction (XRPD) measurement during crystallization of two typical first-generation of chiral racemic compounds (1 and 2) showing an excellent PE phenomenon.

INTRODUCTION

Preferential enrichment (PE) is a unique spontaneous enantiomeric resolution phenomenon applicable to racemic crystals; by repeating simple recrystallization of a certain kind of racemic crystals under nonequilibrium conditions at high supersaturation, alternating enrichment of two enantiomers occurs largely in the mother liquor and slightly in the deposited crystals with full reproducibility (Figure 1a). 1,2 This phenomenon, which contrasts with the well-known preferential crystallization of racemic conglomerates composed of a mixture of homochiral R and S crystals, [START_REF] Jacques | Racemates and Resolution[END_REF][START_REF] Coquerel | [END_REF] originates from regular chiral fluctuation or symmetry breaking after the polymorphic transition of an incipient metastable polymorphic form into a thermodynamically more stable form during crystallization from the highly supersaturated solution of a certain kind of racemic crystals (Figure S1). [5][6][7][8][9][10] Thus, the mechanism of PE has been interpreted in terms of a nonlinear complexity phenomenon including two unique processes: a solvent assisted solid-to-solid polymorphic transition and subsequent selective redissolution into the mother liquor of the excess one enantiomer from the transformed disordered crystals. [START_REF] Coquerel | Disordered Pharmaceutical Materials[END_REF][START_REF] Tamura | [END_REF][13][14][START_REF] Tamura | Enantiomer Separation: Fundamentals and Practical Methods[END_REF][START_REF] Iwama | [END_REF][17][18][19][20][21] With respect to the second process, the kinetic and thermodynamic origin has been unveiled in terms of the unique nonlinear solubility properties of the two enantiomers after the occurrence of polymorphic transition under PE crystallization conditions using high supersaturation, together with the ternary phase diagram that is consistent with the mechanism of PE. 22,23 As to the first process, we have judged the occurrence of polymorphic transition during crystallization indirectly by the in situ ATR-FTIR spectroscopy which can detect a change in the hydrogen bonding mode in small particles (< ca. 1m ) as well as both solute and solvent molecules. By comparing the IR spectrum just before crystallization with that after crystallization, we predicted the occurrence of polymorphic transition. 11- 21 However, some compounds showing an excellent PE phenomenon failed to give an appreciable IR spectral change during crystallization, due most likely to a small change in the hydrogen bonding mode upon the polymorphic transition. Therefore, the use of other in situ techniques to observe polymorphic transition directly during crystallization was indispensable. Here we report the successful direct observation of polymorphic transition during crystallization by the combined use of temperature-controlled video-microscopy (TCVM) 24,25 and time-resolved in situ XRPD measurement 26 under PE conditions with regard to typical first-generation of organic compounds, 1 and 2 (Figure 1b), for which the PE behavior, and the crystal structure and physicochemical properties of the deposited stable  form crystal were fully characterized. 5,6

EXPERIMENTAL SECTION

General. Melting points were measured by differential scanning calorimetry (DSC) at the scanning rate of 5°C/min.

The in situ Fourier transform infrared (FTIR) spectra were recorded in solution and suspension by using the attenuated total reflection (ATR) method on ReactIR 45m. XRPD patterns were recorded at a continuous scanning rate of 1° 2θ/min using Cu Kα radiation (40 kV, 40 mA) with the intensity of the diffracted X-ray being collected at intervals of 0.02° 2θ. A graphite monochromator was used to remove Cu Kβ radiation.

Experimental setup of TCVM. The homemade apparatus (ANACRISMAT) used is shown in Figure 4. We used standard HPLC glass vials as an observation cell containing 1 mL of solvent. The cell was inserted into the temperature controlled unit thermostatted by Peltier elements and observed by a microscope (Nikon Eclipse TE2000-U). Sequential image acquisitions were performed from minutes to hours.

Experimental setup of in situ XRPD measurement. In situ XRPD data were collected using a homemade prototype diffractometer (Figure S3). 27 This apparatus has an original goniometer with an inverted geometry (-/-). Its association with a dedicated reactor, with a transparent bottom to X-rays, provides the possibility to identify solid(s) in suspension, all along a crystallization process (time and temperature dependent), without any sampling. The detector is a LYNXEYE TM (Bruker, Germany) and the beam (Ni filtered) comes from an X-ray tube with a copper anticathode. The in situ XRPD measurement was carried out at a 9 min interval for the first 225 min and at a 240 min interval from 225 to 1300 min.

Preparation of samples. In the case of TCVM: Racemic 1 (100 mg) was dissolved in a mixture of EtOH and 2-PrOH (1/1 v/v, 1 mL) by heating to give a 10 fold supersaturated solution (versus the solubility at 20°C), while 2 (64.9 mg) was dissolved in EtOH (1 mL) by heating to give a 20 fold supersaturated solution. The resulting solution was filtrated through a 0.45 m membrane filter to remove any invisible microcrystals and pour into the observation cell. The cell was set in the temperature control unit preheated at 60°C. The cell temperature was cooled and kept at 20°C until crystallization began.

RESULTS AND DISCUSSION

Crystal structures and melting point diagrams. Racemic crystals of both 1 and 2 are classified as a "mixed crystal" (i.e., solid solution) composed of two enantiomers. Their stable racemic crystal structures are almost identical and are called a form (Figure S2). Only difference between them is that one hydroxy group attached to the asymmetric carbon atom is orientationally disordered for racemic 1 (occupancy factor: 67/33), whereas such a disorder was not noted for racemic 2. The form crystals are most commonly found for the first-generation compounds showing PE. It was proposed that polymorphic transition of the initially formed metastable -form into the stable -form was essential for the occurrence of PE (Figure S2), because the -form crystal structure was very often observed as a stable form for analogous compounds, such as racemic 3 which failed to show PE. 5,[START_REF] Coquerel | Disordered Pharmaceutical Materials[END_REF][START_REF] Tamura | [END_REF][13][14][START_REF] Tamura | Enantiomer Separation: Fundamentals and Practical Methods[END_REF] As an example to support this hypothesis, we observed that the polymorphic transition of a form into another form suddenly occurred to induce PE and since then the -form was never obtained. 7 The melting point diagrams of 1 and 2 were previously constructed by differential scanning calorimetry (DSC) analysis (Figure 2). 5,6 The measured samples were prepared in a kinetically controlled manner so that the obtained crystalline samples should contain both the metastable and stable forms. The diagram for 1 showed two flat convex curves (standing for two couples of liquids and solids; the one at low temperatures corresponding to the metastable phase ); at the initial stage of crystallization two endothermal peaks were observed by DSC analysis for crystals of less than 30% ee and more than 80% ee but not in the range of 40-70% ee. These results indicate that the two curves correspond to a mixed crystal composed of two enantiomers (Figure 2a). On the other hand, for 2, the two phase curves intersect at two points around 35% ee; one steep curve in the range of less than 35% ee corresponded to a non-stoichiometric racemic compound, and the other flat one in the range of more than 35% ee was assigned to a mixed crystal of two enantiomers (Figure 2b). In situ ATR-FTIR measurement. We previously proposed the mode of the -to-polymorphic transition during crystallization for racemic 2 by means of the in situ ATR-FTIR measurement (Figure 3a). 5 In the present study, we measured the in situ ATR-FTIR spectra for racemic 1 in the range of 1260-1160 cm -1 corresponding to S-O stretching vibrations during crystallization and the obtained spectra were compared with those of racemic 2. The spectra of 2 changed for the first 90 min, whereas no appreciable spectral change was noted for 1 during the whole crystallization process (Figure 3b). We suggested that polymorphic transition, which was induced by rearrangement of hydrogen bonds, occurred fast enough to be detected by the in situ IR spectroscopy for racemic 2, whereas that of racemic 1 was too slow to be detected and/or involved no drastic change around the sulfonate ion during crystallization. These results may be associated with the ordered or fairly disordered R and S molecular arrangement in the resulting deposited crystals of 2 or 1, respectively. Nevertheless, since the stable crystal structures of racemic 1 and 2 were almost identical, the polymorphic transition of 1 is most likely to be a -to- type, too. To confirm whether racemic 1 and 2 show a similar type of polymorphic transition during crystallization, we performed the TCVM and time-resolved in situ XRPD observation for their crystallization.

Temperature-controlled video-microscopy (TCVM).

By using the optical video microscope connected to a temperature controller between 0 and 80 °C (Figure 4), we could continually observe not only the crystallization process by cooling but also the dissolution behavior of the deposited crystals by heating.

In the case of 1, the racemic crystals were dissolved in a mixture of ethanol and 2-propanol (1/1 v/v) by heating to give a 10 fold supersaturated solution (vs its solubility at 20 °C). The resulting hot solution (60 ~ 70 °C) was filtered through 0.45 m membrane filter and the filtrate was poured into the observation cell pre-warmed at 60 °C (Figure 4c). Then the cell was cooled and kept at 20 °C until crystallization began. First transparentleaf-shaped crystals were deposited, and second crystals with a dendritic habit appeared by heterogeneous nucleation onto the leaf-shaped crystals and gradually grew (Figure 5a, b and Movie S1). To examine whether the first and second generation of crystals are the same polymorph or not, the suspension containing both crystal habits was gradually heated up to 70 °C (Figure 5c and t = 37 s in Movie S1). At the first stage, the leaf-shaped crystals changed from clear to opaque (the transition lasts for 2 s from t = 37 s in Movie S1), keeping its habit and size, due to a solid-state phase transition (SSPT). If it were a solution-mediated phase transition (SMPT), we should observe the dissolution of the metastable phase and nucleation of the stable phase. However, crystals keep their leaf habit and the speed of the phase transition is in agreement with the assumption that SSPT is faster than SMPT. At the second stage, all crystals with leaf and dendritic habits dissolved (t = 66 s in Movie S1). This crystallization and dissolution process was reproducible. Thus, it has been confirmed that the first formed clear leaf-shaped crystals are a metastable form and the second dendrites are a stable form. Furthermore, it is noteworthy that the initial formation of metastable leaf-shaped crystals followed by a polymorphic transition was fully reproducible and the stable dendrites never acted as the crystal seed to grow the same form.

In the case of racemic 2, a similar crystallization and dissolution process was reproducibly observed starting from a 20 fold supersaturated ethanol solution at 20 °C. As soon as the metastable tiny leaf-shaped crystals formed, the stable dendrites nucleated and grew onto the surface of the metastable crystals at 20 °C (Figure 6b, Movie S2). When the temperature was increased to 70 °C (Movie S3), we observed the dissolution of leaf-shaped crystals followed by that of dendritic crystals, similar to the case of racemic 1.

Thus, we could confirm the presence of two polymorphs for racemic 1 and 2 by TCVM observation. Both compounds showed a similar crystallization and dissolution behavior. Time-resolved in situ XRPD measurement. We used the in situ XRPD apparatus with a modified optical system to observe the deposited crystals on the bottom of the vessel (Figure 7). A double-jacketed vessel with an observation cell and a thermostat to control the solution temperature by circulating the coolant was used (Figure S3).

The preliminary time-resolved in situ XRPD measurement was carried out for racemic 1 and 2 under the PE crystallization conditions. The experimental results indicated that the deposited crystals of 1 exhibited sufficient diffraction intensities, whereas the diffraction intensities of 2 were too low due to the absence of a heavy atom in the molecule. The optimal crystalli-zation conditions for the in situ XRPD measurement for racemic 1 were determined; starting from a 11 fold supersaturated EtOH/2-PrOH (1/1 v/v) solution at 10 °C without stirring. 8) 9 min after the beginning of the experiment. Between the 4 th and 5 th scans (between 36 and 45 min, arrow 2), a polymorphic transition occurred, which should correspond to the nucleation and crystal growth of stable dendrite crystals onto the metastable leaf-shaped crystals observed by TCVM at 20 °C. Between the 6 th and 9 th XRPD patterns (arrow 3), it is difficult to conclude on the exact nature of the solid phase(s). One has to keep in mind that these XRPD patterns are collected under stagnant conditions. Therefore, strong preferential orientation effects could occur (i.e. modifications of peak intensities). Moreover, maturation of the phase can occur which could additionally change the crystallinity and the shape of the particles. After 90 minutes (10 th scan), the subsequent evolution seems to consist in a simple maturation of the final δ-phase (Figure S4a). This maturation does not show any sign of further evolution after 465 minutes (Figure S4b). After the last measurement, the deposited crystals were separated from the solution and the conventional powder XRPD measurement was carried out. The powder pattern simulated from the single crystal data (-form) of racemic 1 was identical to this experimental pattern (Figure 8). However, due to a limited number of diffraction peaks, the determination of the lattice parameters of the metastable crystal was not successful.

CONCLUSIONS

By using both TCVM and time-resolved in situ XRPD measurement, we have clearly observed the crystal growth and polymorphic transition processes during crystallization for two typical first-generation of chiral racemic organic compounds (racemic 1 and 2) showing PE. The TCVM observation of both compounds indicated that the metastable transparent-leafshaped crystals were formed first, followed by a heterogeneous nucleation of the stable form onto the surface of the first-formed crystals to give the second stable dendritic crystals. We clearly evidenced a solid-state phase transition in one case (racemic 1). This mode of polymorphic transition has turned out to differ from our speculation that the transition might occur inside each crystal. By the time-resolved in situ XRPD measurement for racemic 1, the drastic diffraction pattern change corresponding to the polymorphic transition of the metastable -form to the stable -form was observed in conformity with the TCVM observation.

By the combined use of three real-time methods (in situ ATR-IR spectroscopy, TCVM, and time-resolved in situ XRPD measurement), we can obtain various information on a crystallization process, i.e., changes in intermolecular interaction, morphology and crystal lattice. The in situ ATR-IR spectroscopy is useful in the case that relatively fast polymorphic transition occurs with a drastic rearrangement of hydrogen bonds. The TCVM is a valuable method to observe the dissolution behavior of the deposited crystals at elevated temperatures as well as the crystallization process at lower temperatures. The timeresolved in situ XRPD measurement provides a decisive proof of polymorphic transition. It is highly desired that the determination of the lattice parameters of metastable forms becomes feasible by using the in situ XRPD data.

ASSOCIATED CONTENT

Supporting Information Available. General mechanism of PE (Figure S1); schematic representation of the intermolecular interactions in the and -form crystals (Figure S2); an overview of in situ X-ray diffractometer (Figure S3); in situ XRPD patterns of 1 from 90 min to 1300 min (Figure S4); TCMV movies for crystallization of 1 (Movie S1) and 2 (Movies S2 and S3). The Supporting Information is available free of charge on the ACS Publication website at DOI:

Figure 1 .

 1 Figure 1. (a) Features of preferential enrichment (PE) of a racemic crystal in case of substantial enrichment of (R) enantiomer in the mother liquor after the first recrystallization of chiral racemic organic compounds. (b) Molecular structures of 1 and 2 showing PE, and 3 failing to show it.

Figure 2 .

 2 Figure 2. Melting-point diagrams of (a) 1 and (b) 2. White and black circles represent the temperatures of the beginning and the end of fusion of the stable polymorph, respectively, while white and black squares in panel (a) are those of the metastable polymorph. Data cited from refs 5 and 6.

Figure 3 .

 3 Figure 3. (a) In situ ATR-FTIR spectra during crystallization from a 20 fold supersaturated EtOH solution of racemic 2, cited from ref 3. Comparison of spectra at the beginning (solid line, -) and the end (dotted line, •••) of crystallization with the spectra in the solid state of racemic 2 (dashed line, ---, δ-form crystal structure) and racemic 3 (dashed and dotted line, -•-•-, γ-form crystal structure). Black arrows indicate the absorptions corresponding to the δor γ-form supramolecular structure. (b) In situ ATR-FTIR spectra during crystallization from the 10 fold supersaturated solution of racemic 1 in a mixture of EtOH and 2-PrOH (1/1 v/v). Comparison of spectra in solution (blue broken line, ---) just before crystallization, in the suspensions at the beginning (green solid line, -) and the end (red dotted line,•••) of crystallization, and in the solid-state (black dashed line, ---, -form crystal structure).

Figure 4 .

 4 Figure 4. A temperature controlled video microscope. (a) Overview, (b) a temperature control unit and (c) an observation cell.

Figure 5 .

 5 Figure 5. Time sequence showing microscope images of racemic 1 during crystallization from the 10 fold supersaturated EtOH/2-PrOH (1/1 v/v) solution at 20 °C. (a) t = 4 s in Movie S1; the metastable leaf-shaped and the stable dendrite crystals nucleated onto the leaf-shaped crystals are highlighted as red and blue squares, respectively. The left-side inset indicates the magnification of the right-side crystals. (b) t = 33 s; the stable crystals grow on the surface of the metastable crystals. (c) t = 37 s; the metastable crystals started to dissolve when the solution was heated to 70 °C. d) t = 45 s; dissolution of leaf-shaped and stable dendrite crystals. The white arrow corresponds to 50 m.

Figure 6 .

 6 Figure 6. Time sequence showing microscope images of racemic 2 during crystallization from a 25 fold supersaturated EtOH solution at 20 °C. (a) t = 9 s in Movie S2; the metastable leafshaped crystal is highlighted as a red square. (b) t = 20 s; the stable dendrite crystals nucleated on leaf-shaped crystal are highlighted as a blue circle. (c) t = 47 s; the stable dendritic crystals grown onto the surface of the metastable crystal are highlighted as a blue square . (d) t = 64 s; leaf crystals started to dissolve when the solution was heated to 70 °C. The white arrow corresponds to 50 m.

Figure 7 .

 7 Figure 7. Observation cell (double jacketed vessel) wrapped with insulating material. (a) Side view without insulating material, (b) side view with insulating material, (c) bottom view, (d) top view, and (e) overview. In this experiment, the bottom hole was covered by KAPTON® thin film (7.5 µm). With respect to the time-resolved XRPD patterns for racemic 1, major transition points are denoted by arrows 1-3 in Figure 8. First crystals were detected (arrow 1 in Figure8) 9 min after the beginning of the experiment. Between the 4 th and 5 th scans (between 36 and 45 min, arrow 2), a polymorphic transition occurred, which should correspond to the nucleation and crystal growth of stable dendrite crystals onto the metastable leaf-shaped crystals observed by TCVM at 20 °C. Between the 6 th and 9 th XRPD patterns (arrow 3), it is difficult to conclude on the exact nature of the solid phase(s). One has to keep in mind that these XRPD patterns are collected under stagnant conditions. Therefore, strong preferential orientation effects could occur (i.e. modifications of peak intensities). Moreover, maturation of the phase can occur which could additionally

Figure 8 .

 8 Figure 8. In situ XRPD patterns of 1. Conditions: In EtOH/2-PrOH (1/1 v/v, 11 fold supersaturation), at 10 °C, and time interval 9 min. These in situ diffraction patterns were compared with that of the dry final powder sample and that simulated from the single crystal data (-form) of racemic 1. Since the crystallization should be induced by the KAPTON ® thin film used as a cover of the observation window, crystallization upon this in situ XRPD measurement was faster than that upon the TCVM observation.

Figure S1 .

 S1 Figure S1. General mechanism of PE in the case of crystallization from a highly supersaturated solution of a slightly S-rich first-generation organic compound. Homochiral 1D molecular association a) in solution and b) in the metastable crystal, and heterochiral 2D molecular association c) in the transformed crystal after polymorphic transition, followed by d) selective redissolution of the excess S enantiomer into the solution and e) deposition of nonracemic mixed crystals enriched with the opposite R enantiomer.

Figure S2 .

 S2 Figure S2. Schematic representation of the intermolecular interactions in the and form crystals. The -form crystals were commonly found as a stable form for the compounds showing PE.

Figure S3 .

 S3 Figure S3. Overview of in situ X-ray diffractometer. The right and left sides are an X-ray source and a semiconductor detector, respectively. The observation cell is set in the cell holder.

Figure S4 .

 S4 Figure S4. In situ XRPD patterns of 1. Conditions: In EtOH/2-PrOH (1/1 v/v, 11 fold supersaturation), at 10 °C. (a) From 90 to 170 min with 9 min interval, and (b) 225 to 1300 min with 240 min interval. These in situ diffraction patterns were compared with that of the dry final powder sample and that simulated from the single crystal data (form) of racemic 1.

(Movie S1-S3). This information is available free of charge via the Internet at http://pubs.acs.org/.

AUTHOR INFORMATION

Corresponding Author takahashi.hiroki.2x@kyoto-u.ac.jp

ACKNOWLEDGMENT

This work was supported by JSPS KAKENHI Grant Number 26248024.

Supporting Information Available: General mechanism of PE (Figure S1); schematic representation of the intermolecular interactions in the-and -form crystals (Figure S2); Overview of in situ X-ray diffractometer (Figure S3); in situ XRPD patterns of 1 (Figure S4); TCVM movie for crystallization of 1 and

Table of Contents Graphic Use Only

In situ Observation of Polymorphic Transition during Crystallization of Organic Compounds Showing Preferential Enrichment by Means of Temperature-Controlled Video-Microscopy and Time-Resolved X-Ray Powder Diffraction

Hiroki Takahashi, * Sekai Iwama, Simon Clevers, Stéphane Veesler, Gérard Coquerel, Hirohito Tsue, Rui Tamura

The stable crystals with a dendritic habit appeared by heterogeneous nucleation onto the metastable leaf-shaped crystals. This polymorphic transition, which was observed by means of temperature-controlled video-microscopy (TCVM) and time-resolved in situ Xray powder diffraction (XRPD), would be responsible for the occurrence of Preferential Enrichment.

In situ Observation of Polymorphic Transition during Crystallization of Organic Compounds Showing Preferential Enrichment by Means of Temperature-Controlled Video-Microscopy and Time-Resolved X-Ray Powder Diffraction

Hiroki Takahashi, *1 Sekai Iwama, 2 Simon Clevers, [START_REF] Jacques | Racemates and Resolution[END_REF] Stéphane Veesler, 4 Gérard Coquerel, [START_REF] Jacques | Racemates and Resolution[END_REF] Hirohito Tsue, 1 Rui Tamura


H F Zhang

G Royer

Theoretical and experimental α decay half-lives of the heaviest odd-Z elements and general predictions

Keywords: numbers: 27.90.+b, 23.60.+e, 21.10.Tg

Theoretical α-decay half-lives of the heaviest odd-Z nuclei are calculated using the experimental Qα value. The barriers in the quasi-molecular shape path is determined within a Generalized Liquid Drop Model (GLDM) and the WKB approximation is used. The results are compared with calculations using the Density-Dependent M3Y (DDM3Y) effective interaction and the Viola-Seaborg-Sobiczewski (VSS) formulae. The calculations provide consistent estimates for the half-lives of the α decay chains of these superheavy elements. The experimental data stand between the GLDM calculations and VSS ones in the most time. Predictions are provided for the α decay half-lives of other superheavy nuclei within the GLDM and VSS approaches using the extrapolated Audi's recent Qα, which may be used for future experimental assignment and identification.

The possibility to synthesize superheavy elements by cold or warm fusion reactions [1,2,3] or using radioactive ion beams has renewed interest in investigating the fusion barriers. The only observed decay mode of these heaviest systems is the α emission, and an accurate description of the α decay is required. The pure Coulomb barrier sharply peaked at the touching point alone does not allow to determine correctly the fusion cross sections and the partial α decay half-lives. In the fusion path, the nucleon-nucleon forces act before the formation of a neck between the two quasispherical colliding ions and a proximity energy term must be added in the usual development of the liquid-drop model [4]. It is highly probable that the α decay takes place also in this fusion-like deformation valley where the one-body shape keeps quasi-spherical ends while the transition between one and two-body configurations corresponds to two spherical nuclei in contact. Consequently, the proximity energy term plays also a main role to correctly describe the α decay barrier. The generalized liquid drop model (GLDM) which includes such a proximity energy term has allowed to describe the fusion [5] , fission [6], light nucleus [7] and α emission [8] processes. The formation and alpha decay of superheavy elements have been investigated [9] in taking into account the experimental Q α value or the value provided by the Thomas-Fermi model [10]. The heaviest even-Z nuclei have been studied [11] using the Q α value obtained experimentally or given by the FRDM [12].

Recently, isotopes of the element 115 have been synthesized [13] and the observed decays reveal that the dominant decay mode is the α emission. These new experimental observations of Z=115 have already attracted a lot of theoretical studies [14,15,16,17,18,19,20,21]. Most of the earlier investigations have been devoted to the description of the ground-state properties of superheavy nuclei, we focus on calculating their half-lives following the first work for even-Z nuclei [11]. In Ref. [21], the α decay half-lives of Z = 115 isotopes are calculated with the microscopic density-dependent M3Y (DDM3Y) interaction, and the results are well consistent with the experimental data. The purpose of this work is to determine the partial α decay half-lives of these superheavy elements within the macroscopic GLDM from the experimental Q α values using the WKB approximation and to compare with the experimental data and the calculations of DDM3Y effective interaction [21] and the Viola-Seaborg formulae with Sobiczewski constants (VSS) [22]. Finally predictions within the GLDM and VSS formulae are given for the partial α decay half-lives of the superheavy nuclei using the recent Q α decay energies of Audi et al. [23].

The GLDM energy is widely explained in [11] and not recalled here. The half-life of the parent nucleus decaying via α emission is calculated using the WKB barrier penetration probability. In such a unified fission model, the decay constant of the α emitter is simply defined as λ = ν 0 P where the assault frequency ν 0 has been taken as ν 0 = 10 19 s -1 , P being the barrier penetrability.

The α decay half-lives of the recently produced odd-Z superheavy nuclei calculated with the three approaches and using the experimental Q α values and without considering the rotational contribution are presented in Table 1. The Q α values given in [23] are obtained by extrapolation. Within the GLDM the quantitative agreement with experimental data is visible. The experimental half-lives are reproduced well in six cases ( 288 115, 284 113, 272 107, 287 115, 283 113, 275 109) out of nine nuclei along the decay chains of 288 115 and 287 115. Two results ( 280 111, 276 109) are underestimated about four to five times possibly because the centrifugal barrier required for the spin-parity conservation could not be taken into account due to non availability of the spin-parities of the decay chain nuclei. On the whole, the results agree well with the experimental data indicating that a GLDM taking account the proximity effects, the mass asymmetry, and an accurate nuclear radius is sufficient to reproduce the α decay potential barriers when the experimental Q α value is known. The results obtained with the DDM3Y interaction agree with the experimental data as the GLDM predictions and largely better than the VSS calculations. This shows that a double folding potential obtained using M3Y [24] effective interaction supplemented by a zero-range potential for the single-nucleon exchange is very appropriate because its microscopic nature includes many nuclear features, in particular a potential energy surface is inherently embedded in this description. This double agreement shows that the experimental data themselves seem to be consistent. For most nuclei the predictions of the VSS model largely overestimate the half lives. The blocking effect is probably treated too roughly. TABLE I: Comparison between experimental α decay half-lives [13] and results obtained with the GLDM, the DDM3Y effective interaction [21] and the VSS formulae. Parent Expt.

[23] Expt. One can also find that all calculated half-lives of the 279 111 nucleus are smaller than the experimental ones in table 1. If the contribution of centrifugal barrier is included, the theoretical results will close the experimental data. On the other hand, it is expected that great deviations of a few superheavy nuclei between the data and model may be eliminated by further improvements on the precision of measurements.

DDM3Y GLDM GLDM VSS VSS Nuclei Q[MeV] Q[MeV] T 1/2 T 1/2 (Qex) T 1/2 (Qex) T 1/2 (Q Audi ) T 1/2 (Qex) T 1/2 (Q Audi )
Another noticeable point is that the experimental α decay half-lives are between the close theoretical values given by the GLDM and the ones derived from the VSS formulae. Thus predictions of the α decay half lives with the GLDM and VSS formulae are possible as long as we know the right α decay energies. The ones derived from Audi's recent publication [23] are very close to the experimental data. The most deviation is not more than 0.5 MeV, which is a valuable result for studying correctly the half-lives. The calculations using the α decay energies of Ref. [23] for the nuclei of the 288 115 and 287 115 decay chains by the GLDM and VSS formulae are reasonably consistent with the experimental data. The experimental data stand between the calculations of the GLDM and the results of VSS in six cases for the seven nuclei when experimental uncertainty in the Q value is considered. Thus, predictions of the halflives of superheavy nuclei with the GLDM and VSS formulae are provided for a large number of superheavy elements in Table 2 using the extrapolated Q α values given by [23] or the experimental data indicated by an asterisk. They are an improvement relatively to the values previously given in [9] since these extrapolated Q α values are in better agreement with the experimental data than the ones proposed in [10]. It may be useful for the future experimental assignment and identification.

In conclusion, the half-lives for α-radioactivity have been analyzed in the quasimolecular shape path within a Generalized Liquid Drop Model including the proximity effects between nucleons and the mass and charge asymmetry. The results are in agreement with the experimental data for the alpha decay half-lives along the decay chains of the Z=115 isotopes and close to the ones derived from the DDM3Y effective interaction. The experimental α decay half lives stand between the GLDM calculations and VSS formulae results and the α decay half-lives of some superheavy nuclei have been presented within the GLDM and VSS approaches and Q α adopted from the Audi's recent extrapolated data.

H.F. Zhang is thankful to Prof. U. Lombardo, Junqing Li, Wei Zuo, Baoqiu Chen and Zhong-yu Ma for valuable discussions. H.F. Zhang also would like to express his gratitude for the kind hospitality and excellent working conditions at INFN-LNS. Partial support from the Natural Science Foundation of China (NSFC) under Grant No. 10505016, TABLE II: Predicted α-decay half-lives using the GLDM and the VSS formulae, the α decay energies are taken from the extrapolated data of Audi et al. [23] or the experimental data indicated by an asterisk.

Nuclei Q [ MeV ] T GLDM 1/2 T V SS 1/2 Nuclei Q [ MeV ] T GLDM 1/2 T V SS 1/2 Nuclei Q [ MeV ] T GLDM 1/2 T V