M Hechelski

M Saulais

Marc Danjean

Laurent Lassabatère

Guillaume Attard

Jean-Philippe Bedell

Occurrence, spatial distribution and succession of plants in contaminated sediments in urban stormwater basins

à 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.

Occurrence, spatial distribution and succession of plants in contaminated sediments in urban stormwater basins

Détermination, distribution spatiale et évolution dans le temps des plantes sur des sédiments contaminés de bassins d'infiltration urbains HECHELSKI Marie, SAULAIS Muriel, DANJEAN Marc, LASSABATERE Laurent, ATTARD Guillaume, BEDELL Jean-Philippe LEHNA -Equipe IPE, Ecole Nationale des Travaux Publics de l'Etat, UMR 5023, Université de Lyon, Rue Maurice Audin, Vaulx-en-Velin Cedex, 69518 ; E-mail: bedell@entpe.fr, Tel.:+33-4-72047081

RÉSUMÉ

Nous avons identifié sur 10 ans l'évolution de la flore en regardant la richesse et l'abondance des espèces dominantes sur deux types de bassins végétalisés: Minerve où la végétation a été planté volontairement, l'autre qui s'est naturellement colonisé (appelé Django). Pour Minerve, nous avons observé un gradient de l'amont (près de l'entrée de l'eau) à l'aval (partie la plus éloignée de l'entrée de l'eau) pour les paramètres comme la teneur en matière organique, la hauteur de sédiment déposé et les teneurs en éléments traces. Pour Django, la distribution est moins claire, mais nous avons aussi observé des zones avec des contenus contrastées en matière organique et en métaux lourds. La régression de la communauté des hélophytes sur Minerve au profit d'espèces rudérales pionnières marque l'évolution de ce site. Depuis 1999, plus de 93% des espèces implantées ont disparu, mais la diversité des hélophytes semble s'être stabilisée depuis 2008 avec seulement 12,5% de la variation du nombre de ces espèces. Sur le bassin de Django, le nombre d'espèces est resté stable depuis 2008. Mais la diminution de près de 50% du taux de recouvrement des macrophytes, réduit à Phalaris arundinacea, Typha latifolia et Schoenoplectus tabernamontanii, souligne un changement des débits d'eau entrants. La flore dominante trouvée dans des bassins d'infiltration évolue donc au fil du temps avec des espèces euryéces de pseudo-métallophytes des macrophytes (Typha latifolia, Phragmites australis, Iris pseudacorus, Phalaris arundincea), des espèces rudérales typiques des prairies sèches (Asteraceae, Poaceae, Brassicacées) et des espèces des zones humides (Rumex sp.). Avant tout, ces plantes sont principalement divisés/distribuée dans les bassins en fonction de facteurs hydrologiques, à savoir la disponibilité de l'eau.

ABSTRACT

We identified over a decade the evolution of species richness of flora and its abundance on two types of vegetated basins: Minerve where the vegetation was originally planted, and the other naturally colonized (called Django). In Minerve, we observed a gradient from upstream (near the water entry) to downstream (farthest from this water entry) for parameters like organic content, height of layer and trace element contents. In Django, the distribution is less clear, but we also observed zones with contrasting contents in organic matter and heavy metal. The regression of helophytes communities on Minerve for the benefit of ruderals pioneer species marks the evolution of the site. Since 1999, over 93% of implanted species have disappeared, but the diversity of helophytes seems to have stabilized since 2008 with only 12.5% of variation for the number of such species. On the basin of Django, the number of species remained stable since 2008. But, the decrease of almost 50% of macrophytes recovery rate, reduced to Phalaris arundinacea, Typha latifolia and Schoenoplectus tabernamontanii, emphasizes a change in entering water flow rates. The dominant flora found in infiltration basins evolves over time with euryece species of pseudo-metallophytes macrophyte group (Typha latifolia, Phragmites australis, Iris pseudacorus, Phalaris arundincea), ruderal species typical of dry grasslands (Asteraceae, Poaceae, Brassicaceae) and wetland species (Rumex sp.). But before all, these plants are mainly divided/distributed in theses basins according to hydrological factors, i.e. water availability among others.

KEYWORDS

Infiltration basin, temporal evolution, hydrological factors, helophytes, dominant vegetation

INTRODUCTION

