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

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

Ali El-Barraj
email: elbarraj@cinam.univ-mrs.fr

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

Olivier Pierre-Louis
email: olivier.pierre-louis@univ-lyon1.fr

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

2D nanostructure motion on anisotropic surfaces controlled by electromigration

Introduction

The manipulation of nanostructures is an important challenge for nanoelectronics and nanofabrication technologies. External elds and thermal gradients are natural candidates to move and control the nanostructure position on a surface. Particularly promising is the use of electric elds to induce mass transport in specic directions, phenomenon known as electromigration [START_REF] Tao | Visualizing the electron scattering force in nanostructures[END_REF]. This phenomenon has been widely studied in another context since the sixties: an electric current owing through a circuit, if not controlled, can lead to the formation of hillocks and voids resulting in breakdown of metallic wires [START_REF] Blech | Electromigration in thin Al lms[END_REF]. The elementary process of electromigration of an atom on a surface (adatom) is well known: the electric current biases random diusion inducing mass transport in specic directions. The phenomenon is due to either kinetic momentum exchange between conducting electrons and adatoms (wind force) or because of electrostatic forces (direct force). The atoms are thus described as carrying an eective charge that results from the contribution of these two forces. The drift of adatoms along facets or terraces does not necessarily result in a motion of nanostructures. Indeed, if the number of atoms arriving and attaching to a step edge is equal to the number of atoms leaving the step, then that edge does not move. An imbalance between the two uxes is needed to allow for the motion of the edge. This imbalance can for instance originate from an Ehrlich-Schwoebel barrier (an inequality of atomic attachment and detachment properties at steps [START_REF] Ehrlich | Atomic view of surface self diusion: tungsten on tungsten[END_REF][START_REF] Schwoebel | Step motion on crystal surfaces. II[END_REF]) or a dierent diusion coecient at the two sides of the edge. While it has been neglected in most studies, the crystalline anisotropy of the surface may play an important role on adatom motion and therefore on nanostructure motion.

To investigate the motion of nanostructures under electromigration on anisotropic surfaces we have chosen Si(100) as a model system. The eect of electromigration on the pairing of steps and on step bunching in this system has been studied by several authors [511]. The anisotropy of the diusion coecient has been found to play an important role, however the eect of the coupling between the electric current direction and surface anisotropy on the motion of nanostructures (like one-atom thick islands or one-atom deep holes) still needs to be adressed.

The experimental conditions and instrumental tools required to fully characterize the motion of nanostructures on mesoscopic length are strongly constrained: (i) The material surfaces must be at enough to avoid obstacles on the nanostructure trajectory; (ii) an in-situ and in real-time multiscale experimental technique is necessary to obtain simultaneous information on the kinetics of the nanostructure motion and on the shape evolution; (iii) a large temperature range must be explored in order to get information on the thermodynamic behavior of the moving nanostructures. Following the evolution of scientic equipments, rst in-situ experimental studies on electromigration have used Reection Electron Microscopy (REM), while most recent investigations are based on Scanning Tunneling Microscopy (STM). However both approaches have their own limit. In REM experiments, due to the glancing conditions for the incident electron beam, the images are severely shortened along one axis. Therefore the shape evolution of a moving nanostructure cannot be studied [START_REF] Yagi | Reection electron microscopy: studies of surface structures and surface dynamic processes[END_REF][START_REF] Müller | Low distortion reection electron microscopy for surface studies[END_REF]. Con-2 ventional STM experiments only give access to the low temperature regime (<500 K) and do not allow to study fast moving and/or large nanostructures. In order to overcome these limitations we have used Low Energy Electron Microscopy (LEEM), that allows one to observe surfaces in-situ and in real time with high temporal and spatial resolution, and in a large temperature range (1000-1500 K). Furthermore, we have used lithographically prepared surfaces to obtain atomically at terraces with micrometric sizes on which 2D nanostructures may be displaced on large distances without interactions with steps.

In order to support the experimental results and obtain more insights in the motion mechanism, we have developed a Kinetic Monte Carlo model based on both anisotropic and electromigration-biased diusion at surfaces. Nanostructure motion biased by electromigration on anisotropic surfaces may be understood on the basis of simple schemes, which nevertheless may lead to unexpected motions. For instance, since electromigration aects mass transport, one could expect that a nanostructure moves along the current axis, either in the electric eld direction or opposite to it, depending on the nature of the nanostructure (2D hole or 2D island) and on the main mechanism at work (surface diusion, periphery diusion or attachment/detachment). Counter-intuitively we show that, due to the anisotropy of the surface diusion coecient, nanostructures can also move perpendicularly to the electron ow. We have derived a simple analytical expression to describe the motion of 2D nanostructures under an electric eld in case of anisotropic diusion. Furthermore, since working in a large temperature range, an eective activation energy for the motion of 2D nanostructures has been obtained. This energy is given by two contributions, a diusion energy and an adatom formation energy.

Experimental details

Single-crystal B-doped (p-type) silicon wafers with 500 µm thickness, nominal orientation (100), and a miscut of ±0.01 • were supplied by Siltronix. Circular holes 7.5µm deep, with 300 µm diameter separated by 100 µm were prepared by standard optical lithography using a hexagonal mask. The holes were made to interrupt the regular train of steps and allow formation of large terraces, as suggested by [START_REF] Lee | Formation and stability of large step-free areas on Si(001) and Si(111)[END_REF]. Then, the wafers were cut in laths 3 mm wide, 12 mm long, and cleaned with acetone and ethanol. According to the sample, the long side is aligned along the [011] or the [010] direction. They were mounted on a special sample-holder (designed by Elmitec Gmbh) which allows to apply an electric current through the sample by means of two Mo electrodes clamped to the extremities of the sample. The current direction was parallel to the longer side of the sample and the applied electric current was typically within the range 0.5-3 A. The samples were introduced in Ultra High Vacuum (UHV) and degassed for several hours at about 1100 K. Then they were ashed above 1500 K to remove the surface oxide and the rst Si layers. The temperature was decreased and kept for some minutes at 1300 K, where fast surface diusion, slow Si evaporation and electromigration [START_REF] Latyshev | Peculiarities of step bunching on Si(001) surface induced by DC heating[END_REF] allow to increase the size of the terraces up to some µm. The temperatures were measured with a precision of 30 K, using a disappearinglament pyrometer, calibrated with the 7×7 to 1×1 surface phase transition occurring on Si(111) at 1103 K [START_REF] Ino | Some new techniques in Reection High Energy Electron Diraction (RHEED) application to surface structure studies[END_REF]. The in-situ observations were performed in a Low Energy Electron Microscope (LEEM III, Elmitec Gmbh), in dark eld mode, using electron energies ranging between 4 and 6 eV. The images were acquired at rates of 0.1-1 Hz. Single-step holes have been made by evaporation in large terraces, increasing the current through the sample and therefore increasing the temperature (above 1340 K) and observing in real time by LEEM their formation. When the desired hole size is reached, the temperature is lowered to stop evaporation and the sample is kept at a certain temperature and observed to study the hole motion dynamics. Islands can also be created by depositing in-situ less than 1 mono-layer of Si, up to observe the formation of the islands; their number and size can be tuned by changing the substrate temperature, the ux or the deposition time.

Experimental results

The reconstructed surface of Si(100) is arranged in parallel dimer rows running along the [011] or the [01-1] axis: the direction rotates by 90 • from one terrace to the next separated by a monoatomic step (see Figure 1a). The terraces therefore give dierent diraction spots in Low Energy Electron Diraction (LEED). Thus in dark eld LEEM, successive terraces appear alternatively dark and bright. Figure 1b shows dark ellipses, that correspond to one-atom deep holes, in a Si(100) terrace. The long axis of the elliptic holes is oriented along the dimer rows of the upper terrace (schematically drawn with parallel lines in Figure 1). To obtain these one-atom deep holes on an extended terrace, the sample temperature is rst raised to increase the evaporation rate and to nucleate and grow holes. Then the temperature is decreased to reach the desired working conditions (typical temperatures of 1080-1300 K). We have changed the direction of the current with respect to the direction of the dimer rows and observed the motion of nanostructures, in order to study the combined eect of electromigration and anisotropic diusion.

When an electric current is applied through the sample in the [011] direction, the black holes, elongated perpendicularly to the current, move opposite to the direction of the current, as shown in Figure 1c taken 278 s after 1b. Similarly, Figures 1d and1e show the motion of holes elongated in the current direction (dimer rows of the upper terrace parallel to the current). In this case the holes (bright because their dimer rows are perpendicular to the dimer rows of the holes in Figures 1b and1c) move in the same direction as the current.

