Thomas Jourdan
email: thomas.jourdan@cea.fr

Gabriel Stoltz

Frédéric Legoll

Laurent Monasse

An accurate scheme

Keywords: Rate equations, Cluster dynamics, Fokker-Planck equation, Ostwald ripening

come

Introduction

Nucleation, growth and coarsening of particles are three mechanisms which appear in gases, liquids and solids. For example, in solids the formation and evolution of secondary-phase precipitates or point defect clusters such as voids, dislocation loops and stacking fault tetrahedra can be described by a combination of these three mechanims. The evolution of a system in terms of particle size distributions thus can be modeled by coupling methods specifically designed to describe one particular mechanism [1,2]. An alternative modeling approach relies on the observation that the three processes result from the absorption and emission of mobile particles by other particles (birth and death process).

To avoid the approximations which are needed when methods are coupled, it is therefore interesting to directly model the evolution of particles due to these absorption and emission processes.

Only a few methods allow a seamless description of the nucleation, growth and coarsening processes by simulating the absorption and emission of particles over time scales comparable to experiments. Atomistic and Object kinetic Monte Carlo (kMC) methods are classical methods to perform such tasks [3,4,5]. Their computational cost however limits the size of the systems which can be simulated, so that it is impossible to simulate the whole evolution of a realistic microstructure. Cluster dynamics, or rate equation cluster dynamics (RECD), appears to be an efficient alternative which can accurately describe a locally homogeneous system with a moderate volume fraction of particles. In RECD, the system is modeled as a gas of particles (often called clusters), which can migrate in the solvent (the matrix in case of a solid), associate and dissociate [6,[START_REF] Clouet | Fundamentals of Modeling for Metals Processing[END_REF]. Comparisons with kMC methods have shown that, provided the parameters of the RECD models are carefully chosen, particle size distributions (also called cluster distributions in the following) obtained by the two methods are very similar [START_REF] Jourdan | [END_REF]9]. RECD models are written as a set of ordinary differential equations (ODEs) which correspond to the birth and death process. The variables are the con-centrations of each cluster type (see details below). Despite their simplicity, RECD equations are notably difficult to solve. They are in general very stiff, so that the integration in time should be performed with implicit methods [START_REF] Hairer | Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems[END_REF].

These methods require solving linear systems based on the Jacobian matrix of the ODE set. Since there is one equation per cluster size, the number of equations and thus the size of the Jacobian matrix can be larger than 10 6 to handle clusters of a few tens of nanometers. This can become prohibitively large in terms of computational time and memory use, even when the sparsity of the matrix is taken into account [START_REF] Jourdan | [END_REF]. Alternatively, stochastic methods can be used [12,13,14]: they do not require the storage of the Jacobian matrix but the time step is often much smaller than with deterministic methods.

In order to reduce the number of equations to solve, several methods have been introduced. They can be divided into two main classes, the "grouping" methods and the methods based on a partial differential equation (PDE) obtained in the continuum limit. The former approach, first developed by Kiritani [15], consists in solving equations for the moments of the distribution over groups of equations. It was shown by Golubov et al. [16,17] that using at least the zeroth and first moments of the distribution was necessary to obtain accurate results. The second class of methods is based on the observation that, for large cluster sizes, the master equation can be approximated by a Fokker-Planck equation. By discretizing this equation on a sufficiently coarse mesh, a significant reduction in the number of equations can be achieved. However, assessment of the accuracy of this class of methods remains scarce compared to grouping methods. It is sometimes even argued that Fokker-Planck based methods inherently cannot be accurate [START_REF] Golubov | Comprehensive Nuclear Materials[END_REF].

In this article we show that, by carefully discretizing the Fokker-Planck equation, it is possible to achieve high accuracy in the simulation of cluster distributions. A proper spatial discretization of the advection term turns out to be critical to obtain qualitatively and quantitatively accurate results. As an application, we simulate the coarsening of loops, for which an analytical form of the loop distribution is known under some approximations. These various approximations are discussed in light of our RECD calculations.

Our article is organized as follows. We recall the RECD model and its continuum Fokker-Planck limit in Section 2. Various discretization schemes for the Fokker-Planck equation are discussed in Section 3, with an emphasis on the spatial discretization of the advection term. Numerical results are then presented in Section 4. First, we compare in Section 4.1 the various numerical schemes on a simple test case where reference simulations can be performed, in order to validate the scheme presented in Section 3.3. We next use this scheme in Section 4.2 to assess the relevance of various modeling assumptions to describe loop coarsening in aluminum. Our conclusions are collected in Section 5.

Master equation and the Fokker-Planck approximation

Master equation approach to RECD

For the sake of simplicity, we consider in the following the evolution of cluster concentrations due to the absorption and emission of monomers. Clusters are assumed to contain only one type of species. For secondary-phase precipitates or point defect clusters in solids, this species represents a solute, a vacancy or a self-interstitial. Clusters thus can be referred to by a single index n = 1, . . . , N and their concentration is denoted by C n . More general cases (clusters with multiple species and absorption of mobile clusters) are described in [START_REF] Jourdan | [END_REF]. In the case considered here, the RECD equations to solve are

dC n dt = β n-1 C n-1 C 1 -α n C n -β n C n C 1 + α n+1 C n+1 for any n ≥ 2, ( 1
)
dC 1 dt = -2β 1 C 1 C 1 + 2α 2 C 2 - n≥2 [β n C n C 1 -α n+1 C n+1 ] . (2)
The coefficients β n and α n represent the absorption and emission rates of monomers respectively. For example, if clusters are spherical and if the reactions between clusters are controlled by diffusion, these coefficients read [START_REF] Waite | [END_REF]20]

β n = 4πr n D 1 , (3)
α n = β n-1 V at exp - F b n kT , (4)
where r n is a reaction distance, V at is the atomic volume, F b n is the binding free energy of a monomer to a cluster containing n -1 monomers, k is the Boltzmann constant, T is the temperature and D 1 is the diffusion coefficient of the monomer given by

D 1 = D 1,0 exp (-E m 1 /kT ). (5)
In this expression, D 1,0 is a diffusion prefactor and E m 1 is the migration energy of the monomer.

Using (1)-( 2), it is straightforward to check that the total quantity of matter

S is conserved, ie dS dt = d dt   n≥1 nC n   = 0. (6)

Fokker-Planck approach to RECD

Using a Taylor expansion, it can be shown [21,22] that Eq. (1) can be approximated, for n 1, by a Fokker-Planck equation of the form

∂ ∂t C(x, t) = - ∂ ∂x [F (x, t)C(x, t)] + ∂ 2 ∂x 2 [D(x, t)C(x, t)] , (7)
with

x n = n, C(x n , t) ≈ C n (t) and F (x, t) = β(x)C 1 (t) -α(x), (8)
D(x, t) = 1 2 [β(x)C 1 (t) + α(x)] . (9)
To perform this approximation, it has been assumed that there exist two functions α and β such that α n = α(x n ) and β n = β(x n ). In practice this is the case for large values of n, for which F b n and r n in Eqs. ( 3) and ( 4) can be expressed analytically as functions of n (an example is given in Section 4.2). The transformation to a Fokker-Planck equation is interesting as such, since the two physical ingredients present in the rate equations are clearly highlighted. The advection term, proportional to F (x, t), is related to cluster growth. The diffusive part, proportional to D(x, t), is related to cluster size fluctuations due to random absorption and emission of monomers. It is responsible for a broadening of the cluster distribution. If clusters contain multiple species, a multidimensional Fokker-Planck equation can be written in a very similar way [23]. The method proposed in this work for the advective term can simply be applied componentwise. For the diffusion part, cross-terms of the type ∂ 2 ./∂x∂y appear if clusters containing multiple species are mobile. The treatment of such terms is not addressed in the present article.