Rainwaters, which run off over urban surfaces, are collected in basins and loaded with various pollutants such as pesticides and heavy metals [START_REF] Moura | Méthode d'évaluation des performances des systèmes d'infiltration des eaux de ruissellement en milieu urbain[END_REF][START_REF] Winiarski | The impact of stormwater on a soil profile in an infiltration basin[END_REF]Deschene, 2004). Furthermore, these waters are characterized by a high content of suspended solids and organic matter that will form a sediment layer at the soil surface (Deschene, 2004;[START_REF] Coustumer | Colmatage et rétention des éléments traces métalliques dans les systèmes d'infiltration des eaux pluviales[END_REF][START_REF] Gonzalez-Merchan | Amélioration des connaissances sur le colmatage des systèmes d'infiltration d'eaux pluviales[END_REF]. These basins and their sediments are colonized by vegetation, either artificial or natural [START_REF] Saulais | Colonisation végétale des bassins d'infiltration et de rétention[END_REF][START_REF] Bedell | Influences on the establishment and dominance of vegetation in stormwater infiltration basins[END_REF]. Indeed, the type of vegetation that is taking place in such environment depends on biotic and abiotic conditions as: access to light, oxygenation of the environment, availability of water and nutrients [START_REF] Saulais | Colonisation végétale des bassins d'infiltration et de rétention[END_REF]). In the past, studies on the infiltration basins have shown that the spatial distribution of the vegetation was mainly impacted by three factors: the water content of sediments [START_REF] Bedell | Influences on the establishment and dominance of vegetation in stormwater infiltration basins[END_REF], the sediment thickness [START_REF] Bedell | Influences on the establishment and dominance of vegetation in stormwater infiltration basins[END_REF] and pollutant contents [START_REF] Desjardins | Distribution patterns of spontaneous vegetation and pollution at a former decantation basin in southern Québec, Canada[END_REF]. As part of this work, we identified over a decade the evolution of species richness of flora and its abundance on two types of vegetated basins: one where the vegetation was originally planted and grown (maintenance), the other naturally colonized with no maintenance. From these observations, we tried to establish the temporal and spatial evolution of the vegetation at the scale of the basin and to identify the factors that might explain this trend.

MATERIALS AND METHODS

Site location

Two infiltration basins ("Minerve" and "Django") in the Lyon area were studied. These basins has been monitored over the last decade in the context of long term survey developed by the OTHU (Field Observatory for Urban water Management-www.othu.org).The "Minerve" basin collects rainwater from a catchment area with tertiary activity. Its vegetation is artificial and had been implemented at the beginning of the activity of the basin. The mowing and clippings are regular. Conversely, the basin "Django" presents a spontaneous and natural vegetation, without any maintenance. It receives stormwater from an industrial area.

Sampling of flora and sediments

The physicochemical measurement and floristic surveys were conducted from March to May 2015. Thus, 12 quadrats (Django) and 10 quadrats (Minerve), each of 1m² in area, and were positioned on the sites according to vegetation patches observed. For each quadrat, the number of species and their percentage coverage were noted. An abundance index was calculated for each species using the method of Braun-Blanquet (Wikum and Shanholtzer, 1978). Data in 2015 were compared with the data obtained earlier in 1999 to Minerve and again in 2008 for the two basins [START_REF] Saulais | Colonisation végétale des bassins d'infiltration et de rétention[END_REF]. The sediment of each quadrat for each basin was sampled to measure water content, the height of sediment layer, pH, organic matter and finally heavy metal content [START_REF] Roullier | Chapitre 31, Méthode d'analyses des sols, dans « Constituants et propriétés du sol[END_REF].

Statistical results

The processing of data are carried out using R© software. The normality test and Shapiro Wilk test revealed a part of the measured parameters could not be considered as normal distributions. As a result, the non-parametric Kruskal Wallis had to be considered for the tests of hypotheses. The significance of the results is assessed considering a threshold 0.05 for the p-value.

RESULTS

Sediment

High contents of organic matter (18% to 39%) and heavy metal (Zn Django: 1350 mg/mg MS to 2770 mg/mg MS) are observed in samples. In Minerve, we observed a gradient from upstream (near the water entry) to downstream (farthest from this water entry) for parameters like organic, height of layer and trace element contents. In Django, the distribution is less clear than for Minerve, but we also observed zones with contrasting contents in organic matter and heavy metal. The highest contents were observed as a function of water localization at the surface and/or height of sediment deposit.

Abundance and species richness

For Minerve basin, the floristic surveys from 1999 and 2008 have shown a drastic decline in species richness (Saulais, 2011) (Figure 1a, 1b). Most plants covering the basin were macrophytes species. In 2015, the basin includes 18 species (Figure 1c). Wetland zones are dominated by Acorus calamus, Typha latifolia, Iris speudacorus and Phragmites australis associated with other rare species helophytes: Schoenoplectus tabernaemontani, Solanum dulcamara and Phalaris arundinacea. However, dry grasslands appeared at the western end of the basin (farthest from water entry) and are covered by Asteraceae (including mainly Solidago canadensis), Geranium dissectum, Cardamine hirsuta, Rumex obtusifolius and Poa sp. This grassland area allows the development of a covered dense Asteraceae (Figure 1).

For Django, after the basin rehabilitation in 2002, no vegetation was sown or planted (Figure 2a). Six years later ( 2008), the inventories have highlighted a separate organization of vegetation including 58 species [START_REF] Saulais | Colonisation végétale des bassins d'infiltration et de rétention[END_REF] (Figure 2b). Distinct vegetation groups emerge and reflect mainly 3 types of contrasting ecological environments: disturbed environments, dry and stony zones (Senecio inaequidens, Capsella bursa pastoris, Cardamine hirsuta), wetland zones with vegetation subjected to cycles of flooding and drying up (Phalaris arundinaceae, Polygonum mite, Galium aparine) and areas regularly submerged by water (Typha latifolia, Eleocharis palustris) [START_REF] Saulais | Colonisation végétale des bassins d'infiltration et de rétention[END_REF]. In 2015, the floristic survey shows a relative stability in species richness (57 species with most of them being unchanged) (Figure 2c). The three distinct ecological environments are always present but with different recovery rates. The wetland zone is dominated by Phalaris arundinacea and Typha latifolia, associated with other rare helophytes: Cyperaceae and Solanum dulcamara. An "intermediate" zone sporadically flooded is covered mostly with Rumex sp. This Polygonaceae belt defines the boundary between the wetland and the rest of the basin flooded only during intense rainfall events. Many Asteraceae, Brassicaceae, Poaceae and Lamiaceae characterize the driest areas (Figure 2).

DISCUSSION

The regression of helophytes communities on Minerve for the benefit of ruderals pioneer species marks the evolution of the site. The appearance of a separate grassland area from the wetland is the result of such regression. Since 1999, over 93% of implanted species have disappeared. In contrast, the diversity of helophytes seems to have stabilized since 2008: in 7 years, the number of species varied only by 12.5%.

On the basin of Django, the number of species remained stable since 2008. In addition, the decrease of almost 50% of macrophytes recovery rate, reduced to Phalaris arundinacea, Typha latifolia and Schoenoplectus tabernamontanii, emphasizes a change in entering water flow rates. The installation and maintenance of the species Acorus Calamus, Iris pseudacorus, Phalaris arundinacea, Phragmites australis and Typha latifolia is not surprising because these helophytes have a very broad ecological valence (IBMR Standard AFNOR, 2004). They support also any type of environment: oligotrophic to eutrophic and those contaminated with organic and inorganic pollutants (IBMR Standard AFNOR, 2004). They can be well described as pseudo-metallophytes (Remon, 2009). The pathways of the water at surface seems to have changed over time in the two infiltration basins. In Minerve, water and sediments accumulate preferentially at the base of the bridge, but the eastern end of the basin is always ponded with water. Due to such water ponding, the helophytes (Iris pseudacorus, Typha latifolia) are able to remain in this zone. A resurgence of Phalaris arundinaceae (reduced to a few individuals in 2008) was also observed, which may reveal a natural recolonization of the basin (these species tolerate both alkaline and acid soils). In contrast, the presence of a "dry" area and a wetland highlights a chronic drought episode in Minerve basin. Pioneering grassland species such as Asteraceae, have undoubtedly benefited from the decline of Iridaceae, following the reduction of soil water content. Over time, a steady decline in the number of Iris pseudacorus could be observed. In Django, vegetation also seems to follow stormwater pathways. The helophyte vegetation seems to have moved and been replaced by grassland vegetation. In 2008, 20% of the surface of the basin was naked [START_REF] Saulais | Colonisation végétale des bassins d'infiltration et de rétention[END_REF], whereas today it is only 6%. The increase on sediment deposit at the surface as well as most water availability in several zones can explain the increase of the surface vegetalized. The pioneer vegetation of meadows, typical of dry and rocky areas, may have colonized the area.

CONCLUSION

The dominant flora found in infiltration basins evolves over time. It brings together euryece species of pseudo-metallophytes macrophyte group (Typha latifolia, Phragmites australis, Iris pseudacorus, Phalaris arundincea), ruderal species typical of dry grasslands (Asteraceae, Poaceae, Brassicaceae) and wetland species (Rumex sp.). They can be classified for most of them as helophytes with a very broad ecological valence. Indeed, they support very well any type of environment: oligotrophic to eutrophic and environments with highs loads of organic and inorganic pollutants, which is one of the reason why these are referred to as pseudo-metallophytes. But before all, these plants are mainly divided/distributed in theses basins according to hydrological factors, i.e. water availability among others.

ACKNOWLEDGEMENT

The authors also wish to acknowledge OTHU (www.othu.org) for monitored database for which the Greater Lyon Metropolis and and the Water Agency RMC provides technical and financial supports.


Stefano Curiotto
email: curiotto@cinam.univ-mrs.fr

Pierre Müller
email: muller@cinam.univ-mrs.fr

Fabien Cheynis
email: cheynis@cinam.univ-mrs.fr

Frédéric Leroy
email: leroy@cinam.univ-mrs.fr

Kinetic Monte Carlo simulations of the diusion and shape evolution of single-layer clusters on a hexagonal lattice with and without external force

    Kinetic Monte Carlo simulations of the diffusion and shape evolution of single-layer clusters on a hexagonal lattice with and without external force

Introduction

Surface diusion plays a key role in catalysis, adsorption, growth, reactivity, as it aects the stability and the morphology of material surfaces by means of step rearrangement [START_REF] Burton | Anisotropic diusion in stress elds[END_REF] and cluster diusion [START_REF] Voter | Classically exact overlayer dynamics: Diusion of rhodium clusters on Rh(100)[END_REF]. Surface diusion is the most common surface transport process, so that recently it has been proposed to use it to control the motion of nano-objects on a surface [START_REF] Curiotto | 2D Manipulation of Nanoobjects by Perpendicular Electric Fields: Implications for Nanofabrication[END_REF]. This would open new perspectives in nanofabrication. In this context, understanding cluster diusion on various surfaces is becoming a new challenge. Clusters predominantly migrate by the individual motion of atoms instead of collective motion. In absence of any external force, the motion of the clusters is random (Brownian), whereas in presence of an external force, as for instance electromigration [47], a bias is added to the random motion and thus enables to control the cluster motion. If cluster diusion is now well documented by experiments and analytical approach (at least for large islands) [START_REF] Khare | Brownian motion and shape uctuations of single-layer adatom and vacancy clusters on surfaces: Theory and simulations[END_REF], Kinetic Monte Carlo (KMC) simulations have brought new perspectives especially for studying small clusters for which atomic congurations play an important role [START_REF] Lai | Diverse nanoscale cluster dynamics: Diusion of 2D epitaxial clusters[END_REF][START_REF] Lai | Diusion of twodimensional epitaxial clusters on metal (100) surfaces: Facile versus nucleation-mediated behavior and their merging for larger sizes[END_REF]. However, most of the studies have been performed on square lattices. In this article we use KMC simulations to study the diusion of singlelayer clusters (2D voids and 2D islands) on a hexagonal (or triangular) lattice in presence or in absence of an external force like electromigration. We show that most of the results obtained on square lattices are still valid on hexagonal lattices, even if subtle changes are evidenced.

For the largest clusters, simple scaling-laws expressed in terms of size-dependent diusion coecient do not depend on the lattice symmetry. For small clusters, the shapes and shape uctuations play a major role, leading to the appearance of nite-sizes eects that may dier on square or hexagonal lattices. Our goal in this article is to give a clear description of the mechanisms of cluster diusion on hexagonal lattices with or without an external force. For this purpose, we start by describing some concepts of surface diusion from a microscopic point of view, taking into consideration single atomic jumps with a special focus on how atomic jumps lead to cluster diusion. We pay special attention to small clusters for which simple scaling laws are not valid, leading to specic behaviors for welldened cluster sizes. The step-by-step study of the cluster motion in terms of single atomic jumps aims at providing a detailed description of the cluster diusion properties and shape evolution. Section 2 details the model used in the simulation. Section 3 addresses the diusion of clusters without the eect of external forces. We pay particular attention to small clusters, where discrete eects modify the behaviors predicted by continuous models. Section 4 shows how the clusters displace under a force, their velocity and shape for dierent temperatures and force values. In section 5 we change the symmetry of the lattice from 6-fold to 3-fold and we simulate the shapes of clusters with and without the eect of an external force. The most relevant ndings are then summarized in the conclusions.

Kinetic Monte Carlo model

In our model, we consider a 2D hexagonal (or triangular) lattice of positions that can be either occupied by atoms or empty (see gure 1a). The simulation box has at least 10000 lattice sites and periodic boundary conditions. Atoms at steps A and B, like the red and green atoms, have the same number of neighbors but they have a dierent local arrangement (the red atom has only one neighbor in the bottom layer at the exterior of the step, while the green atom has two neighbors in the bottom layer, see the short black lines). The atoms under the top layer that are neighbors of the red and the green atoms are indexed nsl1, nsl2 and nsl3 (see appendix A).

Atoms can jump from an occupied to an empty neighbor position and cannot jump on top of already occupied positions, therefore atomic jumps are not allowed on top of islands and outside voids (innite Ehrlich-Schwoebel barrier).

Figure 2: Schematics of the energy barrier that an isolated adatom has to overcome to jump from a site i to a nearby position i + 1 or i -1. With a force (red horizontal arrow, right panel), the jump energy is no more symmetric in the two directions: it is lower for a jump towards the left (yellow curved arrow), in the force direction, while it is higher for a jump towards the right (violet curved arrow). A jump in the force direction is thus favored.

The model is based on a previous version [START_REF] Curiotto | Shape changes of two-dimensional atomic islands and vacancy clusters diusing on epitaxial (111) interfaces under the impact of an external force[END_REF]. The jump probability depends only on the energy of the departing site (Arrhenius dynamics), i.e. on the binding energy of the jumping adatom (E bin , see appendix A) and on the external force acting on it (see the schematics of gure 2). The jump rates are proportional to:

exp(- E bin + E cg kT )
k is the Boltzmann constant and T is the temperature.

E cg is the energy change due to the external force. It is positive when the jump direction is against the force and negative when the jump is helped by the force. We have used the expression of E cg given in [START_REF] Pierre-Louis | Electromigration of single-layer clusters[END_REF], but modied for a hexagonal lattice:

E cg = |F | • cos( π 3 • b -δ) • a 2
which takes into account the angle between the force and the atomic jumps. |F | is the external force acting on an atomic jump; a is the lattice separation between two neighbor positions, we take it equal to 1 in all directions; b is an integer that denes the arrival site in anticlockwise order, with values between 1 and 6 ; δ is the angle between the force and the x axis. Atomic jumps are selected according to a rejection-free Kinetic Monte Carlo algorithm. The time needed for an atomic jump is equal to the inverse of the sum of all jump rates [START_REF] Voter | Radiation Eects in Solids[END_REF]. The results of section 5 are obtained with an improved version of the model, modied to represent also the (111) surface plane of cubic structures, where atomic steps can form either 100 or 111 microfacets (A and B steps in the following) with the underlying plane, as shown in gure 1b. B steps are more stable than A steps [START_REF] Michely | Visualizing the electron scattering force in nanostructures[END_REF][START_REF] Rousset | Selfordering of Au(111) vicinal surfaces and application to nanostructure organized growth[END_REF]. It is not possible to reproduce this energy dierence taking into account only the number of bonds of surface atoms with in-plane neighbors. In order to consider the step energy dierence in our simple lattice model, we have used a binding energy that depends on the bonds formed by surface atoms with the subsurface atomic layer, that is considered fully occupied, according to a procedure detailed in appendix A. This model is general for 2D triangular lattices and illustrates qualitatively dierent scenarios of homo-epitaxial cluster diusion biased by a force. It was originally developed to understand the diusion of 2D Si islands on Si(111) (see some experimental results in [START_REF] Leroy | Electric forces on a conned advacancy island[END_REF]), but the energies and temperatures used are not characteristic of a particular system. The model could be adapted to dierent systems by changing the binding energies (see appendix A).

Cluster diusion

The model can be used to simulate random diusion when the force |F | is set to 0. The simulations of diffusion of single adatoms and vacancies, available in the supplementary material, allow to validate the model. In the following we focus on the diusion of clusters. a consequence of atom diusion, the edges of islands and voids are modied and the cluster center of mass displaces. Three dierent mechanisms of cluster motion are usually considered: atomic diusion at the cluster periphery, diusion on the terrace outside and inside clusters, and uncorrelated detachment/attachment of atoms (evaporation-condensation). Simple scaling arguments lead to the conclusion that diusion should depend on the cluster radius R: D C ∝ R -α , with α=1, 2 or 3 for kinetics limited by attachment/detachment, terrace diusion and periphery diusion respectively [START_REF] Morgenstern | Brownian motion of vacancy islands on Ag(111)[END_REF][START_REF] Van Siclen | Single jump mechanisms for large cluster diusion on metal surfaces[END_REF]. Khare and Einstein [START_REF] Khare | Brownian motion and shape uctuations of single-layer adatom and vacancy clusters on surfaces: Theory and simulations[END_REF][START_REF] Khare | Diusion of Monolayer Adatom and Vacancy Clusters: Langevin Analysis and Monte Carlo Simulations of their Brownian Motion[END_REF] used an analytical continuous description of the cluster boundary to describe the diusion of clusters, conrming the above nding. Several authors, the rst of them being Voter [START_REF] Voter | Classically exact overlayer dynamics: Diusion of rhodium clusters on Rh(100)[END_REF], also considered the discrete nature of clusters, to take into account the eect of kinks and facets, and performed Monte Carlo simulations to investigate cluster diusion. Most of them focused mainly on the periphery diusion and evaporation-condensation mechanisms [START_REF] Lai | Diverse nanoscale cluster dynamics: Diusion of 2D epitaxial clusters[END_REF][START_REF] Lai | Diusion of twodimensional epitaxial clusters on metal (100) surfaces: Facile versus nucleation-mediated behavior and their merging for larger sizes[END_REF]2228], nding some discrepancies from the power laws detailed above. Cluster diusion mediated by atoms diusing on the terrace has been less studied, but Scanning Tunneling Microscopy [START_REF] Morgenstern | Brownian motion of vacancy islands on Ag(111)[END_REF] and KMC simulations on a rectangular lattice to describe the Ag(110) surface [START_REF] Morgenstern | Motion of vacancy islands on an anisotropic surface: Theory and kinetic Monte Carlo simulations[END_REF] conrm the power law and α=2.