We have also applied a current in the [010] direction, not aligned along the dimer rows, but at 45 • from them. Figure 1f shows bright holes elongated diagonally, where the dimer rows of the upper dark terrace are oriented from the bottom left towards the top right of the image. Surprisingly, with a current applied as shown in the gure, from left to right, the holes move towards the top (Figure 1g), at 90 • from the electric-eld direction and 135 • from the [011] axis. We have also studied islands instead of holes. As shown in the movie of supplementary material S1, the islands displace opposite to the motion direction of holes with the same dimer row orientation.

The velocity of the motion does not depend on the size of holes or islands. Therefore according to established interpretations [START_REF] Pierre-Louis | Electromigration of single-layer clusters[END_REF], the process of mass transport is limited by diusion of adatoms on terraces and is not inuenced by periphery diusion or attachment detachment kinetics. For identical current-temperature conditions, the velocity of the nanostructure is independent from the orientation of the force (within the experimental error). However, the dependence of the nanostructure motion direction on the current direction is not trivial. In order to explain this dependence we have to take into consideration the atomic structure of the Si surface and its diusion anisotropy. Several authors have found that on Si(100) surface diffusion is much faster along the dimer rows than perpendicularly to them. The ratio between the reported values of parallel and perpendicular diusion coecients ranges between 4 and 1000 [START_REF] Ichikawa | Study of Si(001) 2x1 domain conversion during direct current and radiative heatings[END_REF][START_REF] Mo | Anisotropy in surface migration of Si and Ge on Si(001)[END_REF] depending on the temperature. In the following section we develop a model to understand the nanostructure motion under electromigration on anisotropic surfaces.

Model of nanostructure motion under electromigration

In order to model the behavior of holes and islands under electromigration in Si(100), we use a basis of unit vectors u and w oriented along the [011] (we consider it equal to the direction of the dimer rows inside the nanostructure for both islands and holes) and the [0-11] direction respectively, as shown in Figure 2. The force due to the electric eld has a direction forming an angle θ with u. We consider dimers as mobile species (see the KMC section) and discuss here the case of a one-atom deep hole. Analogous expressions can be derived for one-atom thick islands. As the experimentally measured velocity does not depend on the size of the nanostructure, terrace diusion is considered to be the dominant mechanism of mass transport [START_REF] Pierre-Louis | Electromigration of single-layer clusters[END_REF]. Atoms diuse with a bias due to the sum of two forces: a wind force (momentum exchange between conducting electrons and adatoms) and a direct force (electrostatic force from the applied electric eld on atoms). We assume that the prevailing force F is the direct force, in the same direction of the electric current, as also supposed in [START_REF] Ichikawa | Study of Si(001) 2x1 domain conversion during direct current and radiative heatings[END_REF][START_REF] Doi | Anisotropic diffusion between the step-up and the step-down directions on a Si(001) surface[END_REF] and usually found on semiconductor surfaces [START_REF] Yasunaga | Electromigration on semiconductor surfaces[END_REF]. This force can be written F = Z * • e • |E| • (cosθ u + sinθ w), where Z * is the adatom eective valence, e is the electron charge, and |E| is the magnitude of the applied electric eld.

The surface mass uxes (number of atoms per unit time, per unit length) at the exterior and at the interior of the hole depend on the concentration of diusing species, on their diusion coecient, and on the force acting on them. These uxes can be written as vectors j ext and j int with components (α=u or w):

j i,α = c kT β=u,w F α D i,αβ
where if i=int the diusion tensor reads

D int = D uu = D D uw = 0 D wu = 0 D ww = D ⊥ whereas if i=ext the diusion tensor reads D ext = D ⊥ 0 0 D
c is the dimer concentration, considered constant, k is the Boltzmann constant, T is the temperature. D ⊥ and

j ext = A • (D ⊥ • cosθ • u + D • sinθ • w) j int = A • (D • cosθ • u + D ⊥ • sinθ • w) v h = A • a 2 (D -D ⊥ ) • (-cosθ • u + sinθ • w) v isl = A • a 2 (D -D ⊥ ) • (cosθ • u -sinθ • w)
A = c • Z * •e•|E| kT
D are the surface diusion coecients respectively perpendicular and parallel to the dimer rows. There is no Fick-diusion term proportional to the gradient of the concentration in the mass ux j because we assume that the system is in the diusion-limited regime, where the concentration is xed to its equilibrium value close to atomic steps. This assumption is valid as long as the attachmentdetachment dynamics at the steps is fast enough, a condition which should be observed in our high-temperature experiments where the kink density at steps is large.

If the interior and exterior uxes j are not equal, a net mass ux at the edges of the hole results in an advancement velocity. In steady state regime the hole shape is xed and the hole velocity is enforced by mass conservation:

v h = a 2 • ( j ext -j int ) (1)
where a is an atomic size. If instead of holes we consider one-atom thick islands, the growth direction of the steps is opposite to the previous case, therefore the direction of island motion is opposite to that derived for holes.

The analytical expressions of the j uxes, the hole velocity v h and the island velocity v isl that result from previous equations are reported in table 1. This model emphasizes the dependence between nanostructure (island or hole) motion and diusion anisotropy, since for D = D ⊥ , v h = v isl = 0. Furthermore these equations directly give the angle γ formed by the velocity of the nanostructure with the dimer rows direction u. For this purpose it is sucient to compare the velocity expressions with v = |v| • (cosγ • u + sinγ • w) to extract γ = 180 • -θ for 2D holes and γ = -θ for 2D islands.

The expression of the nanostructure velocity obtained with this model explains all the experimentally observed displacement direction of holes and islands. For instance it correctly predicts that when θ = 45 • holes move perpendicularly to the electric eld (γ -θ = 90 • , see table 2). Figure 3 schematically summarizes the motion directions of holes and islands according to the angle formed between the electric eld and the dimer rows inside the nanostructure. Notice that the model reproduces the experimental observation that the absolute value of the velocity of the center of mass does not depend on the direction of the applied eld. Moreover it does not depend on the nanostructure shapes, when isolated nanostructures are considered (that means nanostructures that do not interact with their local environment e.g. steps). The good agreement between the model and the experimental results suggests that the hypothesis of atoms drifting in the same direction of the electric current is correct.

Eective energy

Figure 4a shows the trajectories of several one-atom deep holes as a function of time. For all the nanostructures the distance from the initial position increases linearly with time (i.e. constant velocity). Only at the end of the motion their velocity decreases because they are disappearing as a consequence of Ostwald ripening (see later). From the above model (see table 1) we can develop the expression of the nanostructure velocity modulus:

|v h | = |v is | = a 2 • c 0 • e -Ec kT • Z * • e • |E| kT • (D 0, • e -E kT -D 0,⊥ • e -E ⊥ kT ) (2)
The dimer concentration c depends on the temperature and on a dimer creation energy E c . The diusion coecients are developed considering the prefactors (D 0, and D 0,⊥ ) and the diusion energies (E and E ⊥ ). If D ⊥ is much lower than D , then we can write:

|v h,is | ∝ E • e -Ec +E kT kT = E • e -E eff
kT kT Measuring the velocity at dierent temperatures, we estimate the eective energy E eff =E c +E (see Figure 4). The slight dierence between the data points obtained from holes and those obtained from islands is negligible. The data points obtained by Métois et al. [START_REF] Métois | Steady-state motion of silicon islands driven by a DC current[END_REF] for island motion at Si vapor pressure equilibrium are also reported for comparison, they are in good agreement with our measurements. Notice that, if E ⊥ >> E , then D ⊥ can be neglected and the slope in Figure 4 corresponds to E eff . Also if E ⊥ ≈ E the slope is close to E , but the intercept with the y axis signicantly changes (therefore the pre-exponential cannot be related to physical constants).