Since we intend to use a Finite Volume formulation in Section 3, it is convenient to rewrite Eq. ( 7) in a conservative form, as a function of the flux J. It also permits to exhibit the advection and diffusion fluxes J a and J d :

∂C ∂t = - ∂J ∂x = - ∂J a ∂x - ∂J d ∂x (10)
with

J a (x, t, C(x, t)) = F (x, t)C(x, t) (11)
J d (x, t, C(x, t)) = - ∂ ∂x [D(x, t)C(x, t)] . (12)
In general no analytical solution can be found for the partial differential equation [START_REF] Clouet | Fundamentals of Modeling for Metals Processing[END_REF], but it can be discretized along x and solved numerically as an ODE set. The advantage with respect to (1)-( 2) is that a coarse mesh size ∆x may often be used. The resulting set of ODEs is thus posed in a space of lower dimension than (1)- (2). Usually, due to the strong stiffness of the ODE set, implicit methods with backward differentiation formula or Runge-Kutta time stepping are used [START_REF] Hairer | Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems[END_REF].

Discretizing and directly solving the Fokker-Planck equation for cluster dynamics is not a very common approach. Among the few attempts that can be found in the literature, one can note the works by Ghoniem and Sharafat [24,25] and, more recently, Turkin and Bakai [26]. In these works, the advective term is taken care of by centered finite differences. This approach has been shown to give satisfactory results when the advective flux is small, such as for the simulation of aging [27,[START_REF] Jourdan | [END_REF]. However, for large advection terms, this scheme is known to be unstable [START_REF] Fletcher | Computational Techniques for Fluid Dynamics 1: Fundamental and General Techniques[END_REF], which manifests itself through spurious oscillations in the cluster distributions as will be shown in Sec. 4.1. To solve this problem, upwind schemes [START_REF] Fletcher | Computational Techniques for Fluid Dynamics 1: Fundamental and General Techniques[END_REF][START_REF] Ozkan | [END_REF] or related schemes such as the Chang-Cooper scheme [30] can be used [31,[START_REF] Jourdan | [END_REF]. Essentially, the Chang-Cooper scheme automatically and seamlessly switches from a centered scheme to an upwind scheme as the advective flux increases, so that the resulting scheme remains stable. The drawback of this method is that it introduces (as any upwind scheme) some numerical diffusion [START_REF] Fletcher | Computational Techniques for Fluid Dynamics 1: Fundamental and General Techniques[END_REF], which can alter the cluster distribution by enhancing the effect of the diffusive term D(x, t) [START_REF] Ozkan | [END_REF][START_REF] Jourdan | [END_REF]. If one is interested only in the zeroth and first moments of the distribution (corresponding to the density and average size of clusters respectively) or if the distribution is not strongly peaked, the Chang-Cooper scheme can be used and yields accurate results [START_REF] Jourdan | [END_REF]. In other cases, this scheme should be used with caution.

Coupling the master equation and the Fokker-Planck limit

Another aspect to consider while using a discretization of the Fokker-Planck equation ( 7) is that it is valid only for n 1. For small values of n, only the master equation ( 1)-( 2) can accurately describe the evolution of clusters. This is especially true given that binding energies F b n (and so the parameters α n through (4)) often strongly vary for small cluster sizes, and that cluster distributions are known to be very sensitive to these parameters [32]. This observation led to the development of hybrid schemes [24,26,[START_REF] Jourdan | [END_REF], which use the master equation for small clusters (in the so-called "discrete region") and the Fokker-Planck equation for larger clusters ("continuous region"). A key problem is the coupling between the discrete and continuous regions. For this purpose, it is interesting to note that in the continuous region, a centered scheme with a mesh size equal to 1 reduces to the original master equation when only monomers are mobile [26]. By imposing a centered scheme and a mesh size equal to one at the boundary between the discrete and continuous regions, it is thus possible to seamlessly couple the ODE (1)-( 2) to the PDE [START_REF] Clouet | Fundamentals of Modeling for Metals Processing[END_REF]. If small clusters are mobile, correction terms should be added to obtain a rigorous coupling which ensures the continuity of cluster fluxes at the interface between the two regions. It has been shown that, with this method, the quantity of matter S is conserved up to machine precision [START_REF] Jourdan | [END_REF]. This type of coupling is therefore used in all calculations performed in the present work.

Discretization of the Fokker-Planck equation

To discretize the Fokker-Planck equation, we decouple time and space. Equations are first discretized in space, and the system of ODEs is next solved using an appropriate ODE solver. Indeed, RECD equations are in general very stiff, so we want to take advantage of the implicit solvers [33,34] which have proved to be efficient to solve such equations [START_REF] Jourdan | [END_REF]. In the following, we use the CVODES package [33], which contains variable-order and variable-step linear multistep methods with backward differentiation formulas. As discussed in the next sections, most spatial discretizations of the advection term are built specifically to ensure total variation diminishing (TVD) properties for the forward Euler scheme, or for any Runge-Kutta scheme which can be decomposed into successive explicit Euler stages. The TVD property ensures that an initially monotone solution remains monotone at future time steps, and does not become oscillatory. Although it cannot be guaranteed that CVODES is as precise and efficient as TVD solvers for the sole advection term, it was shown to give very satisfactory results with much lower computation times than TVD solvers for all the tests we performed.

Let us now focus on the spatial discretization of the Fokker-Planck equation.

We consider the non-uniform mesh of Fig. 1, where C n represents the average of C(x, t) over the cell of width ∆x n (for the sake of clarity, we omit the time dependence in the following). To conserve the flux of clusters, we discretize the conservative form [START_REF] Hairer | Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems[END_REF]:

dC n dt = - J n+1/2 -J n-1/2 ∆x n , (13)
where J n+1/2 is the flux of clusters between classes n and n + 1. It is the sum of the diffusion term J d n+1/2 and the advective term J a n+1/2 .

x x n-5/2 x n-3/2 x n-1/2 x n+1/2 x n+3/2 x n+5/2 x n-2 x n-1 x n x n+1 x n+2 C n-2 C n-1 C n C n+1 C n+2 ∆x n Figure 1: One-dimensional non-uniform mesh
The discretization of the diffusion term ( 12) can be done with the classical centered scheme:

J d n+1/2 = D n+1 C n+1 -D n C n ∆xn 2 + ∆xn+1 2 , (14)
with

D n = D(x n , t).
As mentioned previously, the discretization of the advection flux is more difficult. Various numerical schemes offer a better accuracy than the upwind scheme while avoiding oscillations produced by the centered scheme. Among these schemes, which are generally called "high-resolution schemes" [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF], we only consider the semi-discrete ones. In addition, although J a depends linearly on C(x, t) in the present case, it can be nonlinear in general, since the drift F can also depend on cluster concentrations [START_REF] Brailsford | [END_REF]. It is therefore important to use a scheme that can also be extended to nonlinear equations.

One particularly interesting semi-discrete scheme adapted for nonlinear equations is the central-upwind scheme developed by Kurganov, Noelle and Petrova (KNP) [37]. In the wake of the scheme proposed by Kurganov and Tadmor for nonlinear conservation laws [38], it has the advantage of being a Riemannsolver-free approach, thus leading to simple implementation. With the KNP scheme, the advective flux is given by

J a n+1/2 = a + n+1/2 J a x n+1/2 , C - n+1/2 -a - n+1/2 J a x n+1/2 , C + n+1/2 a + n+1/2 -a - n+1/2 + a + n+1/2 a - n+1/2 a + n+1/2 -a - n+1/2 C + n+1/2 -C - n+1/2 . ( 15
)
In this expression, C - n+1/2 and C + n+1/2 are the left and right values of the reconstruction of C at point x n+1/2 from the values {C j } j∈N . Possible definitions for C - n+1/2 and C + n+1/2 will be detailed in Sections 3.1, 3.2 and 3.3 below. The reconstruction procedure is at the heart of the numerical method.

Before moving on to that reconstruction procedure, we explain how a + n+1/2 and a - n+1/2 are defined. In the current case, C is a scalar. In the general case where C ∈ R N , we denote by λ -

1 < • • • < λ - N the N eigenvalues of the Jacobian matrix ∂J a /∂C(C - n+1/2 ) and λ + 1 < • • • < λ + N the N eigenvalues of the Jacobian matrix ∂J a /∂C(C + n+1/2 ). Then the coefficients a + n+1/2 and a - n+1/2 read a + n+1/2 = max λ - N , λ + N , 0 , (16)
a - n+1/2 = min λ - 1 , λ + 1 , 0 . (17)
It should be noted that for the (linear) problem under consideration, it holds

a + n+1/2 = max F (x n+1/2 ), 0 , (18)
a - n+1/2 = min F (x n+1/2 ), 0 , (19)
so for a positive advective term (F (x) > 0), the flux simply is

J a n+1/2 = J a x n+1/2 , C - n+1/2 = F (x n+1/2 )C - n+1/2 . ( 20
)
Similarly, for a negative advective term, the flux reduces to

J a n+1/2 = J a x n+1/2 , C + n+1/2 = F (x n+1/2 )C + n+1/2 . (21)
In the linear case, this expression is the natural upwind flux expression. Different options are available to compute reconstructed left and right values C - n+1/2

and C + n+1/2 . The goal of reconstruction is to increase the order of the method in smooth regions while retaining stability, by avoiding oscillations near discontinuities or strong gradients. Various reconstructions are presented in the following sections.

Minmod reconstruction

The minmod reconstruction is a simple reconstruction, which guarantees that the KNP scheme is total variation diminishing (TVD) when the forward Euler scheme is used to discretize the left-hand side of (13), see [38]. As a consequence, an initially monotone solution remains monotone at future time steps, and does not become oscillatory.

For a non-uniform mesh, the minmod reconstruction can be written as (see [START_REF] Leveque | Finite-Volume Methods for Hyperbolic Problems[END_REF])

C + n+1/2 = C n+1 - ∆x n+1 2 σ n , (22)
C - n+1/2 = C n + ∆x n 2 σ n , (23)
where σ n is a so-called slope limiter. It is given here by

σ n = minmod C n -C n-1 ∆xn 2 + ∆xn-1 2 , C n+1 -C n ∆xn+1 2 + ∆xn 2 , (24)
where, for M arguments, the minmod operator reads

minmod (r 1 , . . . , r M ) =            min{r 1 , . . . , r M } if r i > 0 for all i = 1, . . . , M , max{r 1 , . . . , r M } if r i < 0 for all i = 1, . . . , M , 0 otherwise. ( 25
)
This simple reconstruction performs in general much better than the upwind scheme (which corresponds to σ n = 0 everywhere), but it reduces the slope rather strongly and therefore tends to smear out solutions. Near extrema, as all TVD schemes, it reduces to a first order scheme (upwind scheme) and the solution is in general smoothed.

Monotonized central-difference reconstruction

The monotonized central-difference (MC) reconstruction follows the slope limiter approach. It is given by Eqs. ( 22)-( 23) where σ n is now given by

σ n = minmod C n+1 -C n-1 ∆xn-1 2 + ∆x n + ∆xn+1 2 , 2 C n -C n-1 ∆xn 2 + ∆xn-1 2 , 2 C n+1 -C n ∆xn+1 2 + ∆xn 2 . ( 26
)
As the minmod reconstruction, it can be shown that MC reconstruction is TVD when the explicit Euler scheme is used for the time integration [38], so there can be loss of accuracy near extrema. However, since the slopes are less reduced than with the minmod reconstruction, it leads in general to more accurate results.

MP5 reconstruction

To avoid smoothing of the solution at extrema, it is necessary to relax the TVD constraints in these regions. This is the case for the MP5 reconstruction of Suresh and Huynh [39], which is monotonicity preserving (MP). In particular initial upper and lower bounds on the solution are preserved at later times. The value C + n+1/2 is defined using the five values C n-2 , C n-1 , C n , C n+1 and C n+2 (Fig. 1). Essentially, the use of five points instead of three points as in the two previous methods makes it possible to distinguish a discontinuity (where slope limiting should be active) from an extremum (where slope limiting leads to excessive smearing and should not be used). Though not TVD, MP5 reconstruction does not produce oscillatory solutions when a forward Euler scheme is used to discretize the left-hand side of (13). Any Runge-Kutta method made of successive applications of explicit Euler integrations in time may also be used.

In [39], this scheme was derived for a uniform mesh. It can however be generalized to non-uniform meshes. We have built a non-uniform variant of the scheme (see Appendix A) and only found improvements in terms of accuracy when using meshes exhibiting very strong variations. For the spatial meshes considered in Sec. 4 to discretize the Fokker-Planck equation, no difference is visible between the results obtained by the two variants of the scheme. It should also be noted that the non-uniform scheme is slightly more complicated in its formulation and slightly less efficient. Therefore, in the following, we use the scheme as derived for uniform meshes although our meshes are not uniform.

Results and discussion

We consider two examples in this section. The first one is the evolution of quenched-in vacancies (see Section 4.1), which provides an interesting setup to test the accuracy of the numerical schemes. In particular, it allows us to validate the quality of the MP5 reconstruction. We next use in Section 4.2 the scheme with MP5 reconstruction to test the sensitivity of cluster distributions to the absorption and emission terms β n and α n . More specifically, many experimental studies have focused on the cluster distributions during coarsening by Ostwald ripening [START_REF] Ardell | Phase transformations '87[END_REF], following the seminal theoretical work by Lifshitz and Slyozov [START_REF] Lifshitz | [END_REF] and Wagner [42] (LSW) for three-dimensional particles and by Kirchner [43] and Burton and Speight [44] (KBS) for two-dimensional dislocation loops. We study here the influence of some parameters on the cluster distribution shapes during the coarsening of dislocation loops, a simulation which is impossible by a direct integration of the original ODE set of RECD.

Clustering of quenched-in vacancies

We consider Eqs. (1)-( 2) with the following parameters:

β n = β for n ≥ 2, (27)
β 1 = ηβ, (28)
α n = 0 for n ≥ 2, ( 29
)
for some η < 1. The initial condition corresponds to a material with quenched-in vacancies: C 1 (t = 0) = C q and C n (t = 0) = 0 for any n ≥ 2.

This test case was first proposed by Koiwa [45] to investigate the validity of the grouping method of Kiritani [15], and was subsequently used to validate the grouping method of Golubov et al. (G-method) [16]. Since clusters cannot dissociate (α n = 0), the ratio of the advective flux to the diffusive flux is high (see Eqs. ( 8) and ( 9)), which makes this test particularly interesting to validate the numerical schemes for the Fokker-Planck equation, and in particular to investigate numerical diffusion. As time increases, the concentration of monomers C 1 decreases due to the absorption reactions, until it becomes zero. At this time, the cluster distribution does not evolve any more.

To perform this test, we simply choose β = a 0 D 1 , where a 0 is the lattice parameter and D 1 is the diffusion coefficient (Eq. ( 5 Table 1: Parameters used for the simulation of clustering of quenched-in vacancies [16].

In Fig. 2 we show the stationary cluster distributions obtained for η = 10 -4

with different numerical schemes: the Chang-Cooper scheme [30], which reduces to the upwind scheme in this case (F n = F (x n )):

J a n+1/2 = F n C n , (30)
the centered scheme: As expected, the upwind method (Fig. 2-a) introduces substantial numerical diffusion and the distribution broadens. It should be noted, however, that the zeroth and first moments of the distribution are well reproduced. The error on the cluster density (n ≥ 2) is about 5 %, while for the average number of vacancies inside clusters it is only 2 %. Owing to its simplicity and high efficiency (see below), this method is attractive if the quantity of interest is the cluster density or the average number of vacancies. However, if a precise description of the cluster distribution is needed, this scheme should be avoided.

J a n+1/2 = 1 2 (F n C n + F n+1 C n+1 ) , (31)
The centered scheme (Fig. 2-b) introduces spurious oscillations in the solution.

It is therefore clearly inappropriate for this test, and more generally for all cases with strong advection terms. Significant improvement is obtained with the simple minmod reconstruction (Fig. 2-d): no oscillations appear and it produces less numerical diffusion than the upwind scheme. Numerical diffusion can be even more reduced by the use of MC reconstruction (Fig. 2-e), although near the peak one can still observe a small deviation with respect to the reference calculation. Among the reconstructions that were tested, MP5 clearly appears to be the best one; it perfectly reproduces the reference calculation (Fig. 2-f).

Incidentally, this calculation also shows that the error due to the transformation of the master equation into a Fokker-Planck equation is completely negligible.

Although the Fokker-Planck approximation is derived from a Taylor expansion of the master equation, which is strictly valid only for n 1, we observe numerically that using this approximation for n as low as 5 does not introduce error in the solution. Concerning grouping methods, the G-method (Fig. 2-c) performs very well also, the only difference with the MP5 reconstruction being the presence of unphysical negative concentrations in the large cluster size tail of the distribution (value for x n around 180 in Fig. 2-c). Such unphysical effects cannot appear with the MP5 reconstruction since the scheme is monotonicity preserving. We conclude from this test that the MP5 scheme is indeed very accurate.

The primary purpose of this work is to compare the accuracy of the discretization schemes. It is however instructive to compare the computation times required to perform the simulations with a given time solver. The wall-clock computation times on a single core processor are given in Tab. 2 using CVODES package. To make the comparison sensible between the different schemes, the Jacobian matrix is evaluated by CVODES using finite differences and the standard dense linear solver included in the package is used. In all cases, the computation time is rather short to reach t = 2 × 10 5 s (less than 10 s), but it varies from one method to another. The shortest times are for the Chang-Cooper and centered schemes. It gradually increases with the complexity of limiters. For larger systems, linear algebra dominates the overall computation time and one can expect the efficiency of the grouping method to decrease, since it uses two equations per group, so approximately twice more equations than the Fokker-Planck based methods. A systematic study of the efficiency is however difficult to perform, since a time solver can be more or less adapted to a spatial discretization. For example, implicit methods as implemented in CVODES are probably not optimal when limiters are used, due to the frequent switches in the discretization stencils. Implicit-explicit methods, treating advective fluxes in the explicit part, may be more efficient. We emphasize however that even with CVODES, the computation time is reasonable with the MP5 scheme. If a large number of calculations are required, as for parameter fitting, it may be advantageous to use Chang-Cooper scheme or minmod limiter to get rough solutions, and then refine the final solution with the MP5 scheme. for the test proposed by Koiwa [45]. Time integration is performed on a 2.80 GHz processor using CVODES, with finite difference calculation of the Jacobian matrix and the dense linear solver included in the package.

Loop coarsening

Loop coarsening by Ostwald ripening was studied theoretically by Burton and Speight [44] and Kirchner [43] (a model known as KBS). They found that the loop distribution in the coarsening regime obeys the following law:

C(ρ, t) = A ρ (2 -ρ) 4 exp - 4 2 -ρ 1 r(t) 2 , (32)
where ρ = r/ r(t) , r is the loop radius and r(t) is the mean loop radius of the cluster distribution. The parameter A is a normalizing constant, which can be expressed, for example, as a function of the total concentration of defects inside loops, which remains constant during the coarsening process. For a pure prismatic loop, the loop radius r is related to the number of particles x by πr 2 b = xV at , where b is the Burgers vector.

To make the problem tractable, in addition to the classical hypotheses of the LSW theory for 3D particles [START_REF] Lifshitz | [END_REF]42], two main approximations are made in the KBS approach in order to obtain Eq. ( 32):

1. The formation energy of unfaulted, circular loops is given by a simple line tension model:

E f (r) = 2πrT , (33)
where T is the line tension, which is considered constant. Another expression is given by Hirth and Lothe [START_REF] Hirth | Theory of dislocations[END_REF] using isotropic elasticity, for an isolated pure prismatic loop:

E f (r) = 2πr µb 2 4π(1 -ν) ln 4r r c -1 , ( 34
)
where µ is the shear modulus, ν is Poisson's ratio and r c is a core cut-off radius. Since 4), the logarithmic dependency may affect the emission rates α n and therefore the shape of the loop distribution.

F b n = E f 1 + E f n-1 -E f n with E f n = E f (r n ) in Eq. (
2. The diffusion flux of defects to a loop can be written as the flux to a spherical cluster of the same radius, which means that β n is given by Eq. (3) with πr 2 n b = x n V at (the loop is assumed to be pure prismatic). Therefore, the toroidal shape of the cluster is not taken into account.

Seeger and Gösele [START_REF] Seeger | [END_REF] have shown that, when the diffusion problem is solved in toroidal geometry, the absorption coefficient reads

β n = 2πr n 2π ln 8rn rp D 1 , (35)
where r p is a core cut-off radius. More accurate formulas take into account the drift-diffusion process due to the elastic interaction between the defect and the loop [48,49].

To test the validity of the formula (32) given by the KBS theory and the sensitivity of cluster distributions to absorption rates β n and to the energetics of loops, RECD calculations were performed with different approximations:

(a) Spherical model for the absorption kinetics (Eq. ( 3)), line tension model for the energetics of loops (Eq. ( 33)), as for the KBS approach.

(b) Spherical model for the absorption kinetics (Eq. ( 3)), formula based on isotropic elasticity for the energetics of loops (Eq. ( 34)).

(c) Toroidal model for the absorption kinetics (Eq. ( 35)), line tension model for the energetics of loops (Eq. ( 33)).

The list of parameters corresponding to loops in aluminum is given in Tab. 3.

In our calculations the line tension is given by T = κµb 2 , where κ has been adjusted to give a continuous evolution for binding energies from di-vacancies, the value of which is imposed, to larger clusters, for which the line tension is used. Likewise, the value of the core cut-off radius r c in Tab. 3 ensures a smooth transition between di-vacancies and larger loops.

Initial distributions are only populated for the class n = 100 (x n ≈ 391 due to the geometric progression of the mesh), with a concentration equal to With parametrization (a), which is based on the same approximations as KBS theory, a good agreement is obtained (Fig. 3-(a)). The profile given by RECD is slightly less peaked and broader than the theoretical profile. We have checked, by varying the mesh size, that numerical diffusion is negligible and cannot explain this small discrepancy. Most probably it comes from the probabilistic nature of the growth law in RECD, whereas in KBS it is purely deterministic. In RECD, growth occurs by successive absorptions of monomers but emissions can also occur, whereas in KBS theory any cluster of size larger than the critical radius will grow. Such an argument has been put forward to explain and given by the KBS theory (gray line).

the discrepancy between LSW theory and RECD calculations of the coarsening of 2D deposits [52]. Broader distributions were also obtained with RECD, but the difference was much more important between LSW and RECD than in the present case. Our results are also in line with the work by Senkov [53], who pointed out the role of fluctuations to explain why experimental distributions are often broader than those predicted by LSW theory.

When Eq. ( 34) based on isotropic elasticity is used, results are almost identical to those obtained with the simple line tension model (Fig. 3-(b)). This result is compatible with dislocation dynamics simulations of the coarsening of loops [START_REF] Bakó | [END_REF]: the authors obtained distributions in good agreement with KBS theory, although no assumption on the energetics of loops as in Eqs. ( 33) and ( 34) is made in this type of simulations. These results can also be compared to experimental results in silicon. Some authors found that the loop distributions followed a KBS law [54], while others found loop distributions with much longer tails [55]. They ascribed this discrepancy to the logarithmic term in Eq. (34).

Although our simulations are performed on a different material, they clearly do not support this interpretation.

Interestingly, a more significant difference appears between RECD results and KBS theory if the absorption coefficient takes into account the toroidal geometry of the loop (Fig. 3-(c)). In this case, the distribution is sharper.

This can be explained by the logarithmic term in Eq. ( 35) which decreases the absorption rate for large clusters with respect to the spherical approximation, thereby lowering their growth rate.

It should be noted that the influence of other loops on the absorption flux due to overlapping diffusion fields and on the energetics of loops is not taken into account in these calculations. More complex models would be necessary to include these effects, see e.g. [56], but this is beyond the scope of the present work. What we want to emphasize here is the importance of using a proper numerical scheme to obtain physical distributions which are not distorted by numerical diffusion, in order to compare the various physical models. Essentially, our conclusions would have been difficult, if not impossible to draw, with other numerical schemes.

Conclusion

In this work we have introduced a numerical scheme to accurately and efficiently solve cluster dynamics equations. It is based on a coupling between the master equation (for small clusters) and the discretization of the Fokker-Planck equation (for large clusters). To avoid the distortion of cluster distributions in the region where the Fokker-Planck equation is used, special care must be taken.

We have shown that high resolution methods must be used for the advection term. More specifically, the KNP scheme coupled with the MP5 reconstruction appears to be a very attractive method. On the cases that have been tested, the results obtained by MP5 reconstruction are in very good agreement with the reference solution even when using rather large mesh sizes. The so-obtained numerical method compares favourably with the grouping methods, since it uses twice less equations and prevents the appearance of negative concentrations. In addition, the use of the Fokker-Planck equation permits to clearly identify the drift and diffusion terms, which can be interesting for physical interpretation.

The method can be easily generalized to simulations containing two species 1. A high order approximation of C - n+1/2 is calculated, based on the values {C j } j∈N ; 2. The value of C - n+1/2 is possibly changed (limiting process) to ensure that if the {C j } j∈N are initially monotone, they remain monotone at the subsequent timestep. This amounts to imposing the value of C - n+1/2 to lie inside some prescribed intervals which depend on the initial values of {C j } j∈N .

Limiting leads to a lower order reconstruction, but guarantees that no oscillations appear. Intervals can be set to impose TVD property, but it can lead to too much smearing near extrema. With MP5 reconstruction the intervals are enlarged near extrema, to keep as much as possible the initial high order approximation of C - n+1/2 : although TVD property is lost in general, monotonicity is preserved and the accuracy is higher than with TVD schemes.

at m - 3 η

 3 )). Numerical values of parameters are given in Tab. 1. The value of β only has an influence on the time at which the distribution reaches its final stationary state[45]. This time is set to 2 × 10 5 s in our calculations. Coefficient for the clustering of two vacancies 10 -4

Figure 2 :

 2 Figure 2: Long time limit cluster size distributions for the test proposed by Koiwa [45] (clustering of quenched-in vacancies) using different approximations (circles): (a) Chang-Cooper (upwind scheme), (b) centered scheme, (c) grouping method (G-method), KNP with (d) minmod, (e) MC and (f) MP5 reconstructions. The exact calculation, taking into account all cluster classes, is shown by the black solid line.

  classes. Calculations are performed with MP5 reconstruction up to a time at which the normalized distributions are stationary. It has been checked that it is the case for all parametrizations at t = 5 s. In Fig.3 (a)-(c) we show the normalized cluster distributions C(ρ, t) = C(ρ, t) r(t) 2 /A obtained by RECD within the Fokker-Planck approach at t = 5 s as well as the one predicted by the KBS theory.

Figure 3 :

 3 Figure 3: Distributions obtained by RECD using various approximations (black line, see text)

Table 2 :

 2 Wall-clock computation times to reach t = 2 × 10 5 s with various numerical schemes,

	Scheme	Computation time (s)
	Chang-Cooper	1.3
	Centered	1.3
	Grouping	3.0
	KNP + minmod	2.8
	KNP + MC	4.2
	KNP + MP5	9.1

Table 3 :

 3 10 23 m -3 . This initial condition permits to determine the normalizing constant Parameters used for the simulation of vacancy loop coarsening in aluminum.

	A in (32). The use of a Fokker-Planck approach enables us to simulate cluster
	distributions with less than 1000 classes, whereas a simulation with all cluster
	classes would require approximately 16 million equations. The first 20 cluster
	classes are described by the master equation, then a geometric progression is
	used with = 0.03 for 100 classes and is decreased to 0.01 for the remaining

by applying the algorithm along each dimension, except if cross-terms in the Fokker-Planck equation are present [START_REF] Jourdan | [END_REF]. Such terms arise if clusters containing several species are mobile. For higher dimensions, stochastic methods are probably more appropriate due to the increase in the number of equations to solve [12,13,14]. It should also be noted that the method is not adapted if large clusters are mobile since it is based on a Taylor expansion, whereas grouping methods are applicable to this case [57].

Using this method, we have simulated the loop coarsening in aluminum and compared our results to the KBS theory. When the same assumptions are taken for the absorption coefficients and energetics of loops, a good agreement is obtained. Cluster distributions are found to be only slightly broader in RECD calculations, which can be ascribed to the probabilistic nature of the growth law in RECD. Additional calculations show that normalized cluster distributions do not depend much on the energetics of loops, but that they are sensitive to the model used for the absorption coefficients.

Acknowledgments

This work was performed within the NEEDS project MathDef (Projet Blanc).

We gratefully acknowledge the financial support of NEEDS. The CRESCENDO code, which was used for this study, is developed in collaboration with EDF R&D within the framework of I3P institute.

Appendix A. MP5 reconstruction

We recall here briefly the MP5 reconstruction by Suresh and Huynh [39], adapted to the case of non-uniform meshes. We assume that F (x) is positive at x n+1/2 , so in view of (20) we need to determine C - n+1/2 . The arguments can be easily adapted to the case F (x) < 0.

The method proceeds in two steps that are successively detailed in Appendix A.1 and Appendix A.2 below:

. Original interface value

The value at the interface C - n+1/2 is defined by C - n+1/2 = Cn (x n+1/2 ), where

Following Ref. [58], the coefficients a i are obtained by imposing that

where I n is the interval corresponding to index n. For a uniform mesh, this leads to

Higher order reconstructions may be used [39], but in practice for cluster dynamics problems this reconstruction appears accurate enough. To determine C - n+1/2 in the non-uniform case, a linear system must be solved at each interface point for the non-uniform scheme, which makes it more complicated to implement and less efficient than the uniform scheme. For meshes with low spatial variation as those considered in the present work, we have found that the two versions lead to nearly identical results.

Appendix A.2. Slope limiter

The following steps are performed to preserve the monotonicity of the solution:

• The following values are computed:

with γ ≥ 2. We adopted γ = 4, in agreement with Ref. [39].

is left unchanged. Otherwise, we must perform the limiting process which is detailed in the next step.

with:

(A.8)

(A.9)

For a uniform mesh and γ = 2, inequality (A.13) leads to ζ n ≤ 9/2, which is in line with the choice ζ n = 4 in Ref. [39] (we recall that we have set γ = 4 in this work). Reducing ζ n activates limiting more often, thereby increasing stability at the cost of accuracy.

Finally, the following timestep restriction is necessary to ensure that the scheme is monotonicity preserving with an explicit Euler integration: ∆t ≤ min


Laurence Conty
email: laurence.conty@ens.fr

David Gimmig

Clément Belletier

Nathalie George

Pascal Huguet

The cost of being watched: Stroop interference increases under concomitant eye contact

Keywords: Eye contact effect, Stroop interference, Social perception, Visual cognition

Current models in social neuroscience advance that eye contact may automatically recruit cognitive resources. Here, we directly tested this hypothesis by evaluating the distracting strength of eye contact on concurrent visual processing in the well-known Stroop's paradigm. As expected, participants showed stronger Stroop interference under concomitant eye contact as compared to closed eyes. Two control experiments allowed ruling out low-level account of this effect as well as non-specific effect of the presence of open eyes. This suggests that refraining from processing eye contact is actually as difficult as refraining from word reading in the Stroop task. Crucially, the eye contact effect was obtained while gaze was not under the direct focus of attention and the participants were faced with another powerful distracter (the incongruent word) in the task at hand. Thus, there is a cost of being watched even in circumstances where the processing of direct gaze is strongly disfavored. The present results emphasize the crucial status of eye contact in human cognition.

The cost of being watched:

Stroop interference increases under concomitant eye contact

Eye contact is one of the first and most important visual cues processed by human infants with normal development [START_REF] Batki | Is there an innate gaze module? Evidence from human neonates[END_REF][START_REF] Farroni | Eye contact detection in humans from birth[END_REF]. Its processing appears currently fundamental for building critical inferences about another's level of interest and state of mind [START_REF] Baron-Cohen | Mindblindness: An essay on Autism and Theory of Mind[END_REF][START_REF] Kleinke | Gaze and eye contact: a research review[END_REF][START_REF] Patterson | Non Verbal Behavior: A functional perspective[END_REF], but also for acquiring self-other consciousness [START_REF] Vasudevi | On being the object of attention: implications for self-other consciousness[END_REF]. For these reasons, eye contact perception is considered as one of the basic building blocks for the development of the so-called "social brain" [START_REF] Adolphs | The social brain: neural basis of social knowledge[END_REF], a network of structures involved in human inter-individual interaction. In line with this view, several authors proposed that eye contact may be processed by a specialized brain module [START_REF] Batki | Is there an innate gaze module? Evidence from human neonates[END_REF][START_REF] Perrett | Understanding the intentions of others from visual signals: Neuropsychological evidence[END_REF][START_REF] Senju | The eye contact effect: mechanisms and development[END_REF]. Furthermore, according to the 'fast-track modulator' model [START_REF] Senju | The eye contact effect: mechanisms and development[END_REF], a brain network involved in eye contact processing i) allows detecting direct gaze automatically in the visual environment, and ii) subsequently induces modulation of concurrent cognitive processing and behavioral response (hereafter referred to as the "eye contact effect" or ECE).

Consistent with Senju and Johnson's view of the ECE, human observers are faster at detecting a face or eyes with direct gaze than at detecting averted gaze in visual search tasks [START_REF] Conty | Searching for asymmetries in the detection of gaze contact versus averted gaze under different head views: a behavioural study[END_REF][START_REF] Senju | Does perceived direct gaze boost detection in adults and children with and without autism? The stare-in-the-crowd effect revisited[END_REF][START_REF] Grünau | The detection of gaze direction: A stare-in-the-crowd effect[END_REF]. Direct gaze also facilitates other face-related tasks such as gender discrimination [START_REF] Macrae | Are you looking at me? Eye gaze and person perception[END_REF] as well as the encoding of individual faces and recognition of face identity [START_REF] George | Seen gaze-direction modulates fusiform activity and its coupling with other brain areas during face processing[END_REF][START_REF] Mason | Look into my eyes: gaze direction and person memory[END_REF]. However, most empirical support for the ECE comes from experiments with tasks focused on face perception and/or recognition. This may have encouraged gaze processing, as the eyes are a key facial feature [START_REF] Yarbus | Eye Movements and Vision[END_REF].

The idea of a brain module involved in eye contact processing should imply more largely that the mere presence of direct gaze automatically calls for processing resources [START_REF] Fodor | The modularity of mind[END_REF], at the expense of any concurrent visual processing outside the facial area. In line with this argument, Senju and Hasegawa (2005) found evidence that peripheral target detection was delayed when the target was preceded by a centrally fixated face gazing directly at the subject, as compared with a face gazing downward or with closed eyes (the detection of a peripheral target asterisk was slowed down at least when the presence of direct gaze overlapped in time with the target). This finding suggests how difficult the attentional disengagement from eye contact may be, with sometimes negative effects on concurrent cognitive processing outside the facial area.

However, in Senju and Hasegawa's experiment, the face was the central stimulus on which the participants had to focus before performing the peripheral-target detection task. Again, this may have favored the observation of an ECE. Here, we make a step further by showing that eye contact is a so powerful distracter that it can burden processing resources even when gaze stimulus is not initially presented under the direct focus of attention and when participants are faced with another powerful distracter on the task at hand: the incongruent word in Stroop's Paradigm.

The Stroop task requires individuals to identify the color of the ink in which a letter string is printed. Typically, the time needed to respond is greater when the string represents a word incongruent with the ink color (e.g., the word RED printed in green), as compared to neutral signs (e.g., XXX printed in green), a robust finding called Stroop interference [START_REF] Stroop | Studies of interference in serial-verbal reaction[END_REF]. A core assumption of virtually every theoretical account of this interference is that skilled readers process automatically the incongruent word (MacLeod, 1991), which is therefore conceived as a powerful distracter. If eye contact is itself a powerful distracter, one may expect that the mere presence of concomitant on-screen direct gaze during the Stroop task enhances processing during word reading and/or imposes further delay to the inhibition of the predominant colour word response. This should result in higher Stroop interference under concomitant direct gaze, as compared for example to closed eyes. Such finding would constitute the first direct behavioral evidence for the distracting strength of eye contact.

Experiment 1

Method

Participants. Thirty four undergraduates (mean age=21.2±0.5 years) volunteered to participate, all reported normal or corrected-to-normal vision and were naïve about the purpose of the experiment (presented as part of a larger research project on colour perception).

Stimuli. Letter string stimuli in the incongruent condition were four colour words (RED, BLUE, GREEN, YELLOW) presented in the three colours that did not match their meaning (e.g. RED printed in blue, green and yellow). In the neutral condition, four strings of neutral signs matched for grapheme-length with the four colour words were used (e.g. XXX for RED). Each letter string subtended 1° of visual angle horizontally and 0.5° vertically. Eye stimuli consisted of the eye region cut from 24 face photographs [START_REF] Conty | Searching for asymmetries in the detection of gaze contact versus averted gaze under different head views: a behavioural study[END_REF]. For each original face, there were four different stimuli in which the eye condition (gazing directly at the observer or closed) and head position (frontal or rotated sideways by 30°) were varied in a 2 × 2 factorial design (Figure 1A). Eye stimuli were mirror imaged resulting in a total of 192 pairs of eyes. All eye stimuli were in grey scale and subtended 6-7° of visual angle horizontally and 1.6° vertically.

Procedure. The participants arrived individually and were met by either a male or a female experimenter (random assignment) who described the Stroop task. Each participant began by a practice period of 128 to 224 trials during which he/she learned to associate the ink colour (red, green, blue, yellow) of neutral words (e.g. WIND) to the appropriate keyboard response (S, F, K, M). Then, the participant was trained to the classical Stroop task with the incongruent color words and the neutral strings presented randomly and in equal proportion over 96 trials. Finally, in the test period of 192 trials, the participants performed the Stroop task but with the eye stimuli now appearing 0.5° of visual angle above the string. The association between eye and string stimuli was counterbalanced within subjects. These stimulus pairs were presented in random order. At the beginning of the test period, participants were told that each string would appear concomitantly with opened or closed eyes and that the eyes were not informative for the task.

They were requested to ignore the presence of the eyes and avoid reading the strings. Thus, they were explicitly focused on the ink color. In every trial, the stimulus display was preceded by a fixation cross centred on the letter string (see Figure 1B). Statistical Analyses. The percentage of correct responses was very high in all conditions (mean = 96.1±0.4%) and thus was not analysed. Outlier reaction times (RT) exceeding two standard deviations above the mean were rejected (<5% per subject and string condition). Only RT data from correct responses were analyzed. We checked that these data showed normal distribution in every condition according to the Kolmogorov-Smirnov test (all p.>15). Thus, a three-way ANCOVA was carried out on RTs of the test period with Gaze Condition (direct/closed), Head Orientation (frontal/deviated), and Type of string (incongruent/neutral) as within-subjects factors, and the interference score (RT difference for incongruent minus neutral strings) of the training period as a covariate (F (1,32) =12.5, p<.001). Removing the covariate did not change the findings. Effect size (η 2 ) are reported together with F and p values for the main effects and interactions. Planned comparisons were performed for the analysis of simple main effects when interactions were observed.

Results and Discussion

Not surprisingly, we found the classical Stroop interference. RTs were longer for the incongruent words than for the neutral signs (F (1,32) = 82.8; p <.001; η 2 =.72). More importantly, the type of eyes appearing above the string clearly modulated this interference (F (1,32) =7.4, p=.01; η 2 =.18): the Stroop effect was stronger in the context of direct gaze (M =87±11 ms; F (1,32) = 85.1, 9 p < .0001) than in the context of closed eyes (M =53±11 ms; F (1,32) = 24.7, p <0.001) (Figure 2A). This influence of eye gaze on Stroop interference was independent of head orientation (F<1, non significant, ns). Furthermore, the modulation of the Stroop interference by gaze was uniquely due to the effect of the Gaze Condition on incongruent trials (planned comparison, F (1,32) =5.8; p<.03) (Figure 2B). No other effects were found (all Fs<1, ns). Despite the image for direct gaze was very different in frontal and deviated views of a face, the present ECE occurred for both head orientations. This suggests that the ECE could not be driven simply by low-level image properties. Nevertheless, whatever head orientation, direct gaze and closed eyes physically differed. This difference is inherent to the strong local contrast between the iris and the sclera characterizing opened eyes [START_REF] Kobayashi | Unique morphology of the human eye[END_REF], which makes these stimulus visually more complex and possibly more salient than closed eyes. Thus, we ran a control experiment where eye stimuli were replaced by control images preserving the mean luminance of the original photographs while showing clear different local contrast and complexity.

Experiment 2

Method

Participants, stimuli, and procedure. A new sample of 34 participants (12 females, mean age =19.8±0.4 years) was exposed to the same procedure as before. However, direct eye gaze stimuli were replaced by vertical gratings of the same size and global contrast as the corresponding eyes, and closed eyes were replaced by uniform gray images. All control images were matched for mean gray level with the original eye stimuli (Figure 1C).

Statistical analyses.

As before, the percentage of correct responses was very high in all conditions (mean = 96.1±0.4%). Outlier RTs were rejected (<5% per subject and string condition). These data showed normal distribution in every condition according to the Kolmogorov-Smirnov test (all p.>15). A two-way ANCOVA was carried out on RTs of the test period with Grating Condition (grating/uniform gray), and Type of string (incongruent/neutral) as within-subjects factors, and the interference score (RT for incongruent minus RT for neutral strings) of the training period as a covariate. Removing the covariate did not change the findings.

Results and Discussion

As before, the participants showed the classical Stroop interference (F (1,32) =73.6, p<.001; η 2 =.69). However, the type of image appearing above the strings had absolutely no effect on the size of the Stroop effect (F<1, ns -mean interference effect for gratings=85±10ms; mean interference effect for uniform grey image=89±13ms). Accordingly, Grating condition modulated the RT neither for the incongruent words nor for the neutral strings (all F<.1, nssee Table 1). This control experiment strengthened our confidence that the enhanced Stroop interference under concomitant eye contact (Experiment 1) was not simply due to low-level visual properties of direct gaze compared to closed eyes. However, at this stage, it cannot be rejected that this result reflected the impact of any gaze direction (i.e. the mere presence of open eyes) rather than a specific effect of direct gaze. We addressed this important issue in a third experiment where direct gaze stimuli were replaced by averted gaze stimuli.

Experiment 3 Method

Participants, stimuli, and procedure. A new sample of 34 participants (21 females, mean age =20.0±0.2 years) was exposed to the same procedure as in Experiment 1. However, direct gaze stimuli were replaced by the corresponding averted gaze configurations. Thus, for each original face, there were four different stimuli in which the eye condition (gazing sideway from the observer or closed) and head position (frontal or rotated sideways by 30°) were varied in a 2 × 2 factorial design (see Figure 1D).

Statistical analyses.

As before, the percentage of correct responses was very high in all conditions (mean = 95.8±0.4%). Outlier RTs were rejected (<5% per subject and string condition). These data showed normal distribution in every condition according to the Kolmogorov-Smirnov test (all p.>20). A three-way ANCOVA was carried out on RTs of the test period with Gaze Condition (averted/closed), Head Orientation (frontal/deviated), and Type of string (incongruent/neutral) as within-subjects factors, and the interference score (RT for incongruent minus RT for neutral strings) of the training period as a covariate. Removing the covariate did not change the findings.

Results and Discussion

As before, the participants showed the classical Stroop interference (F (1,32) =56.7, p<.001; η 2 =.64). However, the type of eyes appearing above the strings had absolutely no effect on the size of the Stroop effect (F<1, ns -mean interference effect for averted gaze=76±12ms; mean interference effect for closed eyes=76±13ms). Accordingly, Gaze condition modulated the RT neither for the incongruent words (F (1,32) =1.1; p=0.30) nor for the neutral strings (F (1,32) =2.7; p=0.11 -see Table 1). This last result strengthened our confidence that the higher Stroop interference under concomitant direct gaze relative to closed eyes (Experiment 1) was specifically related to an eye contact effect.

General Discussion

Preliminary findings in the ECE literature have indicated that, when attention is initially engaged on faces, direct gaze (i.e., eye contact) can hold attention at the expense of concurrent visual processing outside the facial area (Senju & Hasegawa, 2005). Here, we make a step further by directly showing the distracting strength of eye contact in circumstances where its processing was strongly disfavored. The ECE was indeed obtained despite the gaze stimulus was not under the direct focus of attention and the participants were faced with another powerful distracter (the incongruent word) in the task at hand. The instructions discouraged eye processing as well.

Consistent with our expectation that eye contact may interfere above and beyond word reading in Stroop's paradigm, participants displayed significantly more color-word interference in the presence of direct gaze, as compared to closed eyes. Hence, processing eye contact seems actually as difficult to refrain as word reading in the Stroop task. There is ample evidence that word reading is automatic in that readers cannot refrain from accessing word meaning (MacLeod, 1991). If refraining from processing eye contact was not at least as difficult, the presence of direct gaze in the vicinity of the incongruent words would not have influenced the size of the Stroop effect. Instead, this effect was stronger under concomitant eye contact. As indicated by experiment 2, this finding can hardly be explained by differences in local contrast characterizing direct gaze relative to closed eyes. Moreover, the size of the Stroop effect did not change under averted gaze (Experiment 3), indicating that the stronger Stroop interference found under direct gaze in experiment 1 was specifically related to eye contact perception. Thus, our findings can be taken as evidence that eye contact is processed even when initially presented outside the direct attentional focus, which also reinforces the view that this contact involves specific mechanisms of detection [START_REF] Conty | When eye creates the contact! ERP evidence for early dissociation between direct and averted gaze motion processing[END_REF][START_REF] Perrett | Understanding the intentions of others from visual signals: Neuropsychological evidence[END_REF].

The lack of averted gaze effect on the Stroop interference is consistent with [START_REF] Burton | Gaze perception requires focused attention: evidence from an interference task[END_REF] who concluded that averted gaze does not interfere with performance when occurring outside the focus of attention. In their research, participants made directional (left-right) judgments to gazing-face or pointing-hand targets, which were accompanied by a distracter face (only averted gaze was used) or hand. Averted gaze actually never interfered with performance even when distracter sizes were increased to compensate for their peripheral presentation, or when making use of only the eye regions as distracters thereby eliminating a potential difficulty in extracting gaze from a face context. Of particular interest here, our own findings in Stroop's paradigm rather support a different conclusion regarding direct gaze, which seems to be an extremely powerful distracter. Future research might clarify whether this distracter does or does not operate in Burton's own interference paradigm. Importantly, in these authors' paradigm, eye gaze was presented in short display times of 200 ms. In our own research, the stimuli were on display until a response was made, and thus it is possible that the eye stimuli were processed serially (being out the primary focus of attention). Future research therefore might also specify whether the present ECE in Stroop's paradigm can be found under the more drastic condition of short display times.

Of particular interest here, the ECE was actually restricted to the (most demanding) incongruent word condition. This supports the view that eye contact did not interfere with the selection of the letter color cues per se but rather competed with and/or delayed the inhibition of word reading. It may therefore be the case that eye contact consumed common processing resources shared with the inhibition of predominant colour word processing and/or response.

Alternatively, eye contact may have selectively enhanced the processing of the words as suggested by recent studies emphasising the influence of direct gaze on the quality of information processing [START_REF] Fullwood | Effect of gazing at the camera during a video link on recall[END_REF][START_REF] Yoon | Communication-induced memory biases in preverbal infants[END_REF]. Interestingly, the present ECE fits with recent findings (Conty et al., in press) indicating that direct gaze may automatically and specifically increase arousal during the performance of demanding tasks.

Future research shall clarify whether increased Stroop effect in presence of direct gaze is mediated by increased arousal. Furthermore, the absence of ECE in the neutral string condition contrasted with Lavie et al.'s findings (2003) indicating that visual distracters, especially faces, produce interference in visual search task independently of the task perceptual load. However, in Lavie et al.'s (2003) research, the targets (persons' names) were always semantically linked to the distracters (faces). Here, the targets (string color) were completely unrelated to the distracters (gaze). This major difference may help explain the inconsistency between Lavie et al.'s (2003) findings and our own results. It seems that when there is no relationship between the targets and the social distracters, the task must be difficult enough (incongruent string condition as opposed to neutral string condition) for an interference to occur (for a similar argument about direct gaze and arousal, see Conty et al., in press). This point will deserve special attention in future research on the ECE.

In any case, our results show for the first time that direct gaze selectively enhances Stroop interference. Complementary with Senju and Hasegawa's (2005) findings, the present results

indicate how distracting eye contact can be, at least when this contact occurs concomitantly with the task at hand. Perhaps more importantly, our results offer first evidence that eye contact can be distracting even when direct gaze is not presented under the direct focus of attention and competes with another powerful distracter. This is the key contribution of the present research. It reinforces the view that eye contact processing may be mandatory, in agreement with the idea that it is adaptive not to ignore faces looking at us [START_REF] Emery | The eyes have it: the neuroethology, function and evolution of social gaze[END_REF]. Moreover, our findings extend previous research demonstrating that the processing of social stimuli (such as faces) can have automatic discrete influences on our behavioural responses [START_REF] Bateson | Cues of being watched enhance cooperation in a real-world setting[END_REF][START_REF] Csibra | Natural pedagogy[END_REF][START_REF] Doherty-Sneddon | Cognitive demands of face monitoring: evidence for visuospatial overload[END_REF][START_REF] Fergusson | How social perception can automatically influence behavior[END_REF].

In sum, the present results provide the first direct behavioral evidence for the distracting strength of eye contact in circumstances where its processing was strongly disfavored by the experimental context. They bring support to the intuitive feeling that there is a cost of being watched. This cost may be related to the interpretation of direct gaze as an evaluative [START_REF] Droney | Attributions of self-esteem as a function of duration of eye contact[END_REF] or threatening [START_REF] Grünau | The detection of gaze direction: A stare-in-the-crowd effect[END_REF] social signal, or to the adaptive value of paying attention to faces looking at us for the purpose of social interaction [START_REF] Emery | The eyes have it: the neuroethology, function and evolution of social gaze[END_REF]. It would be interesting to test whether the eye contact effect on Stroop interference is or is not observed in individuals suffering from severe social and communicative deficits such as autism, who typically display deficient attention to others' eyes and gaze (Senju et al., 2005;[START_REF] Spezio | Abnormal use of facial information in high-functioning autism[END_REF]. Paradoxically, such individuals may show a benefit of their indifference

Figure 1 :

 1 Figure 1: A) Example of eye stimuli under the 4 experimental conditions in Experiment 1. Direct

Figure 2 :

 2 Figure 2: The Eye Contact Effect (ECE) on Stroop interference. Data were averaged over Head

Table 1 .

 1 Mean reaction times with standard error of the mean obtained in the main experimental conditions of control experiments 2 and 3.

		Experiment 2	Experiment 3
		Grey Image	Grating	Closed eye	Averted gaze
	Incongruent string	742±28	745±27	720±26	730±27
	Neutral string	655±20	653±24	643±18	653±18

to the eyes (relative to control subjects) in the present paradigm. As for now, the present findings emphasize how much eye contact constitutes a crucial visual cue in the cognitive functioning of healthy individuals.