Menghua Zhao

Julien Dervaux

Tetsuharu Narita

Francois Lequeux

Laurent Limat

Matthieu Roché
email: matthieu.roche@univ-paris-diderot.fr

Geometrical control of dissipation during the spreading of liquids on soft solids

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Gel layers bound to a rigid substrate are used in cell culture to control differentiation and migration and to lower the friction and tailor the wetting of solids. Their thickness, often considered a negligible parameter, affects cell mechanosensing or the shape of sessile droplets. Here we show that the adjustment of coating thickness provides control over energy dissipation during the spreading of flowing matter on a gel layer. We combine experiments and theory to provide an analytical description of both the statics and the dynamics of the contact line between the gel, the liquid and the surrounding atmosphere. We extract from this analysis a hitherto unknown scaling law that predicts the dynamic contact angle between the three phases as a function of the properties of the coating and the velocity of the contact line. Finally, we show that droplets moving on vertical substrates coated with gel layers having linear thickness gradients drift towards regions of higher energy dissipation. Thus thickness control opens the opportunity to design a priori the path followed by large droplets moving on gel-coated substrates. Our study shows that thickness is another parameter, besides surface energy and substrate mechanics, to tune the dynamics of liquid spreading and wetting on a compliant coating, with potential applications in dew collection and free-surface flow control.

Wetting | Soft Materials | Interfacial Science | G els are soft complex materials made of a polymeric scaf- fold in which a liquid may be embedded [START_REF] Larson | The Structure and Rheology of Complex Fluids[END_REF]. They can be obtained from a broad catalogue of macromolecules and solvents such as poly(dimethylsiloxane) (PDMS) or water-based solutions of poly(vinyl alcohol) (PVA). In turn, gels offer great flexibility with respect to the tuning of their mechanical and physicochemical properties. These materials find uses in a diverse range of applications, such as cell culture (2), regenerative medicine, drug delivery [START_REF] Jagur-Grodzinski | Polymeric gels and hydrogels for biomedical and pharmaceutical applications[END_REF], or the modification of solid surfaces to control friction, wetting and heterogeneous nucleation [START_REF] Gong | Friction and lubrication of hydrogels -its richness and complexity[END_REF][START_REF] Sokuler | The softer the better: Fast condensation on soft surfaces[END_REF][START_REF] Style | Elastocapillarity: Surface Tension and the Mechanics of Soft Solids[END_REF]. Finally, gels as coatings have attracted renewed interest in recent years due to the opportunity they offer to generalize existing descriptions of the wetting of solids to arbitrary materials [START_REF] Lester | Contact angles of liquids at deformable solid surfaces[END_REF][START_REF] Shanahan | L'arête produite par un coin liquide près de la ligne triple de contact solide/liquide/fluide[END_REF][START_REF] Shanahan | The influence of solid micro-deformation on contact angle equilibrium[END_REF][START_REF] Long | Static and dynamic wetting properties of thin rubber films[END_REF][START_REF] Yu | Role of vertical component of surface tension of the droplet on the elastic deformation of PDMS membrane[END_REF][START_REF] Yu | Elastic deformation of soft membrane with finite thickness induced by a sessile liquid droplet[END_REF][START_REF] Pericet-Camara | Solid-supported thin elastomer films deformed by microdrops[END_REF][START_REF] Jerison | Deformation of an elastic substrate by a threephase contact line[END_REF][START_REF] Limat | Straight contact lines on a soft, incompressible solid[END_REF][START_REF] Style | Static wetting on deformable substrates, from liquids to soft solids[END_REF][START_REF] Style | Universal deformation of soft substrates near a contact line and the direct measurement of solid surface stresses[END_REF][START_REF] Style | Surface tension and contact with soft elastic solids[END_REF][START_REF] Bostwick | Elastocapillary deformations on partially-wetting substrates: rival contact-line models[END_REF].

Most of the reports dealing with wetting involve a sessile droplet sitting on a solid with a large shear modulus µ0 O(GPa) [START_REF] De Gennes | Wetting : statics and dynamics[END_REF][START_REF] Bonn | Wetting and spreading[END_REF]. In this context, the normal component F ⊥ cap of the resulting capillary force at the contact line between the droplet, the rigid substrate and the atmosphere deforms the interfacial region of the solid on molecular scales. As a result, this contribution is neglected in the rationalization of wetting on solids, and the problem is usually described with the Young-Dupré equation. In contrast, gels can be soft enough (1 ≤ µ0 ≤ 10 4 Pa) that F ⊥ cap deforms their surface on micrometer scales, much larger than the size of the polymers and of the solvent molecules. Due to progresses in experimental methods, many recent studies have investigated how the presence of this deformation, known as a ridge, may alter the force balance at the contact line and modify the shape of both the sessile drop and the substrate [START_REF] Pericet-Camara | Solid-supported thin elastomer films deformed by microdrops[END_REF][START_REF] Jerison | Deformation of an elastic substrate by a threephase contact line[END_REF][START_REF] Style | Static wetting on deformable substrates, from liquids to soft solids[END_REF][START_REF] Style | Universal deformation of soft substrates near a contact line and the direct measurement of solid surface stresses[END_REF][START_REF] Style | Surface tension and contact with soft elastic solids[END_REF][START_REF] Marchand | Capillary Pressure and Contact Line Force on a Soft Solid[END_REF].

The presence of a soft coating also modifies wetting dynamics. The velocity V at which a contact line between a solid, a liquid and a gas moves on a rigid substrate coated with a soft elastomeric layer is smaller than the velocity of the contact line on the bare substrate. This velocity reduction, known as viscoelastic braking [START_REF] Carré | Viscoelastic effects in the spreading of liquids[END_REF], results from an additional mechanism of energy dissipation due to the viscoelasticity of the gel. As a consequence, the spreading of liquids and the motion

Significance Statement

Soft solids such as gels are routinely used to coat solid surfaces to control wetting and adhesion. While wetting statics on these coatings has received a lot of attention, the dynamics has been hardly studied. Here we show that energy dissipation during spreading depends not only on the viscoelasticity of the gel but also on the thickness of the coating. We provide both an analytical model and a scaling law that predict the dynamic contact angle, i.e. the signature of dissipation in the dropletsubstrate-atmosphere system, as a function of the mechanical properties of the gel and the thickness of the coating. Finally, we show that droplets moving on coatings having thickness gradients drift towards regions of higher dissipation.

D R A F T

of droplets are slowed down. Nonetheless, this description of the problem is incomplete: Figures 1b-c demonstrate that the velocity of a water droplet moving down a vertical glass slide coated with a soft viscoelastic silicone layer (µ0 = 1.1 kPa, used for all experiments in this study, see Materials for a rheological characterization) also decreases as the elastomer thickness increases. The purpose of this paper is to uncover the link between substrate thickness and the motion of a contact line.

Results

As we want to understand a dynamic wetting observation, we must characterize the reference state, i.e. the statics of wetting. We investigated the properties of the region surrounding the contact line between a sessile droplet and a silicone elastomer layer using both side-view microscopy to measure the apparent contact angle at equilibrium θeq and a non-invasive optical method based on Schlieren photography (24) (Fig. 2a). We find that θeq depends on the contact radius R of the droplet on the elastomer at constant sample thickness h0 but it is independent of h0 for constant R (Fig. 2b-c), accounting for measurement uncertainties. For large droplets, θeq tends towards a saturation at a value θ∞ 106 • . The Schlieren results show that the connection between the ridge and the flat elastomer surface far from the droplet is monotonically decreasing for large thicknesses while the surface profile shows a sub-micrometre surface depression in the vicinity of the ridge at small thicknesses, known as a micro-trough [START_REF] Pericet-Camara | Solid-supported thin elastomer films deformed by microdrops[END_REF][START_REF] Jerison | Deformation of an elastic substrate by a threephase contact line[END_REF]. The micro-trough results from the incompressibility of the gel, which has to accommodate the Poisson effect at small thicknesses under the zero-motion constraint at its basis due to its bonding to a rigid substrate. Its depth ζ increases with a decrease of h0 and an increase of the droplet radius R (see Supplementary Materials).

We interpret these observations with an analytical theoretical description of the statics of elasto-wetting in the limit of linear elasticity developed recently by Dervaux and Limat (25) that we extend to account for finite-depth effects (see Supplementary Materials). The model (and the model used in the dynamic context later) considers a rivulet, i.e. a column of liquid whose axis of symmetry is parallel to the surface of the gel, sitting on the surface of the coating. As the angles that we have measured are rather close to π 2 , we assume that the interfacial tensions of the solid with the liquid and the gas are identical, γSL = γSV = γs, and this quantity is the only fitting parameter. The model is able to capture the surface profiles of the elastomer for all the thicknesses that we investigated (Fig. 2d) as well as the evolution of ζ as a function of h0 and R (see Supplementary Materials). The model indicates that γs = 40 mN m -1 , a value comparable to those reported earlier [START_REF] Style | Universal deformation of soft substrates near a contact line and the direct measurement of solid surface stresses[END_REF][START_REF] Karpitschka | Droplets move over viscoelastic substrates by surfing a ridge[END_REF]. To sum up, the surface profile of the elastomer layer in the vicinity of the contact line depends on layer thickness. In contrast, the apparent equilibrium contact angle θeq is independent of the geometry of the gel layer; it depends only on the surface energies of the materials , a feature that is predicted by the model (See Supplementary Materials). The analysis of the statics of wetting allows us to define the reference state of the study of wetting dynamics as a function of the thickness of the elastomer layer to which we now turn.

Wetting dynamics is characterized by the measurement of the dependence of the difference ∆θ = θeq -θ dyn , the latter being the dynamic contact angle, on the velocity of the contact line V [START_REF] Bonn | Wetting and spreading[END_REF]. As θeq > π 2 , we studied a receding contact line surrounding a droplet with an initial contact radius R >> 1 mm that we deflated at a constant flow rate. We recall that θeq reaches its saturation value θ∞ for such large droplets. For all thicknesses, ∆θ increases with an increase of V (Fig. 3a). The dynamic contact angle differs by almost 30 • from the static equilibrium contact angle for the largest velocities we investigated. The two curves show that ∆θ is also a function of thickness: it grows by 30% with an augmentation of h0 from 10 to around 100 µm (Fig 3b). ∆θ saturates to its velocity-dependent value when h0 > 100 µm. To interpret these observations, we recall that, in general, large energy dissipation induces large values for ∆θ at fixed contact line velocity [START_REF] Bonn | Wetting and spreading[END_REF]. In this context, our results show that the thinner the elastomer is, the more energy dissipation is reduced. Voué et al. reported a similar observation for a silicone oil wetting a silicone elastomer [START_REF] Voué | Dissipation and moving contact lines on non-rigid substrates[END_REF]. As the liquid can permeate through the gel in their system, the wetting dynamics they report are difficult to interpret, especially since the thickness of their elastomer layers is of the order of 10 µm. Permeation does not occur in our system. We now unveil how thickness and dissipation are related in elastowetting.