A linear t through the points gives an eective energy of E eff =2.0±0.1 eV. This value is close to that attributed by Keee et al. [START_REF] Keee | Surface self-diusion on Si from the evolution of periodic atomic step arrays[END_REF] to the sum of the surface diusion plus the creation energy of the diusing species (2.3 eV). A similar value is also reported for the eective energy obtained from another phenomenon involving Si(100), i.e. dewetting from a SiO 2 substrate (2 eV [START_REF] Cheynis | Dewetting dynamics of silicon-on-insulator thin lms[END_REF]). If we take the binding energy found at ambient temperature by Swartzentruber [START_REF] Swartzentruber | Kinetics of Si monomer trapping at steps and islands on Si(001)[END_REF] for units at the end of dimer rows, E c =1.0±0.1 eV, equal to the dimer creation energy, then we obtain E = E eff -E c =1.0±0.2 eV. This value is equal to the Table 2: Angles between the electric eld or the hole motion direction with the [011] direction and corresponding gure. Also the angle between the current and the motion direction is reported (γ -θ). activation energy of about 1 eV for dimer motion below 400 K [START_REF] Swartzentruber | Kinetics of atomic-scale uctuations of steps on Si(001) measured with variable-temperature STM[END_REF], and compares reasonably well to the surface diusion energy found from island densities during Si deposition between 300 and 600 K (0.7 eV) [START_REF] Mo | Activation-energy for surface diusion of Si on Si(001) -A Scanning-Tunneling-Microscopy study[END_REF]. We now evaluate the eective valence Z * . Z * can be obtained from equation 2, using the experimental values summarized in table 3. For the calculation, not only the diusion energies, but also the prefactors must be known. Values for both D and D ⊥ have been given by Doi et al. [START_REF] Doi | Activation energies of Si adsorbate diusion on a Si(001) surface[END_REF] who studied surface diusion on Si(100) by measuring the size of denuded zones (distance between a step and the rst island on the terrace) after Si deposition and annealing. Because of the experimental method they used, we rather suggest they actually found D • c and D ⊥ • c (see table 3). We therefore nd 7 < Z * < 25.

The eective charge we found is of the order of magnitude of that estimated by Stoyanov [START_REF] Stoyanov | Electromigration induced step bunching on Si surfaces -how does it depend on the temperature and heating current direction[END_REF] for Si(100) (Z * =1) and consistent with the results of Saul et al. [START_REF] Saúl | Experimental evidence for an Ehrlich-Schwoebel eect on Si(111)[END_REF] who found Z * >1.3. Notice that, while we are condent on our velocity measurements, the eective charge value strongly depends on the order of magnitude of the diusivity [START_REF] Doi | Activation energies of Si adsorbate diusion on a Si(001) surface[END_REF] (D • c) that we have used.

Nanostructures size and shape

As discussed by Hamers et al. [START_REF] Hamers | Nucleation and growth of epitaxial silicon on Si(001) and Si(111) surfaces by scanning tunneling microscopy[END_REF], the nanostructure shapes on Si(100) are determined by the energy dierence of S A steps and S B steps (see Figure 1a). For instance Chadi [START_REF] Chadi | Stabilities of Single-Layer and Bilayer Steps on Si(001) Surfaces[END_REF] calculated 0.01 and 0.15 eV for S A and S B steps respectively. Without electromigration, the equilibrium shape of one-atom thick nanostructures is an ellipse [START_REF] Bartelt | Step capillary waves and equilibrium island shapes on Si(001)[END_REF]. With increasing temperature the aspect ratio (short axis/long axis) increases (see supplementary material S2), and the nanostructures becomes more rounded as also reported in [START_REF] Bartelt | Step capillary waves and equilibrium island shapes on Si(001)[END_REF].

In the experiments, the holes and islands areas change because of a combination of Ostwald ripening [START_REF] Curiotto | Surface diusion of Au on √ 3 × √ 3Si(111)-Au studied by nucleation-rate and Ostwald-ripening analysis[END_REF] (due to atom exchange with the surrounding terraces) and evaporation, that takes place at high temperature. The behaviour of holes and islands is quite dierent. Indeed while 6 2D islands always shrink, 2D holes shrink when they move in the direction of the electric eld (i.e. when θ is 90 • in Figure 2) but grow when they move opposite to the electric eld (θ=0 • ). The area change versus temperature is shown in Figures 5a and5b for holes and islands respectively. Since experimental growth and shrink rates depend on the nanostructure size and on its close environment, we only show average values (for various sizes) of the area change.

|E| [V/nm] |v

h |[nm/s] D • a 2 • c [s -1 ] D ⊥ • a 2 • c [s -1 ] T [K]
The dierence between holes and islands is strongly associated to the anisotropy of surface diusion. Indeed, since Ostwald ripening of a nanostructure is proportional to the diusion rate on the terraces surrounding the nanostructure [START_REF] Curiotto | Surface diusion of Au on √ 3 × √ 3Si(111)-Au studied by nucleation-rate and Ostwald-ripening analysis[END_REF], it is more important when diusion outside holes/islands is high (i.e. when the dimer rows outside the nanostructure are in the direction of the electric eld, θ=90 • ). When diusion outside the nanostructures is limited (θ=0 • ), Ostwald ripening takes place at a lower rate so that evaporation dominates over ripening. For islands, both Ostwald ripening and evaporation lead to shrinking (see Figure 5b). However, as ripening is faster when diusion on the surrounding terraces is faster (θ=90 • ), islands with dimer rows perpendicular to the electric eld shrink faster than those with dimer rows parallel to the eld.

For holes, evaporation leads to growth and Ostwald ripening to shrinking. When diusion outside the holes is hindered (θ=0 • ), evaporation prevails and the holes grow (full squares in Figure 5a). When Ostwald ripening dominates (θ=0 • ), the holes tend to disappear (empty squares in Figure 5a). However at high temperature the data are very scattered because the two phenomena are competing and the error on the measured area change is large.

Kinetic Monte Carlo model

We have developed a simple model which is able to reproduce our experimental results using a standard kinetic Monte Carlo method (KMC). The model is inspired from that of Ghosh and Ranganathan [START_REF] Ghosh | Submonolayer growth study using a solid-on-solid model for 2 x 1 reconstructed surfaces of diamondlike lattices[END_REF] and from the description of electromigration developed by Pierre-Louis and Einstein [START_REF] Pierre-Louis | Electromigration of single-layer clusters[END_REF]. In order to describe the (100) reconstructed surface formed by dimer rows we use a simple cubic solid-on-solid lattice with periodic boundary conditions. Every position of the lattice represents a column of diusing units and has a certain height. Our LEEM analysis does not give information on the nature of the diusing units on Si(100), which could be atoms, dimers or vacancies. Exchange processes between the surface and the bulk could also take place. Adatoms have been rarely observed by STM because either they diuse very fast or because single adatoms are very unfavorable, as concluded by Huang and Allen [START_REF] Huang | Diusion of Si atoms and dimers on Si(100)[END_REF] (because dimer formation reduces the system energy by 1 eV). Also Bartelt et al. [START_REF] Bartelt | Ostwald ripening of twodimensional islands on Si(001)[END_REF], Tromp and Mankos [START_REF] Tromp | Thermal adatoms on Si(001)[END_REF] and Swartzentruber [START_REF] Swartzentruber | Direct measurement of surface diusion using atom-tracking Scanning Tunneling Microscopy[END_REF] suggest that the diusing unit is a dimer. We thus consider the diusing units to be dimers.

A jump of a diusing unit depends on the neighbors of the unit and on the destination site. It takes place with probability proportional to

exp( -E b -E d -E m kT )
where k is the Boltzmann constant, T is the temperature and the three dierent energies E b , E d and E m are explained in the following paragraphs. (•) The binding energy E b of the diusing unit reads:

E b = E V + n dr • E dr + n pr • E pr
Within this expression the diusing unit is bound (i) to the underlying surface with energy E V (E V = E =1 eV corresponds to the order of magnitude of the diusion energy of isolated units along dimer rows [START_REF] Swartzentruber | Kinetics of atomic-scale uctuations of steps on Si(001) measured with variable-temperature STM[END_REF]); (ii) to up to two other dimers of the same dimer row, (0≤ n dr ≤2), with energy E dr ; and (iii) to up to two other dimers of parallel dimer rows, 0 ≤ n pr ≤ 2, with energy E pr . More details on the binding energies and a schematics of the column positions in the model are reported in the supplementary material S3. (•) E d is an energy introduced to take into account diusion perpendicular to dimer rows. It is added only if the arrival site of the diusing unit is on top of a dierent underlying dimer row. In comparison with the analytical model,

E V + E d = E ⊥ . We have used E d =0.1 eV.
(•) E m is the energy change due to the electric eld. It depends on the arrival site of the diusing unit and on the direction of the electromigration force (see schematics of supplementary material S3) according to the expression