Large clusters

In order to get further insights into the diusion of clusters on a surface with atoms organized in a compact hexagonal pattern, we have simulated the evolution of clusters when atoms detach from the perimeter, diuse on the terrace and re-attach to the cluster. Notice that periphery diusion is not impeded. We have veried that the cluster average square displacement is proportional to the time, so that we can dene the cluster diusion coecient

D C = <r 2 cluster > 4t
. The number of atomic jumps of each simulation ranges between 10 8 (cluster with R=5) and 10 10 (R=45). The convergence of the D C measurement is veried using the procedure detailed in the supplementary material, where control charts are shown.

The black squares in gure 3a show the diusion coecient of large islands for dierent radii in log-log scale (R larger than 8 units, more than 200 atoms in the clusters). A linear t through the points gives the value of α = 1.9 ± 0.1, close to 2. Diusion of adatoms on the terrace controls the cluster diusion. Islands are surrounded by a large terrace (the box size is much larger than that of the island) where adatoms can diuse for a long time before re-attaching, therefore the kinetics for island displacement in that case are limited mainly by terrace diusion.

For voids, we nd that α is dierent for small (α = 1.54 ± 0.04) and large (α = 1.92 ± 0.02) voids (empty circles in gure 3a). As suggested by Sholl et al. [START_REF] Sholl | Diusion of clusters of atoms and vacancies on surfaces and the dynamics of diusion-driven coarsening[END_REF], in small voids the diusion time of an adatom between a detachment and an attachment event is small. The kinetics for the displacement of the cluster depends thus also on the time needed to detach an atom and not only on the time of adatom diusion inside the void. The detachment eect contributes to α that is thus between 1 and 2. At lower temperatures, the exponent α further decreases because the detachment rates decrease more than the adatom diffusion rate. For voids we nd α = 1.54, 1.45 and 1.3 at kT=0.4, 0.3 and 0.25 respectively.