To identify the mechanism relating dissipation to the thickness of the coating, we developed an analytical model based on linear viscoelasticity that is a generalization to large velocities and arbitrary rheology of a framework proposed by Long et al. [START_REF] Long | Static and dynamic wetting properties of thin rubber films[END_REF]. In this model, dissipation in the solid balances the energy that sets the contact line in motion. We assume that the shear relaxation modulus of the soft layer is described by the Chasset-Thirion model, [START_REF] Long | Static and dynamic wetting properties of thin rubber films[END_REF]. We find this assumption to be valid for our material with µ0 = 1085 ± 124 Pa, m = 0.66 ± 0.04 and τ = 15.4 ± 0.4 ms (see Materials). We also use the surface tension of the solid γs = 40 mN m -1 obtained from the statics. The model predicts that the function describing the dependence of the dynamic contact angle on the properties of the substrate and the velocity of the contact line G(θ dyn ) = cos (θeq )-cos (θ dyn ) sin2 (θ dyn )

µ(ω) = G (ω) + iG (ω) = µ0(1 + (iωτ ) m )
is proportional to V m , in excellent agreement with the data (Fig. 3a), with an adjusted exponent m adj = 0.62 ± 0.02 close to the experimental value of m = 0.66 ± 0.04 we obtained from rheology (inset in Fig. 3a). The full model also captures the dependence of G(θ dyn ) on h0 at constant velocity of the contact line (Fig. 3b). Finally, we note that we did not observe effects related to the Laplace pressure inside the droplet. Thus the 2D results obtained with the model are valid in 3D (see Supplementary Materials). We can extract from the model scaling laws that describe the asymptotic behavior of G(θ dyn ) for thick and thin coatings. We assume that the problem is invariant with a translation along the contact line. Our considerations are based on the Fourier transform of the problem. The power per unit length of the contact line injected by the driving capillary force is:

Pin ∝ γV (cos (θ dyn ) -cos (θeq)) [1]
where γ is the interfacial tension of the liquid-gas interface. Power Pin is balanced by dissipation in the gel. Using the Chasset-Thirion model, we estimate the power dissipated in the gel per unit volume: where ω is a typical frequency that we assume to be of the order of V , with a length scale to be determined. The dissipated power per unit length of the contact line is obtained by integrating Eq. 2 over 2 , as this length scale also corresponds to the typical vertical displacement experienced by the solid as well as the horizontal extent of the moving ridge. Taking = γ γs sin (θ dyn ) as a typical scale for strains:

p dis ∝ µ0(ωτ ) m 2 ω [2] D R A F T
P dis ∝ µ0( V τ ) m V γ γs sin (θ dyn ) 2 [3]
Equating Eqs. 1 and 3, we find that:

G(θ dyn ) ∝ γ γ 2 s µ0 V τ m [4]
The length is related to the balance between capillarity and elasticity in the system. Its expression depends on whether the coating has a finite thickness or not. We can express in the limit of thin and thick samples by investigating the asymptotics of the amplitude A of the Fourier transform of the Fourier transform of the surface vertical displacement field ζ(k) (see Eq. 20 in the Supplementary Materials):

A(k) = 1 γsk 2 + µ0K(k) -1 [5]
where k is the wavenumber and:

K(k) = (2k) -1 sinh (2h0k) -2h0k 2h 2 0 k 2 + cosh (2h0k) + 1 . [6]
For large thicknesses, h0 → ∞, and we obtain that K(k) → (2k) -1 . Then:

A(k) → 1 γsk 2 1 1 + 2 µ 0 γsk . [7]
The characteristic length scale in Eq. 7 is = γs µ 0 , i.e. the elastocapillary length. Injecting this value of in Eq. 4, we find the scaling predicted 20 years ago by Long et al. [START_REF] Long | Static and dynamic wetting properties of thin rubber films[END_REF] :

G(θ dyn ) ∝ γ γs µ0V τ γs m . [8]
In the limit of thin coatings, we carry out the Taylor expansion of K(k) around h0 = 0 and we keep the first non-vanishing term:

K(k) ∼ h 3 0 k 2 [9]
We find that the characteristic length scale is =

γsh 3 0 µ 0 1/4
in this case. Injecting in Eq. 3, we obtain the following expression for the dissipated power:

P diss ∝ µ0 2 V (V τ ) m γsh 3 0 µ0 (1-m)/4
[10]

and we obtain a new scaling law reflecting thickness effects on the dynamic contact angle for thin coatings: Equation 11 can be tested by rescaling the experimental data as in figure 3c. We see that the data follow Eq. 11 when h0 << s = γs µ 0 . When h0 >> s, the data is independent of h0 and it follows Eq. 8. We are the first to report the validity of Long et al.'s model [START_REF] Long | Static and dynamic wetting properties of thin rubber films[END_REF] to the best of our knowledge. Thus, our results show that dissipation during the wetting of a soft coating by a liquid can be engineered a predictable way by adjusting the thickness of the layer. Finally, the range of thicknesses over which we observe a dependence of dissipation on h0 corresponds to the typical range of coating thicknesses in applications.

G(θ dyn ) ∝ γLV γs V τ µ0 γs

Discussion

The prediction that the dynamic contact angle at a contact line moving on a soft coating is closer to the value of the static contact angle as the thickness of the coating decreases may seem counter-intuitive at first. To gain more insight, we shall now compare contact line motion for a rivulet on a gel layer and on a liquid substrate of similar thickness h l . In the latter case, continuity of stresses and velocities is enforced at all the interfaces. Momentum diffuses from the droplet to the bulk of the liquid substrate through the droplet-substrate interface and stresses and dissipation do not diverge at the contact line. Energy dissipation occurs in the bulk of both the droplet and the substrate. In the latter, the dissipation power per unit length along the contact line between the rivulet and the gel scales as Pvis ∝ ηV 2 R 0 (t) h l [START_REF] Landau | Fluid Mechanics[END_REF], where R0(t) is the radius of the spreading rivulet; we assume that h l remains constant over time. Thus, dissipation increases with a decrease of h l [START_REF] Craster | On the dynamics of liquid lenses[END_REF][START_REF] Karapetsas | Surfactant-driven dynamics of liquid lenses[END_REF]. In contrast, at the interface between a gel coating and a liquid, we enforce the no-slip boundary condition. Thus velocities are zero everywhere along this interface, and the contact line is the locus of most dissipation. In this context, we have shown that dissipation follows Eq. 10. As m < 1, dissipation decreases with a decrease of the thickness h0 of the sample.

We also computed the prediction for the dependence of G(θ dyn ) using another model proposed recently by Karpitschka et al.(26). This model captures well the data for ∆θ(V ) at large thickness provided we adjust the surface tension of the solid to a value different from the one we obtain in our study of the statics of wetting. Nonetheless, it does not capture the decreasing trend of ∆θ at small thicknesses (see Supplementary Materials). In their model, Karpitschka et al. impose that capillary forces are at equilibrium at the moving contact line at all times. The failure of their model to capture our data suggests that this constraint is not true anymore for droplets moving on thin coatings. The question of the validity of this assumption at all thicknesses remains open as its answer lies in a second-order analysis of the dynamics of elastowetting

D R A F T

that has yet to be performed. We take advantage of our findings to control the trajectory of liquids moving at the surface of solid substrates. We coat a glass slide with a soft elastomeric layer with a thickness gradient (Fig. 4a). We place a water droplet at one extremity of the glass slide above the thinnest part of the coating. Then we incline the glass slide to the vertical, the thickness gradient being perpendicular to gravity. The droplet slides down and it drifts from the thinnest regions towards the thickest ones. Drift results from the dependence of dissipation in the gel on h0: the droplet experiences a larger braking over the thick part of the coating than other the thin region (Fig. 4b). Thus, the droplet is subject to a torque that leads to a global drift towards the thick region. In our system, drift will be observed if the gradient starts from h0 < 100 µm. If the thin part of the thickness gradient is thicker than 100 µm, we expect no drift, as is observed in the experiments. Finally, we observe no spontaneous motion of droplets when the glass slide remains horizontal for a long time. This fact highlights the difference between our results, whose origin lies in the variation of the dynamic contact angle around the droplet, and those reported by Style et al. [START_REF] Style | Patterning droplets with durotaxis[END_REF] concerning droplet durotaxis, which occurs due to changes of the static contact angle along the droplet perimeter.

Our findings opens the possibility to design soft coatings whose mechanical dissipation can be adjusted not only by formulation but also by an appropriate choice of thickness. We envision improvements in applications such as dew collection, where the use of soft layers has been shown to improve the nucleation rate of condensation droplets(5); the use of thin layers will help increase the up-time of the condensation cell. We also anticipate the use of our results to control free-surface flows such as droplet impact or gravity currents passively with thickness-gradient coatings.

Materials and Methods

Gel preparation and characterization. All the experiments were performed with Millipore-Q purified water (γ = 72 mN/m, ) on polydimethylsiloxane (PDMS) substrates made from the commercial elastomer Sylgard 527 (Dow Corning, USA). Sylgard 527 is a twocomponent elastomer kit composed of a base and a cross-linker. The two liquids were mixed at a manufacturer-recommended 1:1 ratio, to respect crosslinking stoichiometry, with a magnetic stirrer for approximately 30 minutes. The mixture is then let to degas under vacuum (P 50 mbar) for at least one hour. We obtained thin flat PDMS gel films (thickness h 0 < 100µm) by spin-coating directly the degassed mixture on glass slides. Thickness was controlled by tuning the spinning speed from to 3400 RPM. For thick flat films (thickness > 100µm), the thickness was tuned by changing the elastomer volume poured onto the glass slide mounted with movable walls.

For the thickness gradient experiment, the Sylgard mixture was poured into a mold containing a titled glass slide. Next, the elastomer was cross-linked at 65 • C in an oven for 24 hours. This protocol produced transparent PDMS films strongly bound to glass slides. The thickness of the films was measured with a 3D Profiler (FOCAL 3D Pilot, FOCAL Nanotech, France) operating in white light scanning mode with a precision of the order of a nanometer.

Figure 5 shows the rheological data obtained for Sylgard 527 and measured with a strain-controlled rheometer (Physica MCR 500, Anton Paar, Austria). The elastomer was cross-linked directly between the Peltier plate and the plate tool (PP20-MRD, d = 20 mmAnton Paar, Austria) at T = 65 • C for 4 hours. The sample was then let to cool down to 25 ± 0.2 • C before running the test. The gap was set at 0.6 mm and the strain was fixed at 1%. A fit of the Chasset-Thirion model, µ(ω) = G (ω) + iG (ω) = µ 0 (1 + (iωτ ) m ), to the experimental data provides the parameters for the frequencydependent mechanical response of the gel. We also characterized