E m = |F | • a • cos( π 2 • b -δ).
a is a lattice length that we take equal to 1 in all directions, |F | is the electromigration force acting on a jump of unit length, previously written as Z * • e • |E|, and b is an integer between 1 and 4 dening the arrival site, in anticlockwise order (b=1 for a jump to the right, b=2 for a jump to the top...); δ is the angle between the electromigration force and the x axis. The order of magnitude of the experimentally applied electric eld |E| is 10 -7

V olt lattice-length . We have evaluated an eective valence Z * of the order of 10 elementary electron charges, therefore we should use a force lower than 10 -6 eV lattice-length . However that value is too low to observe any change due Figure 5: a:Average area change as a function of temperature for holes moving in the direction of the electric eld (red empty squares) and for holes moving opposite to the electric eld (black full squares). b:Average area decrease for islands moving in the direction of the electric eld (full squares) and for islands moving opposite to the electric eld (empty squares). Empty and full squares correspond to experimental results, while lines are a guide for the eyes. The error bars correspond to statistical errors that increase with temperature because of evaporation. to electromigration in the simulations in reasonable time, therefore we have used |F |=10 -4 eV lattice-length . Many details of the (100) Si surface are neglected, such as the dimer buckling, the dierent possible diusion sites (on top of dimers or in between dimer rows), or Ehrlich-Schwoebel barriers. We also consider the same eective charge for isolated dimers and dimers at steps. However the ingredients we use are sucient to reproduce many features of the experimental behavior of islands and holes.

Kinetic Monte Carlo simulations

The shapes of the simulated nanostructures change as a result of the atom diusion, but are approximately elliptical and their aspect ratio depends only on the values of the binding energies E dr and E pr . Isotropic shapes are obtained when E dr =E pr . Simulated nanostructures or terraces have S A steps atter, with less kinks than S B steps, as observed by Scanning Tunneling Microscopy (see for instance [START_REF] Swartzentruber | Kinetics of atomic-scale uctuations of steps on Si(001) measured with variable-temperature STM[END_REF]). Figure 6 shows dierent types of holes and islands displacing in dierent directions under the eect of the same electromigration force (directed from right to left). For instance Figure 6a and6b show two snapshots of a hole that moves from left to right. All the obtained motion directions correspond to those experimentally observed and predicted by the analytical model. Under electromigration but without diusion anisotropy (E d =0) holes and islands do not migrate in dened directions but only diuse randomly. This conrms the previously obtained expression for the displacement of nanostructures, that depends on the dierence of diusion coecients parallel and perpendicular to the dimer rows (equation 2).

Figure 7a and 7b show the displacement of an island when the force is directed at θ=150 • . Figure 7c shows the trajectories of islands for dierent orientations of the electromigration force. These simulations reproduce very well the experiments and the prediction of the analytical model.

The velocities obtained by simulating the motion of islands at dierent temperatures are plotted as a function of 1/(kT) in Figure 7d, with a graph similar to that of Figure 4b. The slope of the linear t provides an activation energy of 1.44±0.02 eV. This value corresponds to the sum of two energies used in the KMC model: (i) the diusion energy in the fast direction (E V = E =1 eV) and (ii) the lowest energy to create adatoms (E c ), that is the binding energy of a dimer in a kink position in the KMC model (E pr +E dr =0.186+0.26 eV). In order to obtain the experimental eective activation energy (2 eV), the kink energy in the simulations must be 1 eV (2-E V =1). If we use E pr =0.42 eV and E dr =0.58 eV, the eective activation energy found with the simulations is equal to the experimental one.

We have shown that our KMC model is able to reproduce the behaviors of hole and island motion. In order to test the predictive character of the model, we now use it for the analysis of more complex geometries in regimes where simple analytic predictions are not available.

Edge instabilities induced by electromigration

The instability of step edges appears as a natural candidate to challenge and test the predictive character of the KMC model. This instability is observed below 1100 K in both one-atom thick islands and holes. The shape of the holes and islands is then found to deviate from ellipsoids, and ngers form on the side where electromigration is down-step (Figure 8a). This eect is more visible on straight steps that develop crests and valleys whose amplitude increases with time. In the absence of electromigration, S A steps are straight, whereas S B steps are rough because of their dierent step stiness [START_REF] Chadi | Stabilities of Single-Layer and Bilayer Steps on Si(001) Surfaces[END_REF]. Under electromigration both types of steps are destabilized. The step on the right in Figure 8a (close to a S B step) develops crests pointing roughly horizontally, whereas the step on the left (close to a S A step) has crests pointing at about 45 • with respect to the dimer row orientation (which are vertical on the bright terrace on the gure).

Such instabilities are reproduced by KMC simulations when both electromigration and anisotropic diusion are implemented (Figure 8b). More precisely Figure 8b shows a terrace (red) with horizontal dimer rows one layer above another terrace (blue) with vertical dimer rows. At the beginning of the simulation (not shown) the steps of the red terrace are S B at steps. Because of an electromigration force oriented from left to right, the terrace moves to the right and the advancing edge (the right one, where the current is in the step down direction) is destabilized. Without electromigration or without diusion anisotropy (E d =0) there is no formation of crests and valleys (Figure 8c shows the simulation of a terrace under electromigration but with isotropic diusion, E d =0).

This instability does not correspond to that described by Bales and Zangwill [START_REF] Bales | Morphological instability of a terrace edge during step-ow growth[END_REF], which originates in an increased ux from the lower side of an advancing step. Indeed, we rather have an increased ux of adatoms to the step from the upper terrace which should lead to step The black terrace at the left of the central grey region is lower than the bright terrace on the left. The step between them is some degrees away from a perfect S A step. The black terrace at the right of the grey central region is higher than the bright terrace on the right. The step between them is only some degrees away from a S B step. The height prole is on the top of the image. b: KMC simulation of step instabilities. The red terrace has dimer rows oriented horizontally and is one monolayer higher than the blue one. The electromigration force is oriented from left to right and has a value of 0.01, while all the other parameters have the values described in the text. The red terrace develops crests and valleys at the advancing edge. c: same conditions as in Figure a, but the diusion is isotropic (E d =0); in this case the edges do not move and do not develop crests and valleys.

attening in the Bales and Zangwill picture. Sato et al. [START_REF] Sato | Step wandering induced by the drift of adatoms in a conserved system[END_REF] studied the linear stability of steps on Si(111) vicinal surfaces under electromigration with isotropic diusion. They found an instability of step edges when the adatoms drift in the step down direction, and when the steps are not permeable. Step permeability is dened as the propensity of diusing species to cross easily the steps. In our experiments and simulations, the instabilities also take place when the drift of dimers is in the step down direction. However, steps on Si(100) are known to be permeable [START_REF] Tanaka | Step permeability and the relaxation of biperiodic gratings on Si(001)[END_REF], and the absence of instability observed in KMC simulations with isotropic diusion also suggests that steps could be permeable in the KMC model. As a consequence, the prediction of Sato et al. does not correspond to our observations. In a dierent study, Sato et al. show that on Si(111) a down-step drift destabilizes the steps if there is a dierence of diusion coecient between the lower and the upper terrace [START_REF] Sato | Drift-induced step instabilities due to the gap in the diusion coecient[END_REF]. The alternated anisotropy of terraces on Si(100) surfaces, which induces a dierence in the mass uxes on the two sides of the step, could play a role analogous to that played by the dierence of diusion coecient. Further analytical modeling along the lines of ref [START_REF] Sato | Step wandering induced by the drift of adatoms in a conserved system[END_REF][START_REF] Sato | Drift-induced step instabilities due to the gap in the diusion coecient[END_REF] would be needed to conrm the analogy. Nevertheless, the simple framework of our KMC simulations allowed us to reveal that diusion anisotropy triggers a novel step instability caused by diusion anisotropy that has not been identied in the previous experimental or theoretical literature. Our results call for further theoretical analysis and simulations to unravel the role of diusion anisotropy on migration-induced step instabilities.

Conclusions

We have shown that under electromigration the diusion anisotropy can aect the motion direction of nanostructures. This concept is demonstrated on Si(100) by insitu and real-time observations performed by LEEM of the motion of one-atom thick islands and holes under an electric eld oriented along dierent directions. For a given θ, angle between the electric eld and the dimer rows inside the nanostructure, one-atom deep holes move at an angle of 180 • -θ while one-atom thick islands move at angle equal to -θ. The motion of nanostructures on Si(100) is possible thanks to the dierence of diusion coecient parallel or perpendicular to dimer rows. From the velocity measurements of the motion of one-atom thick nanostructures we obtain an eective activation energy of 2.0±0.1 eV given by the sum of a diusion energy and a creation energy of the diusing species. Electromigration at low temperatures when the electric eld is mainly perpendicular to the S B steps tends to de-stabilize S B steps, leading to formation of crests and valleys. Finally, we have developed a Kinetic Monte Carlo model to simulate electromigration biased diusion on Si(100). The simulations reproduce all the experimental observations and conrm the mechanism of motion proposed.