Voids displace slower than islands, but the diusion dierence decreases when the cluster size increases. The reason of this is that the step edge of islands is convex and therefore presents more atoms in kink positions, easy to detach from the perimeter, than voids, where the step edge is concave. In other words, on average, atoms at islands edges are less bound than atoms at void edges. This dierence between islands and voids is exacerbated for small clusters. The detachment time of atoms from the periphery of small voids is thus expected to be high, and further accentuates the importance of the detachment kinetics found for small voids.

We now test the expression for the cluster diusion coecient given by Van Siclen [START_REF] Van Siclen | Single jump mechanisms for large cluster diusion on metal surfaces[END_REF], developed for terrace diffusion or correlated attachment-detachment mechanisms (i.e. an atom detaching from the cluster edge has higher probability to re-attach in a position close to the detachment site):

D C = < r 2 cluster > 4t = D ad • c • Ω 2 π • R 2 (1)
Ω is the unit cell area and c the adatom concentration. In the simulations it is possible to measure the average square displacement and also the number of adatoms around an island or inside a void. We can thus measure an adatom concentration. The measured cluster diusivity can be compared with that predicted by the above equation using the measured concentration. Table 1 compares some values of diusivity for islands and voids. The above equation is based on a continuous model, while in our simulations the atomic jumps are discrete and the displacement of a cluster could depend on facet nucleation. However the agreement is reasonable for both voids and islands. The high error on the calculated D C is due to the error on the measured adatom concentration. Our equations are valid for diusion only outside islands and inside voids. The expression for the diusion coecient given by Khare and Einstein [START_REF] Khare | Brownian motion and shape uctuations of single-layer adatom and vacancy clusters on surfaces: Theory and simulations[END_REF] corresponds to that of Van Siclen multiplied by a factor 2 because they allow diusion both outside and inside the clusters (for instance considering also diusion of adatoms on top of the islands). We verify if, in analogy with the diusion of a single adatom, we can write

D C = D C0 •exp -E cluster kT
, where D C0 is a prefactor and E cluster an eective cluster diusion energy. The black squares in gure 3b show an example of cluster diusion coecient in natural logarithm scale as a function of 1/(kT ) for an island, while the empty circles represent the diusivity of a void with 10 atomic units radius. For each set, the points are well aligned, and, using the expression above, a linear t gives an eective energy of E cluster =6.06 ± 0.05 E b for islands and 6.22 ± 0.05 for voids. E cluster does not correspond to the adatom diusion energy (E ad = 3E b ). We have changed the adatom diusion energy independently from other binding energies to check if E cluster changes. As shown by the empty pentagons and stars in gure 3b, decreasing (resp. increasing) E ad , the cluster diusivity D C increases (resp. decreases) only slightly, and the slope of the linear t in gure 3b, i.e. E cluster , does not change. While cluster diusion as a function of size (D C ∝ R -2 ) suggests the importance of the adatom diusion on the terrace, the adatom diusion energy E ad does not play a role in the cluster diusion energy. These ndings are consistent with those of Heinonen et al. [START_REF] Heinonen | Ala-Nissila, Island diusion on metal fcc(100) surfaces[END_REF], who simulated cluster diusion on (100) lattices but without extensive discussions.

To understand E cluster we have to consider in equation 1 also the adatom concentration and not only D ad . The concentration c depends from the balance between atoms that detach from edges and adatoms diusing and jumping towards an edge and attaching to it [30]. In our simulations a cluster is at equilibrium with an adatom gas surrounding it and there is no other structure that can capture the atoms detached from the cluster. Therefore, at long times, the number of atoms detaching from the cluster equals that of atoms attaching to it. On average, the time necessary to detach an atom is equal to the time needed to have an attachment to the cluster boundary. Adatoms on the terrace attach to a cluster with a probability proportional to the adatom concentration, to the adatom diusion coecient and to the length of the cluster perimeter (that can be considered as a cross section). The attachment time, inversely proportional to the attachment probability is then

t attachment ∝ 1 c • D ad • 2πR
The probability to detach an atom from the cluster boundary depends on the number of boundary atoms N boundary and on their binding energy (i.e. on their neighbors). t detachment is inversely proportional to this probability:

t detachment ∝ 1 5 j=1 N j • exp -Ej kT
where j is the number of in-plane neighbors, E j is the binding energy of an atom with j in-plane neighbors, N j is the number of atoms at the edge of the cluster with j neighbors (see appendix A and gure A.14) and

5 j=1 N j = N boundary . Writing 5 j=1 N j • exp -Ej kT N boundary =< exp -E j kT > then t detachment ∝ 1/(N boundary < exp( -E j kT ) >) (2)
This expression cannot be directly handled because the average of the exponentials depends on the N j , that depend on the size of the cluster and on the temperature. We thus re-write the expression using an eective binding energy E bin for all the boundary atoms, that can be found considering the energy needed to remove an entire row of a cluster edge. If the average number of atoms of an edge is N edge , and they are removed one after the other, N edge -1 are kinks, while the last one is more weakly bound (atom adsorbed at a step). Therefore

< E bin >= E kink (N edge -1) + E ad-step N edge
We thus write

t detachment ∝ 1 N boundary exp( -E kink (N edge -1)+E ad-step N edge •kT ) With t attachment = t detachment and considering that N boundary ∝ 2πR, we nd c ∝ exp -E kink kT • exp E ad kT • exp E kink -E ad-step N edge • kT cluster type kT simulated D C calculated D C island 0.45 8.6•10 -9 ± 0.5 • 10 -9 6.6•10 -9 ± 0.7 • 10 -9 island 0.5 3.4•10 -8 ± 0.3 • 10 -9 2.6•10 -8 ± 0.4 • 10 -8 void 0.45 1.5•10 -9 ± 0.2 • 10 -9 1.8•10 -9 ± 0.4 • 10 -9 void 0.5 5.9•10 -9 ± 0.2 • 10 -9 7.1•10 -9 ± 0.4 • 10 -9
Table 1: Comparison of measured(simulated) and calculated (equation 1, Van Siclen expression [START_REF] Van Siclen | Single jump mechanisms for large cluster diusion on metal surfaces[END_REF]) cluster diusion coecients, for islands (R=10) and voids (R=20). The calculated D C has an error due to the error on the simulated concentration.

Notice that this expression includes, in the denominator of the exponential, the dependence of the concentration on the cluster size (N edge ∝ R), as predicted by the Gibbs-Thomson equation. Substituting the expression above in equation 1, we obtain:

D C ∝ exp -E kink kT • exp E kink -E ad-step N edge • kT (3)
The argument of the second exponential term in the above expression is negligible with respect to the argument of the rst one and thus E cluster E kink = 6E b in our simulations. However, changing the cluster radius, we have found that E cluster depends on it, though very slightly, and increases when R increases (For 2D islands E cluster = 5.92, 6.06 and 6.10 E b for R=8, 10 and 12). This is qualitatively in agreement with the eect of the second exponential term in equation 3. Another approach to discuss E cluster that does not imply an eective energy is detailed in appendix B. Finite-size eects have been experimentally observed in nanoclusters of some systems. For instance, the catalytic properties of Pd nanoparticles show non-monotonic variations as a function of their size [START_REF] Sitja | Transition from Molecule to Solid State: Reactivity of Supported Metal Clusters[END_REF][START_REF] Sitja | Particle size eect on the Langmuir-Hinshelwood barrier for CO oxidation on regular arrays of Pd clusters supported on ultrathin alumina lms[END_REF]. Also the diusion of small clusters can depend on their size in a non-trivial way, especially at low temperatures. As shown in the previous section, the diusion coecient decreases when the cluster size increases, and, for islands, it is proportional to R -2 . However, for small islands, this general trend changes. Figure 4 shows that D C decreases with the island size N, i.e. the starting number of atoms in the cluster.

Finite-size eects: 2D islands

For large clusters, N is ∝ R 2 , therefore D C ∝ R -2 ∝ N -1 ;
for small clusters R is not signicant, however we verify if the D C proportionality still holds. The inset of gure 4, shows that log D C is proportional to -1.2(±0.3) log N only for islands larger than about N=29. For smaller sizes D C decreases faster than N -1 and ceases to follow a size dependent power law. Figure 4 also shows that D C does not decrease steadily when N increases, but is punctuated by oscillations.

On square lattices, dierent authors have found that cluster diusion does not follow a monotonic descending trend with the cluster size [START_REF] Lai | Diverse nanoscale cluster dynamics: Diusion of 2D epitaxial clusters[END_REF][START_REF] Heinonen | Ala-Nissila, Island diusion on metal fcc(100) surfaces[END_REF]3335]. This is due to compact or "perfect size" clusters that are slower than others because their boundary is only constituted by strongly bound atoms (with few kinks and no adatoms at steps). The reasons of the oscillations are explained for (100) surfaces by Lai et al. [START_REF] Lai | Diverse nanoscale cluster dynamics: Diusion of 2D epitaxial clusters[END_REF]: clusters with perfect size, without adatoms at steps, diuse slowly because the formation of a new edge allowing the cluster displacement requires a nucleation energy equivalent to the detachment of two kink atoms.

On a hexagonal lattice, clusters with symmetric perfect size are for instance those with N=7, N=19 and N=37, as shown in gure 5. Those congurations have the minimum number of border atoms with the maximum possible inplane neighbors. They have 6 kinks, no adatoms at steps, and thus minimize the cluster step energy. However, other sizes (for instance N=10, 12, 13, 14, 16 and for all clusters with N≥18) also exhibit shapes with 6 kinks and no adatoms at steps (see two examples in gures 5 for N=12 and N=13). To understand the D C oscillations, we consider that a signicant displacement of a cluster requires the detachment and re-attachment of atoms, formation of a new edge with at least two atoms. For instance, starting from a cluster with N=19 shown in gure 5, after the removal of a rst kink (blue in gure 5), two other kinks are created (green in gure 5), and all the boundary atoms still have at least 3 in-plane neighbors. Therefore the formation of a new edge requires at least detachment of two kinks. For N=20, the low energy conguration has 4 kinks that are coupled: removing one of them leaves an adatom at C . For N=21, the compact shape has only 2 coupled kinks, less than N=20, therefore the probability to remove one of the coupled kinks (and to leave one adatom at step with 2 in-plane neighbors) is lower than for N=20, so D 20 C > D 21 C . The compact congurations with N=12 and N=13 have both 6 coupled kinks. However for N=13 two of the kinks are special because removing one of them (for instance the blue one in gure 5) leaves two adatoms at step (green in gure 5) instead of one. Therefore D 13 C > D 12 C . Only two coupled kinks are present for N=14, so D 14 C < D 13 C . Similar reasoning can explain all the other oscillations of D C in the red boxes of gure 4. No compact shapes with only 6 kinks can be found for clusters with N=15 and N=17, they have either one adatom at step on the boundary or a shape with 7 kinks. These two clusters should emit the adatom at step very easily to show a compact shape, and this atom should help the nucleation of new steps and thus should make the island displace faster than the island with size N-1. However this is not observed in the value of D C , i.e. no signicant oscillations of D C are observed between 14 and 17. The increase of D C at N=18 is due to the 6 coupled kinks of the cluster compact shape.

It should also be considered that low energy congurations can have, for dierent sizes N, dierent number of atoms constituting the steps. A conguration with more step atoms has more probability to emit a boundary atom, and thus to make the cluster move, than a conguration with less step atoms. This eect could add to the "coupled kinks" eect explained in the previous paragraph. Another way to discuss the oscillations is to consider the number of possible congurations that can be taken by a diusing cluster. This number increases with the cluster size, but, for xed N, some congurations have lower energy and are thus more probable than others. For each cluster with size N, we have run simulations to allow for the rearrangement of a cluster shape. During the simulation, at xed intervals, the cluster conguration was compared with those already found. At the end of the simulation, we have counted the number of times all the dierent cluster congurations appears. Figures 6a and6b show two histograms of the occurrence of the dierent congurations for the cluster with size N=12 and N=13 (red and dark blue bars respectively). The cluster with N=12 tends to have 2 dierent congurations, while others are much less probable. The island with N=13 can have many dierent congurations (≈ 30), with similar probability, while other congurations are less probable. It can thus be assumed that, for the cluster with size N=12, atoms detach but often re-attach back at the same position, leading to few probable cluster congurations and delaying cluster diusion. The island with N=13 can have dierent congurations with similar probability and should thus diuse faster than that with N=12. The insets in gure 6a and 6b show another example of congurations found for clusters with N=19 and N=20. As expected, the cluster with N=19 has very often the conguration corresponding to the symmetric perfect size (see gure 5). Figures 6c shows the number of most probable congurations as a function of the cluster size, together with the diusivity corrected of the monotonic descending trend. A conguration γ is considered probable if it is found more than 100 times with respect to a low probable conguration and if there are no other congurations that occur 5 times more than γ. The clusters with few probable congurations, i.e. those that tend to show very often the same shape, have lower diusivity. An exception, addressed above considering the number of step atoms in the compact shape, is for N=17 and 18. For clusters larger than N=35, the oscillations are within the error of the measured diusion coecient. The number of congurations is also an argument evoked on square lattices [START_REF] Lai | Diverse nanoscale cluster dynamics: Diusion of 2D epitaxial clusters[END_REF]. The main dierence is that, as hexagonal islands have higher symmetry than square ones, the oscillations in gure 4 have only weak amplitude and they disappear for N>35, that is a size much smaller than that observed for square lattice (N=200/300). The perfect size eect, as also highlighted by Heinonen et al. [START_REF] Heinonen | Ala-Nissila, Island diusion on metal fcc(100) surfaces[END_REF], is only observed at low temperatures, where the dierence of Boltzmann factor among the dierent atomic jumps is large.

Finite-size eects: 2D voids

Like islands, also voids show a nite-size dependence of the diusion coecient from the size. Figure 7 shows D C as a function of the size of voids. Clusters with small size often splits during the long simulations necessary to average the cluster squared displacement, and are thus not considered. The general trend is an increase of D C with the size up to about N=50, where a plateau is found. Then D C decreases at rst slowly and, for N>300, steadily, with D C ∝ N -0.65±0.05 . Notice that N -0.65 corresponds to R -1.3 found in the cluster-diusion section. Most of the data, within the simulation error, follow a continuous curve. However, as in the case of islands, there are some clear deviations of D C from the general trend, particularly for N=7, 10 and 19. The increase of D C with the size for small N is surprising, was not observed for islands, was not investigated by other authors and thus deserves a discussion. As in the case of single adatoms,