D R A F T

In principle, the slit light source is perfectly refocused on the knife blade plane as illustrated by the green light beam drawn in Fig. 2 (a), as long as there is no disturbance on the light path. However, changes in the surface slope of the gel film deflect the light beam and the conjugated image is thus shifted by a distance ∆a in the focal plane as illustrated by the red beam shown in Fig. 2 (a). Under the small angle approximation, this distance is related to the deflection of the light beam through the Schlieren object by the formula ∆a = f where f is the focal length of L3. For an accurate quantification of the surface profile, we carried out a three-step measurement. First, a sequence of pictures was recorded at different vertical blade positions and a map relating the pixel grayscale value to the deflection by the knife cut-off is built. During this step, the Schlieren object is absent. Second, all the camera parameters are kept identical and pictures are taken 5 mins after the drop fully relaxes on the soft substrate (for the static deformation measurement). Then, the grayscale value of the Schlieren picture is compared to the calibration map and we obtain a deflection map. Finally, by combining the deflection relation and Snell's law, the map of the surface slope is reconstructed. Our system has a deflection angle sensitivity of 0.0007 rad, and we could measure angles up to 0.05 rad.

Data acquisition. To study the statics of wetting, we enclosed the sample in a transparent glass chamber (75 × 50 × 8mm) to control humidity and avoid evaporation. Its lid could easily be removed to place drops. For experiments, a reservoir of water was enclosed together with the PDMS slides to saturate the atmosphere of the glass chamber. All the measurements were performed 5 minutes after deposition. For the spreading dynamic experiments, a side view camera (DMK 23UV024 Imaging Source, Germany, C1 in Fig. 2 a) equipped with a lens (TV Lens 50mm, f1:2.8, Ricoh, Japan) and two tube sets (C Mount TV Lens Extension Tube Set 40 mm, JAPAN), was used. The spatial resolution of C1 was 5.5 µm px -1 . Using a Gaussian subpixel detection technique, the error in the contact line velocity determination is around 10 µm s -1 . All the drops were deposited using a syringe which was connected to a pump (PUMP 33, HARVARD APPARATUS, USA) at room temperature. For the sliding and drifting experiments on tilted slides, droplets were deposited using a micropipette (0.5-10 µL, FinnPipette) at ambient temperature. A front view camera (DFK 23UX174, Imaging Source, Germany) equipped with a lens (TV Lens 50mm f1:2.8, Ricoh, Japan) and one tube set (C Mount TV Lens Extension Tube Set 20 mm, JAPAN), with a spatial resolution of 15.9 µm px -1 was used to track the sliding droplets.

Model. We use two models, one for the statics and one for the dynamics. The model describing the statics is built in the framework of linear elasticity and we account for the surface tension of the solid. We assume that the system is translation-invariant and we describe a x-y slice of the substrate subjected to capillary force of a liquid. We model the latter as a concentrated load on an elastic half-space of depth h 0 . We are interested in two situations, the single contact line and the rivulet, i.e. a column of liquid that has two straight contact lines with the substrates. We focus on the far-field vertical displacement for which possible divergences in the displacements at the contact line do not matter [START_REF] Dervaux | Contact lines on soft solids with uniform surface tension : analytical solutions and double transition for increasing deformability[END_REF]. We use Fourier transforms to solve the biharmonic equation that describes the problem, and we compare the data with the outcomes of the description of the rivulet. The static model also accounts for the small hysteresis that we measure in experiments in the form of a force tangential to the surface.

The model for the dynamics is built in the framework of linear viscoelasticity. We use the surface tension of the solid deduced from the study of the statics as an input to the model. As experiments indicate that the properties of the motion of the contact line are independent of the size of the droplet, we focus on the case of a single contact line. We compute the dependence of the vertical displacement on the velocity of the contact line using a similar procedure to the statics resorting to Fourier transforms. We then compare the power injected in the system due to capillarity with the dissipated power in the substrate due to its viscoelasticity.

1 a

 1 Geometrical control of dissipation during the spreading of liquids on soft solidsMenghuaZhao a,b , Julien Dervaux a , Tetsuharu Narita b,c , François Lequeux b , Laurent Limat a , and Matthieu Roché a,Matière et Systèmes Complexes, CNRS UMR 7057, Université Paris Diderot, SPC University, 10 Rue A. Domon et L. Duquet, F-75013 Paris, France; b Laboratoire Sciences et Ingénierie de la Matière Molle, PSL Research University, ESPCI Paris, CNRS, 10 rue Vauquelin, F-75231 Paris cedex 05, France; c Global Station for Soft Matter, Global Institution for Collaborative Research and Education, Hokkaido University, Sapporo, Japan This manuscript was compiled on July 12, 2018

Fig. 1 .

 1 Fig. 1. Sliding of a droplet on a soft layer. (a) A water droplet is deposited on the surface of a glass slide covered with a poly(dimethylsiloxane) layer of thickness h0. The glass slide is inclined to the vertical at t = 0. (b,c) Motion of 2µL droplets on glass slides covered with PDMS layers 9 µm and 270 µm thick. Labels: time in minutes elapsed since inclination. Scale bar: 1mm.

Fig. 2 .

 2 Fig. 2.Experimental set-up and characterization of the static three-phase contact line. (a) The region surrounding the contact line between water, the elastomer and the atmosphere is visualised simultaneously from the side with camera C1 and from the bottom with camera C2 that captures the output of a Schlieren photography set-up. The former allows us to measure the contact angle θ and the contact radius between the droplet and the solid R. The latter is used as a surface profilometer along the yellow dashed line. Scale bar on the images: 0.5 mm. (b,c) Dependence of the equilibrium value of the contact angle θeq on R for four different layer thicknesses and on the thickness h0 of the soft layer for sub and supra-millimetric droplets. In the latter graph, the dashed lines represent a spread of one standard deviation around the average value of θeq over all thicknesses. (d) Dependence of the surface position ζ on the distance to the contact line for two different thicknesses. The origin of the x axis is arbitrary; we use it to note that part of the contact line is hidden to the Schlieren set-up by the droplet as θeq > π/2. Dashed lines: predictions of our model.

m h0µ0 γs 3 ( 1 -m)/ 4 [ 11 ] 3 DV 1 Fig. 3 .

 31411313 Fig. 3. Characterization of the dynamics of the receding three-phase contact line on soft elastic layers. (a) Dependence of the contact angle difference θeq -θ dyn on the velocity of the receding contact line on two layers of thickness h0 = 15 and 900 µm. Dashed lines: predictions based on our model. Inset : log-log plot of the function G(θ dyn ) as a function of the contact-line velocity for the same systems. Dashed line: power-law fit with an exponent m = 0.62 ± 0.02. (b) Dependence of θeq -θ dyn on h0 at constant receding velocity V = 0.1 mm s -1 . Dashed lines: predictions of our model for two values of the solid surface tension γs. (c) Test of Eq. 11, relating the dynamic contact angle to the velocity of the contact line and the thickness of the sample in the small-thickness regime against experimental data. Dashed line: Eq. 11, with m = 0.62 ± 0.02 and a slope of 2. The horizontal line is a guide for the eyes. (d) Schematic of the motion of a point of the surface of the elastomer (red dot) during the passage of the contact line.Time starts at the top left and ends at the bottom right.

Fig. 4 .

 4 Fig. 4. Drifting of a droplet moving down a vertical glass slide covered with a soft layer of increasing thickness. (a) a 2-µL droplet moving over a coating with a thickness gradient ∇h ∼ 70 µm mm -1 . Scale bar: 1 mm. (b) Schematic explaining the mechanism leading to droplet drift during forced motion. Black dashed line is the trajectory of the droplet, white dashed line is the vertical. The orange arrow represents the rotation of the droplet as it moves down the coating and experiences differences in dissipation along its contour, in a fashion analogous to a car experiencing braking only on one of its sides.

5 .

 5 Small-amplitude oscillatory shear measurement of the mechanical properties of the Sylgard 527 elastomer with a strain-controlled rheometer (Physica MCR 500, Anton Paar, Austria). Dashed line: best fit of the Chasset-Thirion model, µ(ω) = G (ω) + iG (ω) = µ0(1 + (iωτ ) m ), to the experimental data. µ0 = 1085 ± 124 Pa, m = 0.66 ± 0.04 and τ = 15.4 ± 0.4 ms the contact angle hysteresis of water on Sylgard 527 and we found it to be of the order of 3-4 • .Quantitative Schlieren optics. We built an observation platform for the measurements of out-of-plane deformations of the soft film based on Schlieren optics. A white LED light source (Luxeon, Lumileds, USA) is focused on a mechanical slit (VA100/M, Thor Labs, USA, size 1 × 13.6 mm) by a condenser L1. The light beam from the slit is then collimated by a lens L2 (diameter d L2 = 25.4 mm; focal length f = 101.6 mm; MPD 149-P01, Thor Labs, USA). The parallel light beam goes through the sample between mirror prisms M1 and M2. It is then collected by a third lens L3 (LB1374-B, Thor Labs, USA). To obtain a Schlieren set-up, a knife blade is placed at the focal point of lens L3 and hides part of the focal region. The blade filters those light beams that have been deviated because of surface deformations of the gel sample due to the presence of the droplet. Images of the filtered focal point are captured with a digital camera (DFK 23UX174, Imaging Source, Germany) equipped with a lens (Avenir TV Zoom Lens F1.8, JAPAN). The CCD sensor of the camera has a dimension of 1280 × 960 pixel, resulting in a lateral spatial resolution of 10.2 µm px -1 .

| www.pnas.org/cgi/doi/10.1073/pnas.XXXXXXXXXX Zhao et al.

| www.pnas.org/cgi/doi/10.1073/pnas.XXXXXXXXXX Zhao et al.

| www.pnas.org/cgi/doi/10.1073/pnas.XXXXXXXXXX Zhao et al.

ACKNOWLEDGMENTS. The authors thank