3

 3

Figure 1 :

 1 Figure 1: a: Illustration of the Si(100) surface structure, with dimer rows running along the [011] or the [0-11] axis: the direction changes by 90 • in alternating terraces. S B steps, perpendicular to dimer rows, are rough, while S A steps, parallel to dimer rows, are atter because of their higher step energy [16]. b-g: dark eld LEEM images of one-atom deep holes in Si(100) (E=4.5eV). Under electromigration (current directed towards the right for all the images of the gure) the holes displace in dierent directions according to their orientation. u and w are shown to understand the model explained in section 4 and the schematics of Figure 2. The dimer rows inside the holes are oriented in the direction of u, while outside the holes they are parallel to w. Depending on the orientation of the dimer rows inside the holes, u and w are either in the [011] and [0-11] directions (b and c) or in the [0-11] and [0-1-1] directions (d and e). b-c: the holes with dimer rows running parallel to the force displace to the left, opposite to the electromigration force; d-e: the holes with dimer rows running perpendicularly to the force displace in the same direction of the force; f-g: when the holes have dimer rows at 45 • from the force direction, they displace perpendicularly to the force. The step bunches on the top right of b, c, d and e and on the left of f and g are the xed reference regions (marked with a red cross). All the elliptical structures shown in this gure are holes but they are dark or bright according to the direction of their dimer rows and according to the diraction spot used in the dark eld image.

Figure 2 :

 2 Figure 2: Schematics of a one-atom deep hole in Si(100). The dimer rows inside the hole are oriented horizontally, in the same direction of the u vector. The force exerted by the electric eld is oriented along the blue arrow and forms an angle θ with the u vector. The hole displaces in a direction forming an angle γ with u.

Figure 3 :

 3 Figure 3: Graphical summary of the velocity direction of holes and islands under the eect of an electromigration force. Arrows with full lines represent directions of the force due to the electric eld. They form an angle θ with u equal to 0 • , 30 • , 45 • , 60 • , and 90 • for red, green, blue and yellow arrows respectively. The dashed arrows represent the velocity directions of holes (left) and islands (right) corresponding to force arrows with the same color.

5

 5

θ

  • γ • γ -θ • experiment 0

Figure 4 :

 4 Figure 4: a: Position as a function of time of seven holes at 1150 K.b: logarithm of the experimentally measured hole and island velocity multiplied by kT and divided by the electric eld as a function of 1/(kT). Black square and green triangles refer to our measurements of island and hole motion respectively. For comparison, we also report the data obtained by Metois et al.[START_REF] Métois | Steady-state motion of silicon islands driven by a DC current[END_REF] as red circles. A linear t through the data gives an eective activation energy of E eff =2.0±0.1 eV.

Figure 6 :

 6 Figure 6: KMC snapshots of successive positions of nanostructures that move under the eect of a force to the left. The black lines represent the direction of the dimer rows in the images, the arrows show the motion directions. For this calculations kT=0.05; (a-b): hole with horizontal dimer rows; (c-d): hole with vertical dimer rows; (e-f): island with vertical dimer rows; (g-h): island with horizontal dimer rows. The height of the terraces decreases from green to red to blue. The green and blue regions have vertical dimer rows, while the red ones have horizontal dimer rows.

8

 8

Figure 7 :

 7 Figure 7: (a-b) KMC snapshots of two successive positions of an island that moves under the eect of a force at θ=150 • . The black parallel lines show the direction of the dimer rows, the blue arrow show the direction of the electromigration force, forming an angle of θ with u. c: trajectories of the center of mass of islands moving because of a force with angle θ. Notice that the trajectories do not correspond to averages of several Monte Carlo simulations, so that the deviations from the straight lines depend on the pseudo-random numbers and have no physical meaning. d: logarithm of the simulated island velocity multiplied by kT as a function of 1/(kT).

Figure 8 :

 8 Figure 8: a: Black elongated holes with edges perturbed by electromigration. The grey regions are constituted by very small terraces.The black terrace at the left of the central grey region is lower than the bright terrace on the left. The step between them is some degrees away from a perfect S A step. The black terrace at the right of the grey central region is higher than the bright terrace on the right. The step between them is only some degrees away from a S B step. The height prole is on the top of the image. b: KMC simulation of step instabilities. The red terrace has dimer rows oriented horizontally and is one monolayer higher than the blue one. The electromigration force is oriented from left to right and has a value of 0.01, while all the other parameters have the values described in the text. The red terrace develops crests and valleys at the advancing edge. c: same conditions as in Figurea, but the diusion is isotropic (E d =0); in this case the edges do not move and do not develop crests and valleys.

Table 1 :

 1 Fluxes and velocities of holes and islands.

Table 3 :

 3 Values used to evaluate the eective charge.

Acknowledgments

We acknowledge the nancial support of the ANR grant HOLOLEEM (ANR-15-CE09-0012). We thank the CINaM-Planete facility and A. Ranguis for the technical support. We are grateful to A. Saul, B. Ranguelov and M. Michailov 10 CONCLUSIONS 10 for stimulating discussions we had in the framework of the PHC Rila project (38663TB).

Declarations of interest: none


J.-C Lizardo-Huerta

B Sirjean

L Verdier

R Fournet

P.-A Glaude

Pierre-Alexandre Glaude
email: pierre-alexandre.glaude@univ-lorraine.fr

Combustion and Pyrolysis Kinetics of Chloropicrin

Chloropicrin (CCl3NO2) is widely used in agriculture as a pesticide, weed-killer, fungicide or nematicide. It has also been used as chemical agent during World War I. The precise understanding of its combustion chemistry for destruction processes or in the event of accidental fire of stored reserves is a major safety issue. A detailed chemical kinetic model for the combustion and pyrolysis of chloropicrin is proposed for the first time. A large number of thermo-kinetic parameters were calculated using quantum chemistry and reaction rate theory.

The model was validated against experimental pyrolysis data available in the literature. It was shown that the degradation of chloropicrin is ruled by the breaking of the C-N bond followed by the oxidation of the trichloromethyl radical by NO2 through the formation of the adduct CCl3ONO, which can decompose to NO, chlorine atom and phosgene. Phosgene is much more stable than chloropicrin and its decomposition starts at much higher temperatures.

Combustion and pyrolysis simulations were also compared and demonstrated that the addition of oxygen has very little effect on the reactivity or product distribution due to the absence of hydrogen atoms in chloropicrin.

Introduction

Chloropicrin (CCl3NO2) is a nitroaliphatic compound, also known as trichloronitromethane. It is a colorless liquid which volatilizes slowly at room temperature and atmospheric pressure.

Chloropicrin has been widely used as herbicide, repellent, fumigant, insecticide, nematicide and disinfectant 1,[START_REF] Raman | Encyclopedia of Toxicology[END_REF] . It was also found to be an efficient chemical agent during World War, as a tear and a vomiting gas, also known as acquinite, G-25, S-1 or PS [START_REF] Pesonen | Chloropicrin-Induced Toxic Responses in Human Lung Epithelial Cells[END_REF][START_REF] Gupta | Handbook of Toxicology of Chemical Warfare Agents[END_REF] . Despite its little prevalence in environment, its high toxicity and environmental concerns have led to the prohibition of agricultural uses in Europe, but chloropicrin is still involved in many treatments in other countries, including the US for soil fumigation [START_REF] Triky-Dotan | Fumigant Dosages below Maximum Label Rate Control Some Soilborne Pathogens[END_REF] . The evaluation of danger in case of accidental fire and the destruction by thermal treatments of excess quantities of pesticides as well as stockpiles or unexploded ammunitions still found every year on old battlefields, requires the development of highly safe processes [START_REF] Bunnett | Arsenic and Old Mustard: Chemical Problems in the Destruction of Old Arsenical and Mustard Munitions[END_REF] . The yields of toxic emissions from incinerations or other thermal destruction processes are mostly controlled by chemical kinetics.