D C = 1 4 Γ C • λ 2 C
, where Γ C is the cluster jump frequency, that depends on the sum of the detachment rates of boundary atoms and on the diusion time of adatoms inside the voids; λ C is the jump length, that decreases slowly when the void size increases. Increasing the void size, the number of boundary atoms that can move increases and their jump rates also increase. For instance, the motion of an isolated vacancy (D C marked with a red circle in gure 7a) depends on the jump rates of the 6 atoms around it, ∝ exp -8E b kT . The "perfect" shape of a void with N=7, that should be very stable, has 6 boundary atoms with jump rate exp -8E b kT and 6 other border atoms with jump rates 2 exp -7E b kT . It thus displaces faster than the isolated vacancy. The sum of the detachment rates of boundary atoms in clusters with nonperfect shape is even higher than that of the void with N=7 and thus they move faster. For instance the void with N=6 in gure 7b has a boundary atom with high detachment rate (3 exp -6E b kT ), and the void with N=8 has seven boundary atoms (instead of six for the cluster N=7) with detachment rate 2 exp -7E b kT . For small voids the diusion time inside the voids is not signicant and thus Γ C increases with the size. Increasing the void size, the cluster jump length decreases and the diusion time of the adatoms in- kT , while atoms in a cluster are more strongly bound and have smaller jump rates, ∝ exp -4E b kT or smaller). Furthermore, the diffusion time outside islands before attachment is high for small islands and decreases when the island size increases (increased cross section and thus increased probability of adatom attachment). Therefore Γ C (and thus D C ) always decreases when the island size increases.

The behavior described above concerning the diusivity of voids with dierent sizes is not typical of hexagonal lattices: In preliminary simulations on a square lattice we have observed the same features (see the inset of gure 7a). Applying a force, the diusion of single atoms and vacancies is biased in the force direction. The grey area in gure 8a shows the probability distribution of the nal position of a vacancy starting from the point marked by X and diusing (without force, E cg =0) for a xed time. The gure is based on 500 simulations, a darker color corresponds to higher probability. Applying a weak force on the atoms towards the left, the vacancy does not diuse anymore completely randomly around the starting position but displaces preferentially towards the right. This is shown by the purple probability distributions of the nal positions of a vacancy for 500 simulations performed with E cg =0.03. The average nal position of the vacancy is displaced towards the right, proportionally to the applied force. For strong forces the distribution is less dispersed, the vacancy nal position is well dened, with only a small deviation from the most probable position. Similar results are obtained for single atoms, the main dierence is that they displace in the same direction of the force and the jump energy is lower (3E b instead of 8E b for vacancies). Considering the time elapsed in the simulations, the vacancy velocity can be measured and is shown in gure 8b as a function of the force. When a force aects a jump, the jump frequencies are modied and the velocity v against the force direction is equal to the dierence between the jump frequencies along and against the force of the atoms surrounding the vacancy. Considering the 6 possible jumps of the vacancy:

Displacement under an external force 4.1. Single atoms/vacancies

v = a • (exp(- E vac -E cg kT ) + 1 2 exp(- E vac -0.5E cg kT ) - 1 2 exp(- E vac + 0.5E cg kT ) -exp(- E vac + E cg kT ) - 1 2 exp(- E vac + 0.5E cg kT ) + 1 2 exp(- E vac -0.5E cg kT )) = a • exp(- E vac kT ) • (exp( E cg kT ) -exp(- E cg kT ) + exp( 0.5E cg kT ) -exp(- 0.5E cg kT ))
The velocity as a function of the force found in the simulations agree perfectly with the equation above, as shown by the t of the dashed line with the black squares in gure 8b. For small forces (and thus small E cg ), exp x = 1 + x and therefore

v = 3a E cg kT • exp(- E vac kT ) = 3a 2 F 2kT • exp(- E vac kT )
because, as shown in gure2, the energy E cg is related to the force acting on the atoms as F • a 2 = E cg . As detailed in the supplementary material, where we study the diusion of single atoms and single vacancies,

D vac = 3 2 • a 2 • exp(-Evac kT ). Therefore, v = F • D vac kT
This expression corresponds to the Nernst-Einstein relation.

Clusters

Under the eect of a force, as a consequence of atomic displacement, clusters also move. This motion can be limited by diusion of atoms at the cluster periphery, by diusion on the terrace, and by attachment/detachment kinetics (this last mechanism becomes important if diusion is easy and fast with respect to the detachment). Within our model all these mechanisms can take place, and the balance among them depends on parameters like the Boltzmann factor and the size of the islands. As atomic jumps depend on the number of neighbors of the starting position and not on that of the nal state (we use KMC in Arrhenius dynamics [START_REF] Voter | Radiation Eects in Solids[END_REF] and not Metropolis Monte Carlo), periphery diusion is not particularly favored, contrary to other works where atoms are constrained to move only at Figure 9: Schematics of the motion of 2D islands and voids as a consequence of forced diusion. 2D voids always move opposite to the force acting on the adatoms, while 2D islands can move in the force direction when the motion is determined by periphery diusion.

Figure 10: a: Velocity of islands (black squares) and voids (red empty circles) displacing under the eect of a force, as a function of the radius. Ecg=0.001, kT=0.3. b: Velocity of an island (R=20) displacing in the direction opposite to the force, as a function of the force. The inset is a magnication of the main gure at low forces. the cluster periphery (as for instance done in [START_REF] Pierre-Louis | Electromigration of single-layer clusters[END_REF] in the section on periphery diusion).

While in our model the atoms diuse preferentially in the force direction, both islands and voids displace opposite to the force direction. As explained in [START_REF] Pierre-Louis | Electromigration of single-layer clusters[END_REF] and summarized in gure 9, this occurs because in our simulations periphery diusion does not play an important role. Clusters displace with constant velocity under the eect of a force. The velocity of islands and voids as a function of their size is shown in gure 10a. While the velocity of islands does not depend on their size, that of voids is proportional to the cluster size for small sizes, and becomes constant after a threshold. According to [START_REF] Pierre-Louis | Electromigration of single-layer clusters[END_REF], if the velocity of clusters moving under a force does not depend on their size, then the process limiting the displacement is terrace diusion, while if the velocity is directly proportional to the size, then attachment and detachment processes play the most important role. Our ndings for clusters subject to a force are thus consistent with those discussed in section 3, where we showed that for islands and voids the displacement is diusion limited, and for small voids also the detachment kinetics plays an important role. Furthermore, the velocity of voids increasing with the void size and reaching a constant value correctly reproduces the experimental observations of 2D voids displacing under electromigration on Si(111) [START_REF] Leroy | Electric forces on a conned advacancy island[END_REF].

As shown in gure 10b, for a xed size the cluster velocity increases linearly with the force for small forces and follows an exponential law for strong forces. Supposing that the equation derived for vacancies in section 4.1 (and also valid for adatoms) can be used for clusters, we can t the data to obtain E cluster (instead of E vac ). We nd E cluster = 6.19 ± 0.02E b (for kT = 0.3E b ) that is close to the diusion energy of a cluster obtained in section 3, without forces (E cluster = 6.06 ± 0.05E b ). The agreement is very good, the slight dierence is probably due to the clusters stability (they have dierent shapes for dierent forces, see section 5.2). For weaker forces, where the cluster shape changes marginally, a t with E cluster = 6.06 reproduces well the data (see the black line in the inset of gure 10b).

Cluster shapes in a lattice with 3-fold symmetry