Bruno Andreotti and Jacco Snoeijer for interesting discussions. Menghua Zhao thanks the China Scholarship Council for the PhD fellowship funding. Agence Nationale de la Recherche and Commissariat à l'Investissement d'Avenir are gratefully acknowledged for their financial support through the MMEMI grant (LabEx SEAM (Science and Engineering for Advanced Materials and devices), ANR-11-LABX-086, ANR-11-IDEX-05-02


Bryan Penin
email: bryan.penin@inria.fr

Paolo Robuffo

Franc ¸ois Chaumette
email: francois.chaumette@inria.fr

F Chaumette

P Robuffo Giordano

Vision-Based Reactive Planning for Aggressive Target Tracking while Avoiding Collisions and Occlusions

Keywords: Reactive and Sensor-Based Planning, Optimization and Optimal Control, Visual-Based Navigation

In this paper we investigate the online generation of optimal trajectories for target tracking with a quadrotor while satisfying a set of image-based and actuation constraints. We consider a quadrotor equipped with a camera (either down or front-looking) with limited field of view. The aim is to follow in a smooth but reactive way a moving target while avoiding obstacles in the environment and occlusions in the image space. We propose vision-based approaches based on multi-objective optimization, especially with the occlusion constraint formulation. We design an online replanning strategy inspired from Model Predictive Control (MPC) that successively solves a non-linear optimization problem. The problem is formulated as a Nonlinear Program (NLP) using differential flatness and finite parametrization with B-Splines. This allows a resolution by Sequential Quadratic Programming (SQP) at a rate of 30Hz. The robustness and reactivity of the replanning algorithm are demonstrated through realistic simulation results. Experiments validating the performance with a real quadrotor are also presented.

I. INTRODUCTION

U NMANNED Aerial Vehicles (UAVs) are pushing away the boundaries for applications covering agile navigation, inspection and surveillance to name a few. To be reactive to changes in the environment, the use of on-board cameras has become a fundamental necessity. Indeed, moderately invasive they come with numerous computer vision algorithms for tracking objects, detecting obstacles but also for giving reliable measurements for non-linear observers (e.g., active vision, state estimation). For these reasons, incorporating visual cues in the loop has become a popular approach in robotic control for designing robust and reactive control laws. A historical technique is visual servoing [START_REF] Chaumette | Visual servo control. I. Basic approaches[END_REF]. However, this method cannot be directly applied to quadrotors due to their complex dynamics and inherent under-actuation that generates rotations of the camera which may conflict with the main servoing task. Although, dynamic visual servoing schemes have been developed for second order or under-actuated systems (e.g., Fig. 1. An optimal trajectory for a quadrotor equipped with a camera is generated for reaching a minimum distance with a target (blue sphere) while avoiding collisions with spherical obstacles (inflated dark spheres) and occlusions of the target from the obstacles (red spheres) [START_REF] Bourquardez | Image-Based Visual Servo Control of the Translation Kinematics of a Quadrotor Aerial Vehicle[END_REF], [START_REF] Jabbari | An Adaptative Scheme for IBVS of an Underactuated UAV[END_REF] or [START_REF] Mebarki | Nonlinear visual control of unmanned aerial vehicles in gps-denied environments[END_REF] for quadrotors), the underlying assumptions fail for high speed manoeuvres and in any case, do not take into account possible loss of visibility or occlusions.

Preserving visibility is a substantial concern when vision is in the loop since losing track of features may lead to a failure of the task. Visual servoing original formulation does not prevent critical configurations related to: 1) field of view limits, i.e., the target may become invisible as the camera moves, especially with underactuated robots such as quadrotors; 2) occlusions, i.e., the image features may be occluded due to obstacles or body parts. 3) physical constraints such as feasibility constraints 4) collision with obstacles or self-collision.

When it comes to vision-based control, potential fields are classically used for designing control laws that repel visual cues from the obstacles in the image [START_REF] Cowan | Visual servoing via navigation functions[END_REF], [START_REF] Mezouar | Path planning for robust image-based control[END_REF] but this may strongly conflict with the nominal servoing task and increase the chance of falling into local minima. Gradient Projection Methods (GPM) use the system redundancy to mitigate the completion of two tasks [START_REF] Marchand | Dynamic sensor planning in visual servoing[END_REF], [START_REF] Mansard | A new redundancy formalism for avoidance in visual servoing[END_REF].

However the redundancy formalism does not appear reasonable when dealing with under-actuated robots such as quadrotors.

In contrast with the aforementioned approaches the present work explicitly considers the full underactuation of the camera which is often problematic for standard visual servoing laws.

Optimization-based planning has motivated a great deal of research with its ability to handle multiple types of constraints and to regulate multiple and possibly conflicting objectives. Nowadays, these numerical approaches have become more and more computationally tractable even for highly non-linear problems. A collection of complex aerial tasks are presented in [START_REF] Geisert | Trajectory generation for quadrotor based systems using numerical optimal control[END_REF]. When complex motion is sought, optimization is often the chosen approach and has proved its efficiency for generating high speed manoeuvres with quadrotors [START_REF] Mellinger | Minimum snap trajectory generation and control for quadrotors[END_REF], navigation in heavily cluttered environments [START_REF] Landry | Aggressive quadrotor flight through cluttered environments using mixed integer programming[END_REF], passing though narrow gaps [START_REF] Falanga | Aggressive quadrotor flight through narrow gaps with onboard sensing and computing using active vision[END_REF] and so on.

More recently, self-collision and simple occlusion avoidance tasks for a humanoid were incorporated in a quadratic optimization problem in [START_REF] Agravante | Visual servoing in an optimization framework for the whole-body control of humanoid robots[END_REF]. Visibility of the feature is handled by using an avoidance task as in [START_REF] Folio | A controller to avoid both occlusions and obstacles during a vision-based navigation task in a cluttered environment[END_REF]. A more precise occlusion avoidance formulation is proposed in [START_REF] Nicolis | Occlusionfree visual servoing for the shared autonomy teleoperation of dual-arm robots[END_REF] but uses a larger set of visual constraints.

In [START_REF] Kazemi | Randomized kinodynamic planning for robust visual servoing[END_REF] the authors explored a randomized kinodynamic hybrid path planning approach applied to a manipulator for finding a feasible path. It satisfies a great deal of constraints both in the image and in the joint space but takes a few minutes to return a solution.

Finally, optimal problems can be solved online and handled in a continuous control approach similar to MPC [START_REF] Michalska | Robust receding horizon control of constrained nonlinear systems[END_REF]. This control framework is very convenient since it merges both reactive control and planning approaches by integrating feedback into the optimal planning. Most relevant works consider optimization for generating effortless collision-free trajectories by predicting the obstacles motion using visual measurements [START_REF] Gao | Quadrotor trajectory generation in dynamic environments using semi-definite relaxation on nonconvex qcqp[END_REF] and for more aggressive autonomous target tracking with field of view constraints in terms of acceleration [START_REF] Thomas | Autonomous flight for detection, localization, and tracking of moving targets with a small quadrotor[END_REF]. Along a similar direction, this work considers a multi-objective optimization formulation with extensions for explicit occlusion and collision avoidance with either a down or a front-looking camera (see Fig. 1). In a previous work [START_REF] Penin | Vision-based minimum-time trajectory generation for a quadrotor uav[END_REF] we presented preliminaries for online minimum-time trajectory replanning under field of view constraints. In this work we will show that reactive manoeuvres can still be achieved without minimizing time and in the presence of more complex constraints. Besides, we improved the accuracy and the numerical stability of the gradients evaluation by using complex-step differentiation to improve the SQP convergence. Our method differs from previous works under three main aspects:

• soft occlusion avoidance formulation in the image space • efficient and reactive online replanning strategy • visual constraints independent of the camera/UAV configuration

The rest of this paper is organized as follows. In Section II, we introduce the system and how the problem is formulated as a Nonlinear Program. In Section III, we give the methodology of the vision-based 3D target tracking problem. In Section IV, we present our occlusion and collision avoidance approach. The implementation of the proposed approach is explained in Section V as well as the replanning strategy. Section VII presents the experimental results performed with a quadrotor. Finally, we discuss the results and draw some conclusions in Section VIII. As well-known (see, e.g., [START_REF] Mahony | Multirotor aerial vehicles: Modeling, estimation, and control of quadrotor[END_REF]), the quadrotor dynamics can then be modelled as

ṗ = v (1a) v = g - f m z B (1b) Ṙ = R[ω] × (1c) ω = J -1 ([J ω] × ω + τ ) ( 1d
)
where m is the robot mass, J ∈ R 3 is the inertia tensor, g ∈ R 3 the (constant) gravity acceleration in the world frame, and [•] × the usual skew-symmetric operator. Finally (f, τ ) ∈ R 4 are the total thrust and torques applied by the propellers, which can be expressed as a set of control inputs u in terms of the individual propeller thrusts u = (f 1 , f 2 , f 3 , f 4 ) ∈ R 4 according to the following linear expression

f τ =     1 1 1 1 0 -l 0 l -l 0 l 0 b -b b -b     u (2)
where l is the distance between the propellers and the COM, and b is the drag constant. The set of admissible control inputs is the box

U = [f m , f M ] 4 with 0 < f m < f M .
Let us also assume the robot to be equipped with an onboard camera whose fixed pose w.r.t. B is known from a preliminary calibration. The spherical projection of a target point in 3D w.r.t. the frame B is the bearing vector

β = R T (r -p) R T (r -p) = m m ∈ S (3)
where r ∈ R 3 is the position of the feature in the world frame and S is the surface of the unit sphere and m = (u, v, 1) is the image measurement from which β is computed. The field of view is then modelled as a cone (see Fig. 1).

B. Differential flatness

A powerful property that can simplify the optimization problem is differential flatness. A non-linear system is termed flat if it is possible to find an analytic expression of the system nominal state and control inputs in terms of a finite number of variables (called flat outputs) and their derivatives up to a finite order. Then, this transformation gives a linear equivalent description of the system. Introduced by Fliess [START_REF] Fliess | Flatness and defect of nonlinear systems: Introductory theory and examples[END_REF] in the nineties it has been widely used since then in optimization especially for underactuated systems such as unicycles and (more recently) quadrotors (see, e.g., [START_REF] Mellinger | Minimum snap trajectory generation and control for quadrotors[END_REF], [START_REF] Richter | Polynomial trajectory planning for quadrotor fiight[END_REF], [START_REF] Landry | Aggressive quadrotor flight through cluttered environments using mixed integer programming[END_REF] or [START_REF] Tognon | Position tracking control for an aerial robot passively tethered to an independently moving platform[END_REF]). Differential flatness can significantly simplify any optimization problem by avoiding any numerical integration of the system dynamics (and the need of introducing the dynamics as differential constraints in the optimization scheme). It is well-known that the quadrotor dynamics (1) are flat with flat outputs γ = (p, ψ) T ∈ R 4 , where ψ is the yaw angle from the usual roll/pitch/yaw decomposition of the rotation matrix R.

Indeed, under the assumption f > 0, one can find an invertible algebraic mapping of the form:

χ = φ χ (p, v, v, v, ψ, ψ) (4a) (f, ḟ , f , τ ) = φ u ( v, v, ... v , ψ, ψ, ψ) (4b)
The complete expression of ( 4) can be found in, e.g., [START_REF] Faessler | Differential flatness of quadrotor dynamics subject to rotor drag for accurate tracking of highspeed trajectories[END_REF], and it is not reported here for lack of space. For simplicity of notation, we indicate with σ = (p, v, v, v, ... v , ψ, ψ, ψ) the vector of all quantities appearing on the right side of (4), and with σ χ = (p, v, v, v, ψ, ψ) only those involved in (4a). Now, any smooth enough trajectory in the flat space is also a feasible trajectory for the system.

C. Trajectory representation

We adopted a finite parametrization of the search space (i.e., of the flat outputs) using piecewise polynomial functions in the B-spline form [START_REF] Boor | A practical guide to splines[END_REF]. Given a vector of control points θ = (p 1 , . . . , p np , ψ 1 , . . . , ψ n ψ ) ∈ R 3np+n ψ , and two (fixed) normalized knot vectors s p ∈ [0, 1] lp , s ψ ∈ [0, 1] l ψ , the flat output trajectories from generic time t and of duration T can be represented as:

           p(s) = np k=1 B dp k s -t T p k ψ(s) = n ψ k=1 B d ψ k s -t T ψ k , ∀s ∈ [t, t + T ], (5)
where B d k is the k-th B-spline positive basis function of degree d -1 with minimal compact support where d is the B-spline order. They can be computed recursively as described in [START_REF] Boor | A practical guide to splines[END_REF]. Given [START_REF] Mebarki | Nonlinear visual control of unmanned aerial vehicles in gps-denied environments[END_REF], in order to ensure state continuity and input boundedness, one has to guarantee Lipschitz continuity of v and ψ (and continuity of lower order derivatives). This condition can be met by using open-uniform distributions of l = n + d ordered knots (i.e. s i = 0, for i = 1, . . . , d, s i = 1, In this section, we address the case of tracking a moving 3D target and ensuring final visibility of the target at the unique image center. We assume that the relative pose of the target can be estimated using vision only as done in [START_REF] Thomas | Autonomous flight for detection, localization, and tracking of moving targets with a small quadrotor[END_REF] with some preliminary knowledge of the target model (e.g., the target radius in case of a sphere). We are also interested in allowing the quadrotor undertaking aggressive manoeuvres for reaching the target in near minimum-time conditions. The final goal is to reach a hovering state such that the target appears at the image center while keeping a safety distance from the target (see Fig. 3).

for i = n + 1, . . . ,

A. Multi-Objective Cost Function

For achieving the aforementioned behaviour we adopt multiobjective programming and minimize the following weighted cost function at running time t with weights w i

J = w 1 t+T t ... v (t) 2 dt + w 2 t+T t ψ(t) 2 dt + w 3 t+T t d(t) -R s 2 dt + w 4 t+T t v(t) dt (6)
The latter is divided into four parts: the first and the second terms minimize the norm of the snap and the norm of the yaw acceleration respectively. They are used for encouraging smoothness while still exploiting the quadrotor aggressiveness [START_REF] Mellinger | Minimum snap trajectory generation and control for quadrotors[END_REF]. The third term minimizes the error between the distance of the target to the camera and a value R s defined as a safety distance: the radius of the sphere centred on the target (see Fig. 3). Naturally, to achieve more aggressive trajectories, one can enforce this cost. The fourth term penalizes the path length. Indeed, for a circular symmetric target the quadrotor may fly around the target indefinitely without changing the image feature position, an effect arising from the residual available degrees of freedom and that one can prevent by adding this fourth term.

In [START_REF] Penin | Vision-based minimum-time trajectory generation for a quadrotor uav[END_REF] we minimized the completion time T to generate aggressive trajectories. Here, we are dealing with more complex constraints and we will show that the optimal trajectories Fig. 4. The visual constraints ( 7) is equivalent as ensuring that the angular position β of an image point P is lower than the field of view angle α/2 can still be aggressive. Moreover, we avoid the penalization of such a complex parameter and the re-evaluation of the Bspline bases at each solver iteration. The completion time T is then a fixed parameter and should be chosen as a rough upper bound of the time required to reach the final pose vicinity. If T is too small, the trajectories might be infeasible. If too high, the system may be less reactive.

Finally, by suitably weighing these costs we are able to trigger the expected behaviour and prevent the objectives from conflicting with each other. Considering a front-looking camera, [START_REF] Sheckells | Optimal Visual Servoing for Differentially Flat Underactuated Systems[END_REF] adopted the minimization of the distance between the target and the camera center which also played the role of limiting visibility loss. However, this approach penalizes the motion aggressiveness by not fully exploiting the image space while, in our case, the target is free to move away from the camera center in order to allow large rotations of the camera and therefore to allow large accelerations of the quadrotor. The adopted visibility constraints are shown below.

B. Visibility constraints

Defining e c as the camera optical axis in the frame B, namely, x B (or y B ) for a front-looking camera and -z B for a down-looking camera, the visibility constraint is written as

β T (t)e c ≥ cos(α/2), ∀s ∈ [t, t + T ] ( 7
)
where α is the angle of view of the camera and β is given by (3). Figure 4 shows that ( 7) is equivalent as constraining the feature bearing angle β but whose numerical evaluation is more complex. Instead of imposing equality constraints on the final position in order to guarantee visual convergence of the target to the camera center, we rather define a terminal constraint such that the feature bearing angle has to belong to a (smaller) angular area at the camera centre (see Fig. 3) with

β T (t)e c ≥ cos(γ/2), ∀s ∈ [t+T v , t+T ], 0 ≤ T v ≤ T (8)
where γ is the angle defining the circular region of convergence in the image. The time T v is a parameter that defines at what time the feature shall enter the vicinity region. It can be tuned to affect the convergence rate towards this region. Namely, a value closer to 0 will demand a longer activation of the constraint. T v will also vary depending on the camera orientation due to the quadrotor dynamics.

With the above definitions, one can encapsulate visibility constraints and visual convergence for any camera orientation in the optimal framework defined in the Section V.

IV. VISION-BASED OBSTACLE AND OCCLUSION AVOIDANCE

In this section, we describe the vision-based optimization scheme for avoiding collisions with static spherical obstacles of radius R occ and the occlusions they might generate with the tracked target. Collision-free trajectories with N static spherical obstacles of inflated radius R col > R occ at position r obs are generated using volumetric constraints such as

p -r obsi 2 > R col 2 i , i = 1, . . . , N (9)
The occlusion constraint can be modelled analogously to [START_REF] Geisert | Trajectory generation for quadrotor based systems using numerical optimal control[END_REF], but in the image space, as follows

β -β obsi 2 > a 2 1 , i = 1, . . . , N (10)
where β and β obs are the spherical projections of the target and the obstacle center respectively, and a 1 is the length of the semi-minor axis of the projected ellipse of the spherical obstacle of radius R occ in the image plane. Assuming knowledge of the size of an obstacle in 3D with coordinates

X = (X o , Y o , Z o ) in the camera frame one has [28] a 2 1 = R 2 occ 4(Z 2 o -R 2 occ ) (11)
With this constraint one seeks to prevent the target from colliding the obstacles in the image space by keeping a minimum distance a 1 that grows as the depth Z o of the obstacle in the camera frame decreases. As the quadrotor moves towards the target the occlusion constraints from the obstacles passing behind the camera are of course discarded. However, when dealing with occlusion avoidance with a quadrotor there exist configurations where strict avoidance is not feasible. Indeed, when the target goes exactly below an obstacle (for a downlooking camera, see Fig. 5) the quadrotor may not have sufficient actuation capability or sufficient space for avoiding any occlusions. These situations may occur for any camera orientation. Therefore, in order to avoid such critical situations and always provide a feasible solution, we introduce a slack variable λ within the occlusion constraint to authorize partial occlusion if necessary. One has

a 2 1 = (R occ -λ) 2 4(Z 2 o -(R occ -λ) 2 ) , s.t 0 ≤ λ ≤ R occ (12)
The λ term plays the role of relaxing a hard constraint when the solver encounters not feasible situations. At the most critical configurations when the target is below an obstacle (Fig. 5), the collision constraint can be reduced or even cancelled by having λ reach the value R occ which is the actual radius of the obstacle (lower than R col ). With equation [START_REF] Falanga | Aggressive quadrotor flight through narrow gaps with onboard sensing and computing using active vision[END_REF] we observed very reactive responses from the system in case of occlusions. This may be due to the fact that small changes of λ induce a strong action on the occlusion constraint [START_REF] Mellinger | Minimum snap trajectory generation and control for quadrotors[END_REF]. Besides, by imposing a straightforward upper bound (R col ) for λ, its action will take effect only in case of violations of constraint [START_REF] Mellinger | Minimum snap trajectory generation and control for quadrotors[END_REF]. In case of strong and sudden occlusions when the target is below an obstacle, the red trajectory cannot be a viable solution since it may violate actuation and/or spatial constraints. We instead allow minimal constraint violations to keep the solver efficient. The blue trajectory is then a relaxed solution where λ hits its limit value Rocc (obstacle radius). The orange trajectory represents the case of a less relaxed occlusion constraint where λ eventually reaches a smaller value and is zero when there are no occlusions (green trajectory) It would be of course possible to use the ellipse semi-major axis defined as:

a 2 2 = R 2 occ (X 2 o + Y 2 o + Z 2 o -R 2 occ ) 4(R 2 occ -Z 2 o ) 2 (13)
However, in practice one has

X 2 o + Y 2 o Z 2
o especially if one considers the limited camera field of view. This is why we consider the semi-minor axis [START_REF] Landry | Aggressive quadrotor flight through cluttered environments using mixed integer programming[END_REF] which is also less complex. Of course, the occlusion constraint and the slack variable λ introduce conservatism to some extent. However, the main objective of λ is to improve stability and continuity of the solution in case of critical configurations.

As it is common in MPC, the trajectories are sampled. The discretization step should be chosen regarding the obstacles resolution and the trajectories smoothness since the motor thrusts can be very sharp with the considered B-Spline minimal orders.

V. OPTIMIZATION PROBLEM DEFINITION Finally, the current optimization problem can be stated as the following static NLP Problem 1: find θ, λ such that:

min θ,λ J + w 5 λ 2 (14a) s.t. σ χ (t) = σ χt , (14b)
σ (i) (t + T ) = σ χ * , i = 1, . . . , 3, (14c) u(s) ∈ U, ∀s ∈ [t, t + T ], (14d)
β T (t)e c ≥ cos(α/2), ∀s ∈ [t, t + T ], (14e)
β T (t)e c ≥ cos(γ/2), ∀s ∈ [t + T v , t + T ] (14f) p(s) -r obs 2 > R 2 col , ∀s ∈ [t, t + T ], (14g)
β(s) -β obs 2 > a 2 1 , ∀s ∈ [t, t + T ], (14h) 0 ≤ λ ≤ R occ (14i)
where σ χt = φ -1 χ (χ t ) and σ χ * = φ -1 χ (χ * ). For a final hovering state χ * one has of course σ (i) (t+T ) = 0, i = 1, . . . , 3. We choose not to explicitly constrain the final cartesian position which is considered as a free parameter to be determined by the optimization algorithm. We seek feasible trajectories with constraint (14d) but, due to unknown target motion, the quadrotor might not have sufficient actuation to always satisfy the visibility constraint (14e).

The costs and constraints gradients in the flat space are computed analytically. However, it is more efficient to estimate the gradients of the visual and inputs constraints due to their complex and heavy analytic formulation. More especially, we use complex-step differentiation [START_REF] Squire | Using complex variables to estimate derivatives of real functions[END_REF]. By using the following first-order approximation

∂f ∂x = Im{f (x + ih)} h + O(h 2 ), h ∈ R ( 15
)
This technique is attractive to determine first derivatives since it only requires a single evaluation of the function and avoids the problem of subtractive cancellation of classic finite approximation (round-off errors). Therefore, it is known to have superior accuracy (close to analytic accuracy) and numerical stability as analysed in [START_REF] Abreu | On the generalization of the complex step method[END_REF]. It is also less intrusive in terms of program transformation than automatic differentiation which can sometimes require large and deep source overloads. On the other hand, complex differentiation requires some mathematical adaptations to be used with complex values (the square root or the absolute value functions for instance) which increases the computational cost. The choice between these two methods then hinges on a trade-off between ease of implementation and execution efficiency and is further discussed in [START_REF] Martins | The connection between the complex-step derivative approximation and algorithmic differentiation[END_REF]. At this point, any off-the-shelf optimization solver can be used to find a numerical solution to Problem 1. Unfortunately, due to the non trivial non-linearities in its formulation (e.g.,(14d), (14h)), Problem 1 cannot be proven to be convex. The optimization will thus, in general, return a local minimum.

A. Trajectory replanning Strategy

In order to be reactive to the target unknown motion and increase robustness to the different sources of disturbance (e.g., noise, miscalibrations, neglected dynamics, and so on) that will make the robot to diverge from the planned trajectory, we take inspiration from MPC. Namely, we perform an online re-planning of an optimal trajectory by solving Problem 1 to incorporate a feedback action.

SQP can solve non-convex problems relatively fast by iterating from a good initial guess. For that, we use the same replanning process detailed in our previous work [START_REF] Penin | Vision-based minimum-time trajectory generation for a quadrotor uav[END_REF] where we discussed how to instantly generate an initial guess by operations on the B-Splines only. This method exploits two important properties of B-Splines. The local support property implies that the shape of the curve in a knot span (s k , s k+1 ) is only determined by a subset of the B-spline control points. The convex hull property guarantees, instead, that in each knot span, the spline curve is locally contained in the convex hull of the same subset of control points.

Using the proposed B-Spline adaptation, we achieve faster convergence to a local optimum and provide more stability and continuity of the initial guess than using random search techniques for example.

B. Down-looking case

In this section, we express the replanning strategy by solving Problem 1 for tracking a target on the ground and considering a down-looking camera with a (more realistic) pyramid-shaped field of view as in [START_REF] Penin | Vision-based minimum-time trajectory generation for a quadrotor uav[END_REF]. The perspective projection of the target is considered

β = R T (r -p) z B T R T (r -p) =   u v 1   ∈ P 2 (16)
where P 2 is the space of 3-D homogeneous vectors.

We adopt here the minimization of the distance error by the penalization of the position error with the target in the XY plane e xy ∈ R 2 .

Moreover, we let the final height z T free but bounded for limiting the motion along z W . This lighter formulation will basically compute trajectories for steering the quadrotor to a final hovering position such that the target appears at the camera center.

J 2 = w 1 T t ... v (t) 2 dt + w 2 T t ψ(t) 2 dt + w 3 T t e xy (t) 2 dt + w 4 z T + w 5 λ 2 (17)
We define the following optimal problem Problem 2: find θ, λ such that:

min θ,λ J 2 (18a)
s.t. σ χ (t) = σ χt , (18b)
σ (i) (t + T ) = 0, i = 1, . . . , 3, (18c)
p xy (t + T ) = p * xy , (18d)
ψ(t + T ) = ψ * , (arbitrary value) (18e) β(s) ∈ Ω, ∀s ∈ [t, t + T ], (18f) u
(s) ∈ U, ∀s ∈ [t, t + T ], (18g)
p(s) -r obs 2 > R 2 col , ∀s ∈ [t, t + T ], (18h)
β(s) -β obs 2 > a 2 1 , ∀s ∈ [t, t + T ], (18i) 0 ≤ λ ≤ R occ , (18j)
Z min ≤ z T ≤ Z max (18k)
Where p * xy is the target position in the XY plane at time t and the (square

) image domain Ω = {β ∈ P 2 s.t. max(|β T x B |, |β T y B |) ≤ tan(α)}.
It should be noted that the gradient of the occlusion constraint (18f) gives three possible directions for the decision variables to satisfy this constraint: 1) increase the distance between β(s) and β obs , 2) reduce the quantity a 1 by increasing Z o , 3) relax the constraint by increasing λ. In order to restrain the second direction that steers the solution towards the local minima of infinite height (a 1 → 0) we put an upper limit to the quadrotor final height z T with an additional constraint (18k). We also minimize z T (in θ) to prevent the camera from staying at the maximum height Z max . The constraint also helps reducing unnecessary and tricky motion towards the opposite direction by considering a minimum height Z min . This formulation allows efficient replanning at the rate of 30Hz. The next section gives more insight on the different system behaviours with respect to the given weights.