Few quantitative experimental studies have been devoted to chloropicrin thermal decomposition [START_REF] Gardner | [END_REF][START_REF] Radulescu | Sur La Décomposition Thermique de La Chloropicrine[END_REF][START_REF] Steacie | The Unimolecular Decomposition of Gaseous Chloropicrin[END_REF][START_REF] Dubikhin | Thermal Decomposition of Chloropicrin[END_REF][START_REF] Battin-Leclerc | Thermal Decomposition of Chloropicrin, Diphosgene and Phosgene between 100 and 530°C[END_REF] . High toxicity of this compound induces a high level of complexity for any experiments. Consequently, there is a considerable lack of knowledge on the behavior of this chemical under high-temperature conditions with or without oxygen. Chloropicrin is an active chemical species, which manifests low bond dissociation energies (BDE) compared to hydrocarbons. Our calculations at the CBS-QB3 level of theory show that the BDE of the C-Cl and C-N bonds are 68.9 and 42.1 kcal mol -1 respectively (Fig. 1). Note that chloropicrin, which does not contain any H-atom, is known as non-combustible, but enhances the reactivity of combustible mixtures and broaden the flammability limits of gas mixtures thanks to the oxidant role of the NO2 function [START_REF] Ashmore | A Study of Sensitized Explosions. VIII. Experimental Work on the Hydrogen-Oxygen Reaction Sensitized by Chloropicrin[END_REF][START_REF] Ashmore | A Study of Sensitized Explosions. X. The Kinetics of Decomposition of Chloropicrin and of the Hydrogen-Oxygen and Hydrogen-Chlorine Reactions Sensitized by Chloropicrin[END_REF] . The first analytical attempts demonstrated that the decomposition of chloropicrin starts at its boiling temperature (~112°C) and that the main products of decomposition are phosgene (COCl2) and nitrosyl chloride (ClNO) [START_REF] Gardner | [END_REF] . Radulescu and Zemfirecu 8 later measured the global rate of decomposition in the temperature range 130-160°C and found an activation energy around 40 kcal mol -1 , which is consistent with a reaction occurring mainly by the unimolecular breaking of the C-N bond, which BDE equals 42.1 kcal mol -1 (Fig. 1). In 1938, Steacie and Smith 9 pyrolyzed chloropicrin between 138 and 170°C. Their study confirmed that the reaction occurs by a unimolecular decomposition producing phosgene and nitrosyl chloride. The authors obtained a high pressure limit rate constant of 4.90×10 15 exp - 18958/T s -1 . They also observed that, under these temperature conditions, ClNO decomposed slowly to NO and Cl and concluded that the main products of chloropicrin decomposition are COCl2 and an equilibrium mixture of ClNO, NO and Cl2. Experimental work by Dubikhin et al. [START_REF] Dubikhin | Thermal Decomposition of Chloropicrin[END_REF] in 1971 on the same temperature range led to rate constant of 2.24×10 15 exp -18837/T s -1 , in very good agreement with previous data [START_REF] Radulescu | Sur La Décomposition Thermique de La Chloropicrine[END_REF][START_REF] Steacie | The Unimolecular Decomposition of Gaseous Chloropicrin[END_REF] . The authors observed that the reaction depends on pressure below 300 Torr. A more recent experimental work was carried out by Battin-Leclerc et al. [START_REF] Battin-Leclerc | Thermal Decomposition of Chloropicrin, Diphosgene and Phosgene between 100 and 530°C[END_REF] . In this work, the decomposition of chloropicrin under pyrolysis conditions was studied for temperatures ranging from 100 to 530°C in an isothermal closed reactor at an initial pressure of 25 Torr (0.033 bar) and for reaction times ranging between 10 and 120 min.

According to this work, phosgene is the main carbon containing reaction product; the other quantified products are majorly NO and Cl2 and the less important are CO and CCl4.

These studies show potential conditions for the destruction of chloropicrin, however to our knowledge there is no detailed kinetic model for the description of chloropicrin chemistry.

Following a methodology used for other chemical agents such as yperite (sulfur mustard) and HN-3 (nitrogen mustard) [START_REF] Sirjean | Kinetic Modeling of the Thermal Destruction of Mustard Gas[END_REF][START_REF] Lizardo-Huerta | Kinetic Modeling of the Thermal Destruction of Nitrogen Mustard Gas[END_REF] and more recently phosgene and diphosgene [START_REF] Lizardo-Huerta | Thermal Decomposition of Phosgene and Diphosgene[END_REF] , a comprehensive detailed kinetic model for the thermal decomposition and the combustion of chloropicrin has been developed in this work and validated against available experimental data.

Computational methods

Calculations from quantum chemistry were used to evaluate the high-pressure kinetic constants and thermodynamic data (∆fH°298 K, S°298 K, C°P(T)) associated with elementary processes involved in the primary mechanism of chloropicrin thermal decomposition, that are not available in the literature. Note that barrierless reactions, such as unimolecular initiations and radical combinations, which do not involve a defined transition state, were not treated theoretically. The calculations of the electronic structure of the reactants, products and transition states have been computed at the CBS-QB3 level of theory. This composite method involves a geometry optimization at the B3LYP/6-311G(2d,d,p) level of calculation and energy computations using several levels of theory. A basis set extrapolation allows to estimate a "complete basis set" (CBS) energy [START_REF] Montgomery | A Complete Basis Set Model Chemistry. VI. Use of Density Functional Geometries and Frequencies[END_REF][START_REF] Montgomery | A Complete Basis Set Model Chemistry. VII. Use of the Minimum Population Localization Method[END_REF] . In addition, an analysis of vibrational frequencies was systematically performed to confirm the nature of the transition states (TS).

In the cases where the imaginary frequency is not non-ambiguously related to the reaction coordinate, the Intrinsic Reaction Coordinate (IRC) approach [START_REF] Malick | Transition States for Chemical Reactions I. Geometry and Classical Barrier Height[END_REF] has been used to ensure the correct connection between a given transition state structure and the reactant and product. In this study, all the calculations have been performed using the GAUSSIAN 09 Rev. B.01 software 20 . For all the studied molecules, the enthalpies of formation have been evaluated using atomization reaction energies [START_REF] Petersson | Calibration and Comparison of the Gaussian-2, Complete Basis Set, and Density Functional Methods for Computational Thermochemistry[END_REF] . The references for atomization energies of atoms were taken from the CODATA database [START_REF] Cox | CODATA Key Values for Thermodynamics[END_REF] . Spin-orbit corrections have been taken into account [START_REF] Johnson | NIST Computational Chemistry Comparison and Benchmark Database NIST Standard Reference Database Number 101[END_REF] .

This methodology has proven its precision with molecules containing hetero-atoms [START_REF] Lizardo-Huerta | Kinetic Modeling of the Thermal Destruction of Nitrogen Mustard Gas[END_REF][START_REF] Lizardo-Huerta | Thermal Decomposition of Phosgene and Diphosgene[END_REF][START_REF] Khalfa | Quantum Chemical Study of the Thermochemical Properties of Organophosphorous Compounds[END_REF] . In addition, the low frequency vibration modes, corresponding to internal rotations, were treated as hindered rotors instead of harmonic oscillators, using THERMROT [START_REF] Lizardo-Huerta | Intramolecular Effects on the Kinetics of Unimolecular Reactions of β-HOROO˙ and HOQ˙OOH Radicals[END_REF] . This in-house code allows a one-dimensional hindered rotor (1D-HR) treatment of internal rotations and avoids the coupling of torsional modes encountered in the frequency analysis made by Gassian09, using the method proposed by Vansteenkiste et al. [START_REF] Vansteenkiste | An Extended Hindered-Rotor Model with Incorporation of Coriolis and Vibrational-Rotational Coupling for Calculating Partition Functions and Derived Quantities[END_REF] . 1D-HR treatment was used with torsional potential calculated from relaxed scans computed at the B3LYP/6-311G(2d,d,p) level of theory. Rate constant were eventually calculated as a function of temperature thanks to the Transition State Theory. The high pressure limit kinetic rate constants were fitted using the three parameter modified Arrhenius equation (A×T n ×exp(-Ea/RT)) over the temperature range 300-2000 K.

Chemical kinetic model development

The detailed kinetic model for the thermal decomposition of chloropicrin is composed of three sub-mechanisms: a first part corresponding to the decomposition of nitrogen species, another one dedicated to chlorinated species, and a mechanism specific to the decomposition and oxidation of chloropicrin and derived products.

Reactions of nitrogen species are taken from the mechanism proposed by Konnov et al. [START_REF] Konnov | The Pseudo-Catalytic Promotion of Nitric Oxide Oxidation by Ethane at Low Temperatures[END_REF] .