The shape of clusters uctuates around an equilibrium shape that depends on the temperature, on the binding energy and on the external force applied. For hexagonal clusters with 6-fold symmetry, when k nsl = 0, the shape changes have been shown in [START_REF] Curiotto | Shape changes of two-dimensional atomic islands and vacancy clusters diusing on epitaxial (111) interfaces under the impact of an external force[END_REF]. Here, we discuss the shapes of clusters on a lattice with 3-fold symmetry. Figure 11 shows the shape of an island at dierent temperatures, for k nsl = 0.18 (a parameter dening the relative stability of A and B steps, see appendix A). For each temperature, the simulations start from a round cluster; then atomic jumps are allowed, so that the shape rearranges; when the shape does not signicantly evolve anymore, an average shape can be calculated. At low temperature corner and edges of the island are well dened; increasing the temperature the corners become rounder and the length of less stable edges decreases; at higher temperature the whole cluster is roundish, and straight edges In absence of any force, clusters are at thermodynamic equilibrium with an adatom gas and their shape uctuates around an equilibrium shape. Under external forces, a ux of adatoms leads to the cluster motion. Clusters are thus no more at equilibrium. Their shape can reach a steady state, or can uctuate among dierent congurations. Kuhn, Krug, Rusanen and Dasgupta et al. have carefully studied this behavior with analytical models [3640]. We have studied the shapes of clusters moving under strong forces (E cg larger than 10 -4 E b ), however we have not carefully investigated possible uctuations between dierent congurations. Notice that average forces found in electromigration phenomena are of the order of 10 -6 eV/lattice length [START_REF] Curiotto | 2D nanostructure motion on anisotropic surfaces controlled by electromigration[END_REF], that corresponds to E cg ≈ 10 -6 E b , with E b of the order of 1 eV. The 6-fold symmetry is broken when the simulations are performed with k nsl = 0 and the shape of holes and islands has a three-fold symmetry (see gure 11). We have studied the eects of forces applied in three dierent directions on the nanostructure shapes. Figure 12 shows snapshots of the simulations of islands and holes moving under the eect of these forces. Each snapshot corresponds to a representative island or void shape moving under a constant force. We consider that a cluster has reached a steady state shape when three shape descriptors do not signicantly evolve anymore, as detailed in the supplementary material. We at rst discuss the shapes of voids. When the force is towards the [1-10] direction (towards left in gure 12, rst panel, second line), the shapes seem pushed towards the facet of the right, i.e. towards the advancing edge. The advancing edge is well faceted and its length increases with the force value. For very strong forces the voids are elongated and attened on the right facet. Loosely bound and kink atoms are easily removed from the front under the eect of the force. The unzipping of the atomic rows by removal of the kinks (detachment of a kink atom leaves another kink atom that will also be detached) leads to a displacement of the front edge. This unzipping is only stopped at the corners of the cluster. The clusters tend to keep a convex shape to maximize the total number of bonds between atoms. However the force leads to elongation of the advancing edge and thus of the whole shape, because also atoms that do not belong to the front edge can detach. When this happens, atomic rows of that edge unzip; when the unzipping reach the front edge, this results elongated of one unit length towards the edge that has unzipped. Because of the force towards the back, the detached atoms accumulate at the back edge and only rarely they move against the force, towards the front edge. The cluster thus keep an elongated front edge because of a kinetic eect. A similar elongation, perpendicular to the motion, is observed when the force is in the direction (third panel in gure 12). In this case the elongation is along the top facet, that corresponds to the advancing edge. Applying a force in the [-1-12] direction (second panel in gure 12), the elongation perpendicular to the motion is no more observed. Because of the triangular symmetry, an edge at the advancing front of the void would correspond to a A step, less stable than B steps. Under the eect of the force, atoms are removed more easily from A steps than from B steps. The lateral facets of the void (B steps) in the gure are thus enhanced. The void moves by unzipping of the atomic rows from the lateral facets as previously explained, the detached atoms diuse towards the back edge where they accumulate, and the void shows a strongly faceted triangular shape. Latz et al. [START_REF] Latz | Simulation of electromigration eects on voids in monocrystalline Ag lms[END_REF] simulated the shape of 3D holes in Ag(111). They found elongation of the holes along the force axis and formation of a tip at the hole advancing front. The dierence of their results with ours is due to the dierent mechanism of atomic motion: in their simulations atoms diuse along the edges, by periphery diusion, while in our case atoms diuse inside the hole, by terrace diusion.

Cluster shape dependence on temperature

Shapes of clusters under strong forces

For islands the mechanisms are similar to those described for voids, but in this case atoms are removed from the back edge and accumulate to the front edge that thus advances. Under very strong forces, islands can also break down in smaller islands that advance together. While in voids atoms arriving at the back are stuck because of the force, in islands the atoms arriving to the front are pushed by the force towards the back moving along the periphery of the cluster. The less stable facets decrease in size or are rough, because atoms can be removed easily from them by the force. The most stable facets tend to be attened because of the unzipping of atomic rows previously discussed for voids. Like voids, also islands tend to elongate along the most stable facets of the front edge. However, increasing the temperature, i.e. decreasing the step energy anisotropy, the island shape tends to elongate parallel to the motion, as shown in gure 13 for a triangular islands moving in the [-1-12] direction at three dierent temperatures. The elongation parallel to the force is more dicult for islands with the front facet perpendicular to the motion direction (those of the second panel in gure 12) and it is never observed for voids. The increased temperature decreases the anisotropy of the cluster shape, the tendency to form facets decreases and atoms also diuse in the force direction along the island periphery, that gives a shape elongated in the force direction. In voids, the atoms do not have to displace along the void edges to follow the force direction and thus the elongation due to a periphery diusion eect is absent.

Conclusions

We have used a KMC model to highlight the diusion behaviors of voids and islands in a hexagonal lattice. The diusivity of small clusters does not depend on the size in a continuous way and a perfect size eect is clearly observed. Most of the results obtained on square lattices are also valid on hexagonal lattices, however subtle changes are evidenced. For instance, since hexagonal islands have higher symmetry than square ones, the oscillations of the diusion coecient that exist for the smallest sizes are weakened in hexagonal lattices. We also evidence the peculiar behavior of the diusion coecient of small voids: when the void size increases, the diusivity at rst increases then reaches a plateau and nally decreases. We show that this behavior is not specic to hexagonal lattices but also exists (but to the best of our knowledge has not yet been described) on square lattices. We also show how electromigration adds a bias to the atomic jumps, how the cluster velocity depends on the force strength and on the cluster size and how the cluster shapes change with the temperature and the force direction. Finally, we establish approximated analytical expressions of the diusion coecient of clusters that gives values in quite good agreement with the values extracted from the KMC simulations. These developments lead to a better understanding of the atomic events at the origin of the behavior of diusing clusters and thus to the capability of controlling the motion of nano-objects on a surface. This appendix completes the details of the model used in the simulations. In particular we explain the binding energy of the dierent kind of atoms considered. Every surface atom has three neighbors in the subsurface layer (nsl 1 , nsl 2 , nsl 3 ) (see gure 1b). For each nsl we count the number of neighbors (nn nsl ) in the surface, value variable between 1 and 3. Thus, nn is variable between 3 (for an isolated adatom with 3 nearest neighbors in the underlying layer) and 9 (surface atoms surrounded by 6 in-plane nearest neighbors and 3 atoms in the underlying layer). The binding energy of each surface atom is:

E bin = E b • (nn nsl1 -k nsl • nn 2 nsl1 + nn nsl2 -k nsl • nn 2 nsl2 + nn nsl3 -k nsl • nn 2 nsl3 ) (A.1)
Where k nsl is a parameter to be adjusted. This expression allows to dierentiate the energy of atoms according to their environment and in particular reproduce the difference between A and B steps. When k nsl is 0, atoms at steps A and B have the same binding energy and the clusters have hexagonal shape. In the text, energy values are expressed as multiples of E b , for instance an isolated adatom, orange in gure A.14, has E bin = 3E b . Atoms in kink positions, grey in gure A.14, have E bin = 6E b .

Atoms surrounding a vacancy have E bin = 8E b (values given for k nsl = 0). In the simulations, with the aim to reproduce general behaviors of clusters diusing on a hexagonal lattice, we have used E b = 1 (eV). For specic systems, this value could be changed. For instance, for Au, a value of Eb=0.64 eV (see [START_REF] Curiotto | Surface diusion of Au on √ 3 × √ 3Si(111)-Au studied by nucleation-rate and Ostwald-ripening analysis[END_REF]) could be used. However, as the model is very simple, it cannot be used for accurate quantitative predictions. To obtain shapes with a three fold symmetry (like on Pt(111) [START_REF] Michely | Visualizing the electron scattering force in nanostructures[END_REF] or Si(111) [START_REF] Akutsu | Statistical mechanical calculation of anisotropic step stiness of a two-dimensional hexagonal latticegas model with next-nearest-neighbour interactions: application to Si(111) surface[END_REF]) we have used k nsl values between 0.05 and 0.18. The ratio between the binding energies of atoms at B and A steps results less than 0.9, a reasonable value for (111) fcc and diamond systems [START_REF] Michely | Visualizing the electron scattering force in nanostructures[END_REF]44].

Appendix B. Cluster energy