VI. SIMULATION RESULTS

The presented approach was validated in a physically realistic environment. The inertial parameters of the quadrotor were slightly biased on purpose to introduce model uncertainties. Then, since the controller cannot perfectly track the desired trajectories the replanning strategy played an important role for re-adapting the reference trajectory online. The quadrotor dynamics were simulated using V-Rep1 at 150Hz. The replanning strategy described in Sect. V-A was implemented in C++ within ROS and we exploited the SQP optimization routine from NLOPT [START_REF] Johnson | The nlopt nonlinear-optimization package[END_REF] to compute the optimal solutions. The generated trajectories of the flat outputs were sent to TeleKyb [START_REF] Grabe | The TeleKyb Framework for a Modular and Extendible ROS-based Quadrotor Control[END_REF] which then computed the actual control inputs using controller [START_REF] Lee | Geometric tracking control of a quadrotor UAV on SE(3)[END_REF].

The tracked target shown in Fig. 5 is manually commanded in velocity saturated at 1.0 m/s. We considered three spherical obstacles of radius R occ = 0.15m and of inflated radius R col = 0.4m (see Figure 8).

We invite the reader to refer to the attached video for a better insight of the results and to a higher-resolution video at https://www.youtube.com/watch?v=mvvF1I72HM8

Table I shows the solver parameters used for computing the optimal trajectories from Problem 2, respectively: the number of computed trajectories, the duration T max for the stopping criterion, the planning horizon T , the range of SQP iterations and the mean iteration, the step and constraint tolerance and finally the percentage of total optimization instances that failed to return a solution within T max . When it occurs the previous solution is returned. The evolution of λ and the the final height z T during the simulation is show in Fig. 6 and Fig. 7 respectively. The In the proposed formulations, the weights in the cost function are free parameters and need to be tuned empirically in order to generate the expected behaviour. Since the optimal problem is complex and highly nonlinear, any set of weights can generate a different response. However, from our experience the tuning triggers consistent responses and does not necessitate a particularly fine analysis. Generally, a higher w 1 will produce smoother trajectories and less aggressive motion. A higher w 3 will increase the convergence rate towards the target and a lower w 5 will tend to produce more relaxed occlusion constraints. Table II shows the numerical parameters used for the simulation. Moreover, we have used five different sets of weights (Table III) an discuss their consistent effect on the generated trajectories shown in Fig. 8. III) in our simulation environment. The camera view is shown in the upper right corner with the image trajectories and a view from above the target is shown below. The orange trajectory uses the chosen values (in Table II) with a higher w 3 . One can observe a faster convergence towards the image center (with an overshoot). The blue one has a higher w 1 resulting in a smoother trajectory than the orange one. The green one is even smoother and is more compliant to occlusions (less jerky) since w 5 is lower. The yellow trajectory is very sharp (low w 1 ) because the occlusion constraints are more respected (high w 5 ) and the final height is less penalized (low w 4 ). The pink solution gives a different and smoother path (high w 1 ) that benefits less from the constraints relaxation (since w 5 is high).

VII. EXPERIMENTAL RESULTS

This section reports the experiment meant to validate the proposed replanning strategy applied to the tracking of a target with the visibility constraints detailed in Sect. III-B (the reader is also referred to the attached video). For the experiment illustrated in this section we used a MK-Quadro from MiKroKopter equipped with a front-looking camera with a field of view of 90 • . The setup included an on-board ODROID-XU4 Linux Computer running ROS and the TeleKyb framework for interfacing the replanning algorithm ran on a standard desktop PC (Quadcore Intel i7 CPU@2.6 GHz). A Vicon motion capture system was employed for giving state measurements of the quadrotor and position measurements of the tracked target. We used AprilTags as a generic target attached on top of a other MK-Quadro controlled remotely (Fig. 9).

VIII. CONCLUSION AND FUTURE WORK

In this work, we adapted our previous approach [START_REF] Penin | Vision-based minimum-time trajectory generation for a quadrotor uav[END_REF] on vision-based optimal trajectory generation to a wider context by considering reactive target tracking and both occlusion and obstacle avoidance for either a front or a down looking camera (or any other camera orientation).

The quadrotor trajectories are mainly driven by vision while seeking aggressive but smooth trajectories that respect actuation and sensor limits for any camera orientations. Then, starting from a good initial guess the solver was able to return an optimal solution about 94% of the time within 33ms allowing an online replanning strategy capable of absorbing noise, disturbances and any non-modelled effect for long duration flights. The same strategy was applied during an experiment using a real quadrotor for the case of a frontlooking camera.

In contrast with the previous literature, we proposed a method that explicitly handles visibility constraints and occlu-sion avoidance within a fast online replanning strategy. The occlusion constraint is expressed as a single constraint per object in contrast with [START_REF] Nicolis | Occlusionfree visual servoing for the shared autonomy teleoperation of dual-arm robots[END_REF] and could be applied to non-spherical objects to some extent by finding the enclosing sphere. For elongated objects such as bars, one could cover the object main axis with a finite number of spheres as it is done sometimes for collision avoidance with robot manipulators (known as "sphere swept models" see e.g., [START_REF] Sugiura | Real-time collision avoidance with whole body motion control for humanoid robots[END_REF]). Thus, there is a tradeoff between conservatism/accuracy and ease of computation. Finally, we coped with the issue that the discontinuity of occlusion constraints can generate by introducing a minimal relaxation in Sect. IV.