Since chloropicrin just involves C, N, and Cl atoms, the species containing hydrogen atoms in their chemical structure have been removed. In the same way, the sub-mechanism for chlorinated species has been adapted from the mechanism of Leylegian et al. [START_REF] Leylegian | Experiments and Numerical Simulation on the Laminar Flame Speeds of Dichloromethane and Trichloromethane[END_REF] for the reactions of chloromethanes where hydrogen containing species have been removed.

Recently, a new mechanism for HCl combustion has been proposed by Pelucchi et al. [START_REF] Pelucchi | High-Temperature Chemistry of HCl and Cl2[END_REF] which notably updates the reactions of light chlorinated compounds and re-examines the mechanism of flame inhibition involving the HCl/Cl2 system. Thus, we added to the mechanism of Leylegian et al. all the missing reactions and updated some kinetic data of reactions that appeared in both mechanisms. Some reactions of polychloromethanes not included in the previous mechanisms were also taken into account from the kinetic model proposed by Drakon and Eremin [START_REF] Drakon | Promotion of Methane Ignition by the Fire Suppressants CCl4 and CF3H. Combust[END_REF] .

A mechanism specific to chloropicrin has been developed in this work, in which were included unimolecular initiations, chlorine-atom abstractions, decomposition of the radicals formed, chlorine-nitrogen coupling reactions, as well as some oxidation reactions not included in kinetic models of chlorinated and nitrogenous species. The thermochemical properties were taken from the literature in the case of light species involved in the nitrogen and chlorine sub-mechanisms [START_REF] Konnov | The Pseudo-Catalytic Promotion of Nitric Oxide Oxidation by Ethane at Low Temperatures[END_REF][START_REF] Leylegian | Experiments and Numerical Simulation on the Laminar Flame Speeds of Dichloromethane and Trichloromethane[END_REF][START_REF] Pelucchi | High-Temperature Chemistry of HCl and Cl2[END_REF] , and theoretically calculated for chloropicrin and related products. Main species theoretically studied are presented in Table 1. The kinetic constants were either taken from the literature or theoretically obtained using the CBS-QB3 level of theory or estimated by analogy with similar reactions of the chlorinated species. The details of this mechanism will be presented in the following.

Mechanism specific to chloropicrin

Possible unimolecular initiation reactions involve the rupture of the C-N and C-Cl bonds, with the first being the easiest (BDE 42.1 vs. 68.9 kcal mol -1 , Fig. 1). Reactions were written in the reverse direction as barrierless radical combinations (reactions (1) and (2) in Table 2).

The kinetic constant of the first reaction was estimated from the correlations involved in the EXGAS software [START_REF] Buda | Progress toward a Unified Detailed Kinetic Model for the Autoignition of Alkanes from C4 to C10 between 600 and 1200[END_REF] and that of the second one by analogy with the reaction CCl3+Cl 32 which produces a C-Cl bond with a close BDE. Decomposition reactions of the produced radicals were included, based on analogies with species featuring similar chemical structures, like CH3NO2 [START_REF] Zhu | CH3NO2 Decomposition/Isomerization Mechanism and Product Branching Ratios: An Ab Initio Chemical Kinetic Study[END_REF] , CHCl2NO2 34 and C2Cl3NO2 [START_REF] Liu | Reaction Mechanisms of C2Cl3 + NO2 via Nitro and Nitrite Adducts[END_REF] . The authors showed that NO2 and the carbon radicals (CH3, CHCl2 or C2Cl3) react by combination to form either the reactant (addition of the N side of NO2) or a R-ONO species (addition of the O side of NO2). These combinations are barrierless. To better understand this behavior, the potential energy surface for the primary decomposition of chloropicrin was calculated (Fig. 2).

Figure 2.Energy profile at 0 K of the primary decomposition of chloropicrin, computed at the CBS-QB3 level of theory (in kcal mol -1 ).

In Fig. 2, the combination of CCl3 with NO2on the O-NO side produced Cl3CONO, which lies about 9 kcal mol -1 below chloropicrin (CCl3NO2). The latter can directly interconvert into Cl3CONO (TS1 and reaction 3 in Table 2) but the energy barrier is high (67.3 kcal mol -1 ) and remains much above the exit channel yielding trichloromethyl and NO2. Meanwhile, Cl3CONO decomposes easily through a cyclic transition state (TS2 and reaction 4 in Table 2), in which a Cl-atom is transferred on the N-atom, while the O-N bond breaks to produce

ClNO and phosgene COCl2. The geometry of TS2 is displayed in Fig. 3. The energy barrier is only 11.1 kcal mol -1 at 0 K. To verify the connection between the reactant and the products, an IRC calculation at the B3LYP/6-311G(2d,d,p) level of calculation was performed.

Chloropicrin will then readily react by an initial C-N bond breaking and the combination of produced radicals will further produce Cl3CONO (reaction 5). The latter can decompose by Est. [START_REF] Buda | Progress toward a Unified Detailed Kinetic Model for the Autoignition of Alkanes from C4 to C10 between 600 and 1200[END_REF] (2) Cl2CNO2+Cl  CCl3NO2 3.5210 13 0.0 0.0 3.510 [START_REF] Ashmore | A Study of Sensitized Explosions. X. The Kinetics of Decomposition of Chloropicrin and of the Hydrogen-Oxygen and Hydrogen-Chlorine Reactions Sensitized by Chloropicrin[END_REF] Est. [START_REF] Kumaran | Ab Initio Calculations and Three Different Applications of Unimolecular Rate Theory for the Dissociations of CCl4, CFCl3, CF2Cl2, and CF3Cl[END_REF] (3) CCl3NO2  Cl3CONO 3.8610 [START_REF] Battin-Leclerc | Thermal Decomposition of Chloropicrin, Diphosgene and Phosgene between 100 and 530°C[END_REF] The chlorine atom abstraction reaction was also considered (reactions 8). The kinetic parameters were estimated by analogy with the Cl-atom abstraction from CCl4 calculated theoretically by Brudnik et al. [START_REF] Brudnik | Theoretical Study of the Kinetics of Chlorine Atom Abstraction from Chloromethanes by Atomic Chlorine[END_REF] . The values were recalculated from the high-pressure limit rate constant given in the paper. The radical Cl2CNO2 can decompose by a rearrangement reaction, in which an oxygen atom from the nitro group is transferred to the carbon atom (reaction 9). The C-N and N-O bonds break simultaneously forming COCl2 and NO as products (Fig. 4). The combination of Cl2CNO2 with NO2 to yield Cl2C(ONO)NO2 (reaction 10) and the subsequent decompositions were also written by analogy with Cl3CONO and the rate of the subsequent decompositions calculated theoretically (reactions 11-12). The decomposition reactions of the primary products were also included. Reactions of phosgene COCl2 were taken from a previous study [START_REF] Lizardo-Huerta | Thermal Decomposition of Phosgene and Diphosgene[END_REF] . Reactions of nitrosyl chloride ClNO, and reactions of coupling between light chlorine and nitrogen containing species were taken from the literature or by analogy with equivalent reactions [START_REF] Zhu | Ab Initio Studies of ClOx Reactions: Prediction of the Rate Constants of ClO+NO for the Forward and Reverse Processes[END_REF][START_REF] Zhu | Ab Initio Studies of ClOx Reactions: Prediction of the Rate Constants of ClO+NO2 for the Forward and Reverse Processes[END_REF][START_REF] Manion | NIST Chemical Kinetics Database, NIST Standard Reference Database 17[END_REF] .

Results and discussion

The detailed kinetic model for the thermal decomposition of chloropicrin is composed of 54 species and 207 reactions and is available as Supplemental material. The Senkin program of Chemkin II [START_REF] Kee | A Fortran Chemical Kinetics Package for the Analysis of Gas-Phase Chemical Kinetics[END_REF] was used to perform simulations of the decomposition of chloropicrin under pyrolysis conditions in a close vessel. The simulations were done for temperatures ranging from 100 to 500°C, the reactor volume was 350 cm 3 , the initial pressure was set at 25 Torr (0.033 bar) and the residence time up to 10 minutes. These conditions correspond to the experimental study on the thermal decomposition of chloropicrin carried out by Battin-Leclerc et al. [START_REF] Battin-Leclerc | Thermal Decomposition of Chloropicrin, Diphosgene and Phosgene between 100 and 530°C[END_REF] .