Expression 2 used in section 3.1 does not directly lead to E cluster = E kink but can be useful to further discuss our results. Clusters tend to be compact, and the boundary atoms have j=2 (E = 5E b , atoms adsorbed at steps, green in gure A.14), j=3 (E = 6E b , kink atoms, grey in gure A.14) or j=4 (E = 7E b , step atoms, red in gure A.14). The cluster energy we nd (≈ 6E b ) does not correspond to the energy of atoms adsorbed at the cluster steps (5E b ) because the displacement of these atoms implies only a uctuation of the position of the cluster around a certain position. An important displacement of the cluster can instead be obtained when kinks are removed (E j = 6E b ), as they lead to the unzipping of entire atomic rows. The boundary atoms constituting the cluster edges (j=4) are many but they are strongly bound (7E b ) and, because of the exponential term in the concentration expression, their contribution to the adatom concentration or to the eective cluster energy is small. However in voids, where the boundary is concave, they could play a role and could explain why the cluster energy found in the simulation is slightly higher than the kink energy.

In order to check a practical example, we have counted the dierent N j of a cluster with R=10 at kT=0. [START_REF] Salo | Island migration caused by the motion of the atoms at the border: size and temperature dependence of the diusion coecient[END_REF]

Figure 1 :

 1 Figure 1: a: Schematics of lattice positions of a single plane that can be either empty (white) or occupied by atoms (yellow). Possible jumps are shown with arrows, the red dashed arrows represent jumps at the box edges with boundary conditions. b: side view and top view of atoms in our model. The small yellow circles represent atoms of the top layer, while large grey circles are atoms under the top layer. Atoms at steps A and B, like the red and green atoms, have the same number of neighbors but they have a dierent local arrangement (the red atom has only one neighbor in the bottom layer at the exterior of the step, while the green atom has two neighbors in the bottom layer, see the short black lines). The atoms under the top layer that are neighbors of the red and the green atoms are indexed nsl1, nsl2 and nsl3 (see appendix A).

  2D voids and islands are clusters of vacancies and adatoms respectively. Their edges uctuate with time, because atoms can detach and re-attach somewhere else. The diusion of clusters on a surface has been investigated by dierent authors experimentally (some prominent examples are [1619]) and theoretically. As

Figure 3 :

 3 Figure 3: a: Diusivity of islands (black square) and voids (red empty circles) as a function of their radius, in logarithmic scale (kT=0.4 E b ). b: Logarithm of the cluster diusion coecient as a function of 1/(kT) for dierent islands and a void, as dened in the gure. The simulations of cluster diusion are carried out with k nsl =0 (see appendix A). Simulations with dierent E ad have been carried out to evaluate its eect on D C .

Figure 4 :

 4 Figure 4: Cluster diusion coecient as a function of the number of atoms in simulations with kT=0.2. The inset shows the same graph in log-log scale to highlight that D ∝ N -β holds only for large N.

Figure 5 :

 5 Figure 5: Low-energy shapes of clusters with a xed number of atoms. Clusters with a compact, and thus very likely, shape can be drawn for dierent sizes, not only for N=7, 19 and 37 (that have compact, 6-fold symmetric shape). Blue and green atoms are kinks, with 3 in-plane neighbors; orange atoms (belonging to edges) have 4 in-plane neighbors; the red atom (adatom at a step) has only 2 in-plane neighbors. For N=20 coupled kinks are light blue, normal kinks are dark blue.

  For instance, the compact, low energy conguration of the cluster with N=12 has 3 step atoms (highlighted in orange in gure 5); a cluster with N=13 is larger and should thus have lower D C , but its compact conguration has 4 step atoms and can thus account for D 13 C > D 12 C . The change of number of step atoms in compact shapes explain almost all the observed oscillations. An exception is the compact conguration of the cluster with N=14, that has the same number of step atoms of the cluster with N=13. Thus the large decrease of D C between N=13 and N=14 is due to the general decrease of D C with the size and to the "coupled kinks" eect.

Figure 6 :

 6 Figure 6: a and b: Histograms of the occurrence of dierent congurations for clusters with N=12 and N=13. The insets show the histograms for N=19 and N=20. c: The red circles show the deviation of D C from a curve t of the D C data in gure 4. The aim of these points is to highlight the oscillations of D C due to the nite-size eect. The blue squares represent the number of probable congurations for each cluster dimension.

Figure 7 :

 7 Figure 7: The dark squares show the cluster diusion coecient D C as a function of the cluster size N (number of empty positions) in log-log scale. The diusion coecient of a vacancy is also shown as a red circle. Notice that both the x and the y axis have a break. D C generally increases with the cluster size up to a plateau around N=50 and then decreases. Certain sizes show signicant deviation from the general trend. Three compact congurations for voids with N=6, 7 and 8 are shown in the bottom of the gure. The boundary atoms (yellow) have dierent jump energies (values written in E b units).

4 DISPLACEMENT

 4 increases, reducing Γ C . Therefore for large voids D C decreases when the size increases. For comparison, in very small islands, the jump rates of border atoms decrease when the island size increases (for instance an adatom has high jump frequency, exp -3E b

Figure 8 :

 8 Figure 8: a: Distribution of the nal positions of a diusing vacancy, after 500 simulations. The grey and purple areas correspond to Ecg=0, and 0.03 respectively. The purple and grey probability distributions are partially overlapped. The starting position of the vacancy is the same for the two sets of simulations and marked by X. b: Velocity of a vacancy for dierent force (and thus Ecg) values. The y axis is in logarithmic scale, the velocity depends exponentially from Ecg. The inset shows that for small Ecg the vacancy velocity depends linearly from the force. kT=0.1E b .

Figure 11 :

 11 Figure 11: Shapes of a stationary 2D island (white) on a (111) surface at dierent temperatures. At low temperatures the images shown correspond to the average of several shape congurations.

5 CLUSTER

 5 SHAPES IN A LATTICE WITH 3-FOLD SYMMETRY 10 are no more observed; at very high temperatures the cluster disintegrates in dispersed atoms. A similar behavior is observed for 2D voids.

Figure 12 :

 12 Figure 12: Shapes of a 2D island (rst line) and a 2D hole (second line) on a (111) surface under a force that adds a bias to diusion. The force is towards the left, the top and the bottom of the images for the three panels respectively. The force acting on the shapes implies an energy change from left to right Ecg=0.0001, 0.001, 0.01, 0.1 E b . kT=0.1 E b .

Figure 13 :

 13 Figure 13: Shapes of a 2D island on a hexagonal lattice under a force implying an energy change Ecg= 0.01 E b to the atomic jumps, that adds a bias to the atomic motion in the [11-2] directions. The temperatures are kT=0.05, 0.1 and 0.15 E b for the shapes from left to right.

Figure A. 14 :

 14 Figure A.14: The atoms at an island edge can have j=1 (blue), 2 (green), 3(grey), 4(red) or 5(black) in-plane neighbors. Within our model, when k nsl = 0, their binding energy is respectively 4, 5, 6, 7 and 8 E b . An isolated adatom is shown in orange (binding energy 3E b ). Other atoms of the island are yellow, while empty positions are white.

Table 2 :

 2 . Averaged over many cluster congurations, N 2 =2 N 3 =18 N 4 =32. With these values, we nd where 6E b corresponds to the value of E Cluster that we have found. Cluster cluster diusion energy N number of atoms in the island or of empty positions in the void N edge average number of atoms of a cluster edge N boundary total number of atoms of a cluster perimeter v velocity of a vacancy displacing under the eect of a force nsl 1 , nsl 2 , nsl 3 indexes identifying three dierent atoms of the substrate that are neighbors of a certain atom in the surface nn nsl1 number of atoms in the surface that are neighbors of the substrate atom nsl 1 summary of the symbols used in the text

	< exp(	-E j kT	) >= exp(-6E b /kT )

quères, Diusion Rates of Cu adatoms on Cu(111) in the presence of an adisland nucleated at fcc or hcp sites, Physical Review B 72 (2005) 115402.

Acknowledgments

We thank B. Ranguelov and M. Michailov for helpful suggestions. We also thank A. Saul and O. Pierre-Louis for stimulating discussions. We acknowledge nancial support from the French Region Sud Provence-Alpes-Cote d'Azur (MELOPEE project).

Declarations of interest: none