To be successful, our path-planning approach requires a complete knowledge of the environment and robot model. These requirements can be limiting in many real applications. The need for such exact knowledge could be relaxed by accounting for modelling and calibration uncertainties at planning stage. Finally, future work includes validating the method with real (and maybe dynamic) obstacles and the use of vision only for estimating the target relative position. In this case the relaxation term will play an even more decisive role.

Fig. 2 .

 2 Fig. 2. Quadrotor model II. PRELIMINARIES A. Quadrotor dynamics With reference to Fig. 2, let us define a world frame W and a body frame B with origin O B attached to the center of mass (COM) and axis z B parallel to the propeller rotational axes. Let us also assume, without loss of generality, that the robot COM corresponds to the barycentre of the propellers. The robot state is χ = (p, R, v, ω) ∈ SE(3) × R 6 where p ∈ R 3 is the position of the robot COM in W, R ∈ SO(3) is the rotation matrix from W to B, v the COM linear velocity expressed in W and ω the angular velocity expressed in B.As well-known (see, e.g.,[START_REF] Mahony | Multirotor aerial vehicles: Modeling, estimation, and control of quadrotor[END_REF]), the quadrotor dynamics can then be modelled as

  l, and s d , . . . , s n+1 equally spaced in [0, 1]) and by taking d = d p = 4 for p and d = d ψ = 2 for ψ.

Fig. 3 .

 3 Fig.3. The quadrotor has to follow the target in red while keeping a safety distance represented by the light blue sphere. The terminal constraint in the image space is represented by the green circle. The blue circle represents the field of view inside which the feature trajectory must lie. Here, the optimal trajectory in blue steers the quadrotor towards a final position (white dot) and the resulting image feature trajectory is the one in white

Fig. 5 .

 5 Fig.[START_REF] Cowan | Visual servoing via navigation functions[END_REF]. In case of strong and sudden occlusions when the target is below an obstacle, the red trajectory cannot be a viable solution since it may violate actuation and/or spatial constraints. We instead allow minimal constraint violations to keep the solver efficient. The blue trajectory is then a relaxed solution where λ hits its limit value Rocc (obstacle radius). The orange trajectory represents the case of a less relaxed occlusion constraint where λ eventually reaches a smaller value and is zero when there are no occlusions (green trajectory)

Fig. 6 .

 6 Fig.6. Evolution of λ for the simulation. The peaks reveal the presence of sudden occlusions. The flat section (from around instance 1100 to instance 1700) respresents a flight period far from the obstacles. Therefore, λ is very close to 0. Here λ did not hit the upper limit Rocc = 0.15m.

Fig. 7 .

 7 Fig. 7. Evolution of z T for the simulation. The final height was bounded between 1.9m and 3.5m. Again, z T reaches its minimum value when no collisions and occlusions occur.

Fig. 8 .

 8 Fig. 8. Different solution trajectories returned with different sets of weights (in TableIII) in our simulation environment. The camera view is shown in the upper right corner with the image trajectories and a view from above the target is shown below. The orange trajectory uses the chosen values (in TableII) with a higher w 3 . One can observe a faster convergence towards the image center (with an overshoot). The blue one has a higher w 1 resulting in a smoother trajectory than the orange one. The green one is even smoother and is more compliant to occlusions (less jerky) since w 5 is lower. The yellow trajectory is very sharp (low w 1 ) because the occlusion constraints are more respected (high w 5 ) and the final height is less penalized (low w 4 ). The pink solution gives a different and smoother path (high w 1 ) that benefits less from the constraints relaxation (since w 5 is high).

Fig. 9 .

 9 Fig. 9. The center of the box (yellow dot) attached on a second quadrotor and covered by four ApriTags had to remain inside the field of view of the camera (red circle) during the entire flight and appear inside the blue circle in the end of the computed trajectories

TABLE II CONSIDERED

 II PARAMETERS FOR THE SIMULATION

TABLE III VALUES

 III OF THE WEIGHTS USED FOR GENERATING THE TRAJECTORIES SHOWN IN FIG. 8

http://www.coppeliarobotics.com/

' comments. This work was supported by the French Direction Générale de l'Armement (DGA).


M Querre

J Tranchant

B Corraze

Stéphane Cordier

V Bouquet

Stephanie Deputier

Maryline Guilloux-Viry

M. -P Besland

Etienne Janod

Laurent Cario
email: laurent.cario@cnrs-imn.fr

M Querré

S Députier

M.-P Besland