Figure 5 shows the total pressure profiles as a function of residence time at different temperatures. Simulations are in very good agreement with the experimental results, without any fitting of rate constants, considering the experimental uncertainties due to the initial injection of the liquid reactant [START_REF] Battin-Leclerc | Thermal Decomposition of Chloropicrin, Diphosgene and Phosgene between 100 and 530°C[END_REF] . The constant pressure observed at 100°C implies that liquid chloropicrin is quickly vaporized after its injection into the reactor (the saturating vapor pressure of chloropicrin at 20°C is 20.25 Torr) and that no further reaction occurs thereafter.

The increase of pressure with time at 200°C indicates the onset of a slow chloropicrin decomposition. At 400°C, the pressure increases almost instantly after injection up to 66 Torr experimentally indicating that the chloropicrin is completely decomposed, but the very little change of the pressure over the next 10 min denotes a little reactivity of first reaction products at this temperature. The further increase in pressure at 500°C implies the occurrence of secondary decomposition reactions of the primary products. The model reproduces well these different behaviors at the different temperatures. Figure7 shows the evolution of the partial pressures of experimentally quantified products as a function of the temperature, for a residence time of 10 min. Results of simulations agree well for the main products, i.e. NO and phosgene. For CCl4, model predictions are correct for temperatures below 300°C, but underestimate them at higher temperatures. Note that wallpromoted combinations of CCl3 radicals with Cl-atoms cannot be disregarded. In the case of NO, the simulated profile displayed in Fig. 7 is the sum of the calculated mole fraction of NO and NOCl. This latter was indeed not detected in the experiments by Battin-Leclerc et al. [START_REF] Battin-Leclerc | Thermal Decomposition of Chloropicrin, Diphosgene and Phosgene between 100 and 530°C[END_REF] , whereas it was found to be the primary product of chloropicrin decomposition in previous experimental studies [START_REF] Gardner | [END_REF][START_REF] Steacie | The Unimolecular Decomposition of Gaseous Chloropicrin[END_REF][START_REF] Dubikhin | Thermal Decomposition of Chloropicrin[END_REF] . The very low Cl-NO BDE (37.2 kcal mol -1 ) makes almost impossible the quantification of the species by gas-phase chromatography, and we assumed that nitrosyl chloride decomposed completely into NO and Cl during the analytical process.

The model slightly overestimates the formation of NOCl at 225°C. However, for temperatures above 300°C a good agreement is obtained. With regard to phosgene, model predictions are in good agreement with the experimental results. This compound is a primary product of chloropicrin decomposition, which shows a much stronger thermal stability than the parent molecule. The decomposition of the phosgene starts to be noticeable around 370°C in these conditions [START_REF] Lizardo-Huerta | Thermal Decomposition of Phosgene and Diphosgene[END_REF] . In order to compare pyrolysis and combustion behavior of chloropicrin in incineration or accidental fire conditions, simulations were performed in a perfectly stirred reactor for mixtures of 5% of chloropicrin in air and in nitrogen under atmospheric pressure from 400 K to 1500 K for a residence time of 1 s. Chloropicrin profiles are identical in both conditions since reactant consumption occurs uniquely by the initial bond breaking rather than any radical bimolecular reaction. Product distribution and consumption are very similar. Figure 8 depicts the reaction pathways in pyrolysis and combustion at 1000 K. Chloropicrin reacts via the same reaction paths under both conditions sets by decomposition to NO2 and CCl3, followed by the combination to CCl3ONO and the formation of phosgene and nitrosyl chloride. The main difference appears in the products reactions. In these high temperature conditions, phosgene reacts by Cl-atom abstractions by NO and Cl in both conditions, but the role of NO decreases slightly in air compared to pyrolysis. In presence of O2, the NO2/NO ratio is indeed strongly modified by the oxidation of NO to NO2. Reaction pathways of ClNO are affected: reaction with Cl-atoms remains the main consumption channel, but reaction with NO2 accounts for 18% of ClNO consumption in air, but is almost negligible in pyrolysis.

Conclusions

The first detailed kinetic model of the thermal decomposition and the combustion of chloropicrin has been developed and validated. The mechanism is composed of submechanisms from the literature describing the reactions of the nitrogen containing and chlorinated species. The comprehensive mechanism for the decomposition of chloropicrin includes unimolecular reactions, Cl-atom abstractions, decomposition of the radicals formed, chlorine-nitrogen coupling reactions, as well as some oxidation reactions not included in the sub-mechanism of chlorinated and nitrogenous species. The kinetic constants were either taken from the literature or theoretically computed using the CBS-QB3 level of theory or evaluated by analogy with reactions included in the mechanism of the chlorinated species.

Simulations are in good agreement with experimental results from the literature obtained in pyrolysis in an isothermal closed reactor. Chloropicrin decomposes easily at low temperature to reach a complete decomposition in few minutes in a vessel at temperatures as low as 200°C. Its main products are phosgene, which is much more stable and reacts only at a higher temperature, above 370°C, and nitrosyl chloride, which decomposes to NO, Cl2 and Cl atoms.

Simulations predict that chloropicrin reacts identically in pyrolysis and combustion, because of the absence of H-atom and of little interaction of oxygen with decomposition intermediates. The effect of O2 is mainly a change in the NO2/NO ratio, which changes slightly reactions routes of phosgene and ClNO, but does not significantly change the overall kinetics.

Supplemental material

The detailed chemical kinetic model is available as supplemental material.

Figure 1 .

 1 Figure 1. Bond dissociation energies at 298 K in chloropicrin computed at the CBS-QB3

reaction 4

 4 through TS2, or by the O-NO bond breaking (bond dissociation energy 40.1 kcal mol -1 ) yielding CCl3O radicals and NO (reaction 6), which rate constant was set equals to that of CH3ONO[START_REF] Fernández-Ramos | Direct Dynamics Study of the Dissociation and Elimination Channels in the Thermal Decomposition of Methyl Nitrite[END_REF] . The rate parameter of the decomposition of CCl3O to phosgene and Cl-atom was calculated theoretically (reaction 7).

Figure 3 .

 3 Figure 3. CBS-QB3 geometry of chloropicrin transition state TS2. Bond length in Ǻ, angle in

Figure 4 .

 4 Figure 4. Reaction of rearrangement-elimination of radical CCl2NO2.

Figure 5 .

 5 Figure 5.Profiles of pressure vs. residence time for the thermal decomposition of

Figure 6 .

 6 Figure 6. Conversion of chloropicrin vs. temperature. Initial pressure, 25 Torr; residence

Figure 7 .

 7 Figure 7. Products quantified during the thermal decomposition of chloropicrin. Initial

Figure 8 .

 8 Figure 8. Flux analysis for the reaction of 5% of chloropicrin in nitrogen and in air at 1000K

Table 1 :

 1 Thermochemical properties of chloropicrin and related species determined at the CBS-QB3 level of theory.

		∆ 𝑓 𝐻 298 ° (kcal mol -1 )	𝑆 298 ° (calmol -1 K -1 )	𝐶 𝑝 °(𝑇)(calmol -1 K -1 ) 300 K 400 K 500 K 600 K 800 K 1000 K 1500 K
	CCl3NO2	-19.81	88.11	24.29	27.45	29.72	31.31	33.19 34.22	34.39
	Cl2CNO2	19.30	82.42	20.72	23.24	25.17	26.61	28.47 29.47	29.59
	Cl3CONO	-28.79	91.99	27.60	30.08	31.75	32.83	34.00 34.69	34.89
	Cl2C(ONO)NO2 -23.93	99.17	33.42	37.84	40.99	43.16	45.56 46.57	46.74
	CCl3O	-10.36	75.88	20.01	21.80	22.94	23.66	24.50 24.94	25.40
	CCl3OO	-7.60	84.70	23.35	25.77	27.34	28.31	29.21 29.63	29.79
	ClNO	12.88	62.64	10.62	11.20	11.65	11.99	12.47 12.82	12.88
	ClONO	14.32	68.43	14.29	15.82	17.05	18.01	19.35 20.10	20.24
	ClONO2	3.86	74.07	16.64	18.82	20.46	21.67	23.20 24.01	24.13
	COClNO2	-40.68	78.94	19.14	21.48	23.30	24.69	26.60 27.78	27.91

Table 2 .

 2 Reactions of decomposition of chloropicrin and its primary products. Units are cm 3 ,

	mol, s, cal with k = A T n exp(-E/RT).
	Reaction	A	n	E	k500K
	(1) CCl3+NO2  CCl3NO2	2.010 12	0.0	0.0	2.10 12

Acknowledgements

This work was supported by DGA Maîtrise NRBC. This work was granted access to the HPC resources of IDRIS under the allocation 2017-A0010807249 made by GENCI.