Non-volatile resistive switching in the Mott insulator (V1-xCrx

Keywords: Mott insulators, V2O3, Mott transition, resistive switching, non-volatile memory, ReRAM

Introduction

Resistive Random Access Memories (ReRAM) are considered as one of the most promising technologies for future non-volatile memories. They are based on resistive switching mainly related to two physical mechanisms: ionic electromigration (anionic for OxRAM and cationic for CBRAM) and thermochemical reactions (typically in NiO) [1,2,3]. But recently the insulator to metal transition (IMT) in Mott insulators has focused much attention due to its serious potential as an alternative mechanism for ReRAM [4]. Mott insulating state can be usually destabilized by pressure and doping into a correlated metal phase via an isostructural electronic phase transition, as shown for example in the temperature-pressure phase diagram of V2O3 in Fig. 1(a) [5]. However recent works [4,6,7,8] have shown that electric field can also destabilize the Mott insulating state and trigger a non-volatile resistive switching in these materials. Studies performed on the chalcogenide compounds AM4Q8 (A=Ga, Ge; M=Nb, V, Ta; Q=S, Se) have demonstrated that this resistive switching (RS) originates from an avalanche breakdown phenomenon [7] which causes the collapse of the paramagnetic Mott insulator state at the nanoscale [10] and results in the creation of granular conductive filaments [11].

This discovery on the chalcogenide AM4Q8 compounds raises the question of the universal nature of this new type of non-volatile resistive switching. We have therefore sought such resistive switching in the most famous canonical narrow gap Mott insulator, namely (V1-xCrx)2O3. We report here that this oxide Mott insulator exhibits a non-volatile resistive switching of the same physical nature as the AM4Q8 compounds. Moreover we show that the performances of a (V1-xCrx)2O3 memory cell obtained on a device mostly compatible with microelectronics standards are very promising.

Material and methods

(V1-xCrx)2O3 crystals synthesis

V2O5 (Aldrich, > 99.6 %) powder was first dried at 400°C , then mixed with Cr2O3 (Prolabo, 99.9 %) in the 97.5 : 2.5 molar ratio, and then heat-treated at 900°C for 10h under a 95% Ar -5% H2 gas flow. Half a gram of the obtained (V0.975Cr0.025)2O3 powder was then placed in a silica tube together with 40 mg of sulfur as a vapor phase transport agent. The tube was vacuum sealed, heated up to 1050°C in a temperature gradient (≈ 10°C/cm) furnace and then cooled down at -2°C/h to 900°C, before a faster cooling (-300°C/h) to room temperature. Energy Dispersive X-ray Spectroscopy analyses carried out on the obtained crystals (typical size 100 ´ 50 ´ 300 µm) indicate that the V / Cr ratio is close to the nominal value and show the absence of sulfur impurities in the crystals.

Thin film synthesis and device preparation

Thin films of (V0.95Cr0.05)2O3 were deposited by reactive co-sputtering of V and Cr targets in Ar/O2 discharge [12], leading to thicknesses in the 100-880 nm range. Post deposition annealing at 500°C in a reducing atmosphere yields crystallized and single-phased (V0.95Cr0.05)2O3 layers as checked by X-Ray diffraction (XRD). In V2O3, deviation from the ideal V/O ratio by 1.75 % lowers the magnetic ordering temperature (paramagnetic metal PM to antiferromagnetic insulator AFI transition, see Fig. 1(a)) from 155 K to 11 K [13]. But our (V0.95Cr0.05)2O3 thin films deposited on SiO2/Si substrates and measured with four probes in a planar configuration exhibit a paramagnetic Mott insulator PMI -AFI transition temperature (158 K) very close to the expected one (175K), which indicates that they are almost stoichiometric. This prevents spurious electromigration processes associated to ionic diffusion and excludes potential OxRAM/CBRAM mechanisms. Devices based on a symetric TiN / (V0.95Cr0.05)2O3 / TiN MIM structure were subsequently realized with via diameters in the 330 -1600 nm range.

Electrical measurements

(V0.975Cr0.025)2O3 crystals were contacted with 17 μm gold wires using a carbon paste (Electrodag PR-406) and then annealed in vacuum at 150 °C for 20 min to achieve a low contact resistance. The low-bias resistance was measured using a high-impedance Keithley 236 source-measure unit. Thin film devices and crystals were connected in series with a load resistance Rload playing the role of current limiter and ranging in the 10-80 % of the device or crystal resistance. The write-erase cycling measurements were performed using an Agilent 81150A pulse generator.

Results and Discussion

Resistive switching in (V1-xCrx)2O3 crystals

Figure 1(b) shows the resistance versus temperature (blue curve) measured on a (V0.975Cr0.025)2O3 single crystal grown by the sulfur vapor transport method. Above 190K, the crystal shows an insulating behavior consistent with the expected PMI state. The transition from the PMI to the AFI phase occurs around 180 K which is again in good agreement with the value expected from the phase diagram (see Figure 1(a)). The crystal in series with a load resistance was subsequently subjected to an electric pulse (75V, 100µs at 155K). As shown in Figure 1(b), the application of this electric pulse induces a non-volatile decrease of the resistance R(T) by several orders of magnitude at low temperature (red curve). This observation reveals that the application of short electric pulses can induce a non-volatile metal-insulator transition in the Mott insulator (V0.975Cr0.025)2O3 which is reminiscent of the behavior of AM4Q8 compounds [14].

Previous in-depth studies of GaTa4Se8 have unveiled important features of the non-volatile resistive switching of the AM4Q8 family of compounds. This narrow gap Mott insulator exhibits indeed an insulator to metal transition under pressure (at about 5 GPa) and becomes superconducting below TC= 4-7K above 11 GPa [11,12]. In GaTa4Se8, Scanning Tunneling Microscopy (STM) measurements have revealed that the non-volatile resistive switching is related to an electronic phase separation at the nanoscale [10,11]. The surface topography of pristine GaTa4Se8 crystals is structureless, but after resistive switching some granular filamentary structures qualitatively oriented along the direction of the electric pulses appear on the surface. Spectroscopic (STS) experiments have confirmed the existence of electric-field-induced metallic nano-domains of about 30-70 nm within these granular filamentary paths. Moreover the presence of granular superconductivity at low temperature, reminiscent of the one observed under pressure, directly suggests the presence of compressed metallic nano-domains in transited crystals of GaTa4Se8 [11].

This hypothesis was further tested in the case of (V1-xCrx)2O3 using the transition towards the AFI phase as a probe. The phase diagram of this compound (see Fig. 1a) shows indeed that the transition temperature from the Mott insulating phase (PMI) to the AFI phase does not change much with composition or pressure. Conversely the huge variation observed under composition/pressure for the transition temperature from the metallic (PM) to the AFI phase may be used to characterized the formation under electric pulse of compressed metallic domains. conducting path (light blue regions in Fig. 1 (c)). However a new transition is clearly observable around 40K in the R(T) curve of the switched (V0.975Cr0.025)2O3 crystal (in red). According to the phase diagram of this system, this may be viewed as the PM ® AFI transition of (V0.975Cr0.025)2O3 under a pressure of 30 kbar [17]. This is a clear indication of the creation of compressed metallic domains under electric field. After resistive switching the metallic domains would therefore correspond to (V0.975Cr0.025)2O3 domains under pressure, as indicated by the red dotted line in Fig. 1(a). From these experiments, we may therefore infer, as sketched in Figure 1(c), that the Mott insulator compound (V1-xCrx)2O3 that has experienced a resistive switching under electric field contains some conducting granular filamentary paths embedded in a pristine like matrix. Moreover these granular filamentary paths include some compressed regions that have undergone a Mott insulator to metal transition. . It shows that reversible resistive switching may be obtained on crystal at room temperature using a series of voltage pulses (+90V for the SET and -50V for RESET). This observation is consistent with data obtained on AM4Q8 [6]. The RESET transition was explained on the basis of a modeling work by the transition back of the metallic domains to the insulating state thanks to a thermal activation process enabled by self-Joule heating of the filament during the pulse [7,16]. This reversible resistive switching opens the way to the use of the Mott insulator (V1-xCrx)2O3 as an active material of a non-volatile ReRAM memory. However the ROFF/RON ratio on crystals is quite poor (< 1.05), ROFF and RON being defined respectively as the resistance values of high (OFF) and low (ON) resistance states. This is related to the large section of the crystal compared to the size of the conducting filamentary path embedded in this crystal. To enhance the ROFF/RON ratio it is therefore important to design new devices with smaller electrode areas.

Reversible resistive switching in (V1-xCrx)2O3 crystals

Electrical performances of (V1-xCrx)2O3 thin film devices

To further explore the potential for memory application, some devices in Metal-Insulator-Metal (MIM) structure with small electrode areas (vias) were prepared from well crystallized thin film of (V1-xCrx)2O3 Mott insulators (see Figure 3). As displayed in Fig. 4(a) in the case of a 880 nm thick (V0.95Cr0.05)2O3 film, the application of 3.5 V / 1.5 µs (SET transition from the high to the low resistance state) and 4.5 V / 250 ns (RESET transition) pulses to such device with 500 nm vias engenders resistive switching cycles with an average ROFF/RON close to 20. This ratio evolves linearly with the inverse of the electrode surface 1/S, since vias of 330 nm and 1.6 µm diameter lead to average ROFF/RON ratios of 50 and 2, respectively. Moreover, no significant effect of layer thickness on the ROFF/RON ratio was noticed, as ROFF and RON values increase in the same proportions with film thickness. This scaling of ROFF/RON ratio with 1/S is conform to previous observations on GaV4S8 thin films [4] and is consistent with the quasi-dissolution of a granular metallic filament in a Mott insulating matrix. Higher ROFF/RON ratios are then expected for electrode sizes lower than 330 nm. It must be noted that the resistive switching cycles obtained on a very thick layer (880nm, see Fig. 4 confirms previous measurements on GaV4S8 thin films [2,17]. It has to be noted that resistive switching cycles could be obtained on the device after 12 days of retention test. Lastly, preliminary endurance results show that these devices can be cycled more than 1000 times.

Conclusions

We have demonstrated non-volatile resistive switching in the canonical Mott insulator (V0.95Cr0.05)2O3. Experiments conducted on crystals reveal that this resistive switching is consistent with the formation of granular conducting filaments under electric field. These filaments contain metallic domains closely related to the correlated metallic phase observed under high pressure in the P-T phase diagram of (V1-xCrx)2O3. It confirms therefore that the resistive switching in (V1-xCrx)2O3 is related to an electric field driven Mott transition at the local scale as observed previously in chalcogenide Mott insulators.

Memory cell devices based on the Mott insulator oxide (V0.95Cr0.05)2O3 were fabricated with an excellent crystalline quality of the layers. The demonstrated performances are competitive in terms of programming voltage / time. Moreover the increase of ROFF/RON ratio with miniaturization, and the preliminary retention / endurance tests are very promising. These results definitely set (V1-xCrx)2O3 based Mott memories as promising candidates among ReRAM solutions.

Figure 1 (

 1 Figure1(b) shows the resistance versus temperature (blue curve) measured on a (V0.975Cr0.025)2O3 single crystal grown by the sulfur vapor transport method. Above 190K, the crystal shows an insulating behavior consistent with the expected PMI state. The transition from the PMI to the AFI phase occurs around 180 K which is again in good agreement with the value expected from the phase diagram (see Figure1(a)). The crystal in series with a load resistance was subsequently subjected to an electric pulse (75V, 100µs at 155K). As shown in Figure1(b), the application of this electric pulse induces a non-volatile decrease of the resistance R(T) by several orders of magnitude at low temperature (red curve). This observation reveals that the application of short electric pulses can induce a non-volatile metal-insulator transition in the Mott insulator (V0.975Cr0.025)2O3 which is reminiscent of the behavior of AM4Q8 compounds[14].Previous in-depth studies of GaTa4Se8 have unveiled important features of the non-volatile resistive switching of the AM4Q8 family of compounds. This narrow gap Mott insulator exhibits indeed an insulator to metal transition under pressure (at about 5 GPa) and becomes superconducting below TC= 4-7K above 11 GPa[11,12]. In GaTa4Se8, Scanning Tunneling Microscopy (STM) measurements have revealed that the non-volatile resistive switching is related to an electronic phase separation at the nanoscale[10,11]. The surface topography of pristine GaTa4Se8 crystals is structureless, but after resistive switching some granular filamentary structures qualitatively oriented along the direction of the electric pulses appear on the surface. Spectroscopic (STS) experiments have confirmed the existence of electric-field-induced metallic nano-domains of about 30-70 nm within these granular filamentary paths. Moreover the presence of granular superconductivity at low temperature, reminiscent of the one observed under pressure, directly suggests the presence of compressed metallic nano-domains in transited crystals of GaTa4Se8[11].This hypothesis was further tested in the case of (V1-xCrx)2O3 using the transition towards the AFI phase as a probe. The phase diagram of this compound (see Fig.1a) shows indeed that the transition temperature from the Mott insulating phase (PMI) to the AFI phase does not change much with composition or pressure. Conversely the huge variation observed under composition/pressure for the transition temperature from the metallic (PM) to the AFI phase may be used to characterized the formation under electric pulse of compressed metallic domains. Figure 1(b) compare the R(T) curves of the (V0.975Cr0.025)2O3 crystal in the pristine state (blue curve) and after resistive switching (red curve). Both curves show the signature of the PMI-AFI transition around 180K; in the switched crystal, it corresponds to the transition of the (V0.975Cr0.025)2O3 matrix embedding the filamentary

Fig. 1 :

 1 Fig. 1 : (a) Theoretical phase diagram of (V1-xCrx)2O3 [5]. (b) Evolution of resistance vs temperature of a (V0.975Cr0.025)2O3 crystal in its pristine state and in a transited state. As expected the PMI-AFI phase transition occurs around 180 K in the pristine state (see blue dotted line in (a) and blue arrow in (b)). In the transited state a new transition occurs at 40 K suggesting the formation of compressed metallic domains (see red dotted line in (a) and red arrow in (b)). (c) Schematical drawing of the formation of a granular filamentary path in the transited state of a (V0.95Cr0.05)2O3 crystal containing some compressed domains with a PM-AFI phase transition around 40 K as indicated by the red dotted line in (a).

Fig. 2 :

 2 Fig. 2: Resistive switching cycles obtained at room temperature on large crystals of (V0.925Cr0.075)2O3 by applying a series of pulses (+90V for the SET and -50V for RESET).

Figure 2

 2 Figure 2 reveals another interesting feature of the non-volatile resistive switching in (V1-xCrx)2O3. It shows that reversible resistive switching may be obtained on crystal at room temperature using a series of voltage pulses (+90V for the SET and -50V for RESET). This observation is consistent with data obtained on AM4Q8[6]. The RESET transition was explained on the basis of a modeling work by the transition back of the metallic domains to the insulating state thanks to a thermal activation process enabled by self-Joule heating of the filament during the pulse[7,16]. This reversible resistive switching opens the way to the use of the Mott insulator (V1-xCrx)2O3 as an active material of a non-volatile ReRAM memory. However the ROFF/RON ratio on crystals is quite poor (< 1.05), ROFF and RON being defined respectively as the resistance values of high (OFF) and low (ON) resistance states. This is related to the large section of the crystal compared to the size of the conducting filamentary path embedded in this crystal. To enhance the ROFF/RON ratio it is therefore important to design new devices with smaller electrode areas.

Fig. 3 :

 3 Fig. 3: (a) SEM cross-sectional picture and (b) XRD diagram of a 880 nm thick (V0.95Cr0.05)2O3 film. (c) Top view of a substrate used for deposition of (V0.95Cr0.05)2O3 thin films and (d) SEM picture of a bottom electrode with a via of 500 nm in diameter. (e) Cross-sectional schematic view of the device based on (V0.95Cr0.05)2O3 and (f) corresponding back scattered electron SEM picture of the multi-layer.

  (a)) without any initial forming definitely exclude any potential OxRAM mechanism, since the diffusion of oxygen vacancies on such distances is not realistic. Temporal retention of OFF and ON resistance states presented in Fig. 4(b) shows a very good stability of both states at room temperature. The slight increase of ROFF/RON ratio with time

Fig. 4 :

 4 Fig.4: (a) Resistive switching cycles obtained at room temperature on devices with 880 nm thick films of (V0.95Cr0.05)2O3 and 500 nm vias, by applying 3.5V / 1.5 µs pulses for the SET and 4.5V / 250 ns pulses for RESET. (b) Data retention at room temperature, extrapolated to 10 years, for OFF (red circles) and ON (blue squares) states on a 880 nm thick film of (V0.95Cr0.05)2O3.

Acknowledgements :

The authors acknowledge SATT Ouest-Valorisation, CNRS and "Region Pays de la Loire" for funding.