Magali Benoit

Jonathan Amodeo

Ségolène Combettes

Ibrahim Khaled

Aurélien Roux

Julien Lam

Measuring transferability issues in machine-learning force fields: The

I. INTRODUCTION

Atomistic modeling is often divided in two different types of simulations. On the one hand, quantum methods including Hartree-Fock and DFT approaches are considered the most accurate and are employed for virtually any types of chemical species [START_REF] Marx | Ab Initio Molecular Dynamics : Basic Theory and Advanced Methods[END_REF][START_REF] Martin | Electronic structure : basic theory and practical methods[END_REF]. On the other hand, classical force fields are used to perform large-scale and long-time simulations with less accuracy [START_REF] Frenkel | Understanding Molecular Simulation[END_REF][START_REF] Stone | The theory of intermolecular forces[END_REF]. However, it is still difficult to connect both approaches and until now, one can hardly perform a simulation involving millions of atoms for nanoseconds while retaining the accuracy of quantum methods.

In this context, machine-learning interaction potentials (MLIP) have been proposed in the recent years and have shown great potentials to achieve such simulations [START_REF] Behler | [END_REF][6][7]. Numerous approaches are currently considered including Artificial Neural Networks [8], Gaussian approximation methods [9], Linearized potentials [10,11], Spectral Neighbor Analysis Potential [12], Symmetric Gradient Domain Machine learning [13,14] and Moment Tensor Potentials [15]. The success of these techniques is recognized by the large variety of materials that were successfully tackled: pure metals [16][17][18][19][20], organic molecules [21][22][23][24], oxides [25,26], water [27][28][29][30][31], amorphous materials [32][33][34][35][36][37] and hybrid perovskites [38].

For all of these techniques, the main procedure consists in using a very universal analytical formulation for the force field which is then parameterized to match a database of DFT calculations including total energy, forces and stress tensors. However, it is admitted that MLIP can sometimes show poor transferability towards systems that are not included in the learning database. In the worst scenario, the MLIP is so-well fitted to its learning database that non-physical behaviors may be observed outside of it. In order to fix this issue, the main proposal is to regularly check the accuracy of the potential as the machine-learning molecular dynamics simulations are carried out and to improve the MLIP "on the fly" [38][39][40]. Yet, to the best of our knowledge, such flaw of the approach has never been quantitatively investigated while being acknowledged by both users and developers.

For our case study, we choose interactions in gold-iron nanoparticles. In principle, such system can be found concurently in three different chemical orderings namely alloy, Janus and core-shell. Yet, recent experiments have shown that the synthesized Au-Fe nanoparticles are made of an iron core wrapped in a gold shell and that the shape of the iron core depends strongly on the amount of surrounding gold [41][42][43][44][45]. These nanoparticles have potential biomedical applications as iron is known for its intrinsic ferromagnetism and gold capping can protect the iron core from oxidation. However, rationalizing the results of the synthesis along with predicting the material properties would require numerical simulations which are sparse for gold-iron nanoparticles [45][46][47][48][49]. Indeed, while full quantum calculations can not be employed to study clusters of more than tens of atoms, the empirical modeling of gold-iron nanoparticles is a also very difficult case because these two metals are non miscible at room temperature on a very large domain of the phase diagram. There are therefore no iron-gold alloy crystal structures on which to adjust the parameters of an empirical potential. Moreover, iron is magnetic and crystallizes in a bcc structure while gold crystallizes in an fcc structure, which makes the development of a potential capable of capturing all the properties of this alloy even more complex. Previous attempts have shown their limits by stabilizing metastable alloys [50] or by not being able to find the most stable Fe/Au interfaces [46,47], leading us to develop potentials specifically dedicated to a particular problem, and hence highly non-transferable [49].

In this article, we begin by describing the methodology including linearized machine-learning potential and a penalizing regression scheme. In the results section, we first studied the influence of the descriptor functions. Then, we showed that the methodology allows one to quickly obtain MLIPs with different degrees of complexity. Afterwards, the transferability of these different potentials was tested on forces in untrained chemical orderings namely Janus and Core-shell. While the error should decrease monotonically when increasing the MLIP complexity, we observed a surprising non-monotonic behavior thus illustrating that more complexity does not necessarily lead to a better MLIP overall. Such transferability issue was reduced by using a more diverse set of descriptors. Finally, we measured some properties of the bulk and investigate the possibilities and the limitations of the obtained MLIP towards bulk simulations even if it was trained on nanostructures.

II. METHODS

A. The Φ-Lassolars machine-learning interaction potential

For our MLIP, we employed the analytical formulation originally put forward by Seko et al. [10,17,18,51,52]. In this method, the total potential energy of a configuration made of N atomic positions is first given by E tot = N i=0 E i where E i is the atomic energy. For E i , we considered a weighted linear combination of descriptors indexed by n:

E (i) = n ω n X (i) n (1)
where ω n is the linear coefficient associated with the descriptor X (i) n . Until now, moment tensors [16], grouptheoretical high-order rotational invariants [52] and bispectrum components [11,12,53] were previously proposed as descriptors for such linearized potentials. In this work, we favored a simpler formulation which consists in developing the descriptor space in explicit twobody, three-body and N-body interactions:

[2B] (i) n = j f n (R ij ) × f c (R ij ) (2)
[3B]

(i) (n,l) = j k f n (R ij )f c (R ij )f n (R ik )f c (R ik ) cos l (θ ijk ) (3)
[NB]

(i) (n,m) =   j f n (R ij ) × f c (R ij ) × f s (R ij )   m (4)
where R ij is the distance between atoms i and j, θ ijk is the angle centered around the atom i, and l and m are two positive integers. For the cut-off function, we chose what was originally proposed by Behler and Parrinello [8]:

f c = 1 2 (1 + cos(π(R ij /R cut ))
) with R cut being set at 6 Å. The switch function denoted f s is employed in order to prevent from non-physical behavior of the N-body contribution at short distances. To do so, two distances r 1 and r 2 are first defined as respectively 95% and 105% of the minimum of the dimer interactions and then a function is constructed to smoothly go to from 0 to 1 in the range of [r 1 : r 2 ]: f s (u) = 6u 5 -15u 4 + 11 3 where u = (R ij -r 1 )/(r 2 -r 1 ) [54]. Altogether, these expressions allow for a direct computation of the forces as well as the stress tensors by differentiating with respect to the positions. [10] Regarding the basis of functions f n , there are no physical restrictions. In particular, for the twobody interactions, one can tune f n to mimic traditional interatomic potentials as for example Morse, Lennard-Jones, Buckingham or Yukawa potentials or use simple functions like Gaussians, Lorentzian or Asymmetric lognormal functions. Likewise, for the three-body interactions, the current formulation is very similar to what is done in the Stillinger-Webber potential [55]. Finally, the N-body interaction is a generalized form of the EAM potential where the embedding function is a polynomial of the atomic density [56]. In this case, the integer m corresponds directly to the degree of N-body order. The difference between our formulation and the most recent ones proposed by Seko et al. is that, in the N-body interactions, we did not include any explicit angular dependence and did not mix different forms of the atomic density.

For the fitting procedure of such a linear model, previous studies have proposed the use of genetic algorithm [12,53], weighted ordinary least squares [11], Bayesian linear regression [39], ridge regression [18,51,52] and Least absolute shrinkage and selection operator (Lasso) [10,17]. In order to construct a simpler MLIP, we employed the Lasso regression with the Least Angle Regression Scheme (Lars, together denoted LassoLars) [57]. In practice, along with the ordinary least square objective function, χ 2 OLS , the Lasso scheme adds a penalty on the sum over the absolute value of the coefficients ω n and the employed error function is therefore given by:

χ 2 = χ 2 OLS + α n |ω (n) | ( 5
)
where α is a parameter that controls the degree of penalty. The penalty on the absolute value of the coefficients enforces lots of the linear coefficients to be exactly 0. Additional, using Lars allows us to select the most relevant descriptors by measuring their correlation to the target. Using LassoLars as a regression scheme is at the expense of accuracy and flexibility for the MLIP but it allows for a considerable reduction in the complexity of the potential [Please see the Appendix for further motivation regarding the choice of LassoLars and for additional details on the Φ-LassoLars implementation].

B. DFT database

Building a general purpose potential for Au-Fe nanocrystals would require to model the atomic interactions not only in different structures (crystal polymorph, interfaces and liquid for instance) but also in different chemical orderings including alloyed, Janus and core-shell nanoparticles. Because the focus of this work is to measure the issues related to transferability, we purposely employed an incomplete database made of only three types of nanoparticles: (1) Alloys with almost equimolar compositions, (2) Pure iron in the body-centered cubic (bcc) phase, (3) Pure gold in the face centered cubic (fcc) phase[see Fig. 1]. Then, molecular dynamics (MD) simulations by means of a house code were carried out to melt the constructed nanoparticles. Simulations were performed in the NVT ensemble obtained with Andersen thermostat at 1400 K during 500 ps using a timestep of 1 fs. We used simple pair-wise potentials made of Lennard-Jones and Morse interactions for gold and iron respectively and of Lennard-Jones interaction for the gold-iron cross-interaction. The employed Morse parameterization is the one of Hung et al. [58] while the two Lennard-Jones potentials for gold-gold and iron-gold were simply parameterized in order to match the bulk lattice parameters and the cohesive energy. Along the melting path, we extracted configurations that are representative of the solid to liquid transition. In addition, each initial structure was also manually compressed. The distances are reduced by a factor of 75% and along the compression, we extracted ten configurations in order to sample structures of higher density and better reproduce the repulsion at short distances. For the same reason, diatomic molecules FeFe, AuFe and AuAu with distances down to 1 Åwere also added in the database. For each of these configurations, forces were finally computed at the DFTlevel using single-point calculations. Spin-polarized DFT calculations were performed with the VASP code [59], using PAW type pseudopotentials for iron and gold [60], a plane wave cutoff of 650 eV and a Methfessel-Paxton smearing parameter σ of 0.01 eV. All calculations were done at the Γ-point of the Brillouin zone. Altogether, the database is made of 181653 atomic configurations with an almost equal proportion in the three types of nanoparticles (34% of alloy, 34% of pure gold and 32% of pure iron). We note that the database sampling was made with classical interaction potentials which is less satisfying than performing ab initio MD. Yet, the advantage is that it allows us to quickly sample configurations that remain physically valid and it prevents from running very computationally demanding simulations for those relatively large nanoparticles (more than 50 atoms and more than 100 for respectively the alloy and the pure nanoparticles). In addition, by employing configurations that span from the crystalline to the liquid regime, we assure a large variety of atomic neighborhoods with forces ranging from 10 -4 to 5 eV/ Å. Finally, during the manual compression of the nanoparticle and when computing the forces for diatomic molcules, it may happen that some atoms are very close thus leading to very large forces. In the χ 2 , those large forces would contribute a lot while being very unlikely to emerge in a realistic dynamics. Therefore, before performing the fitting procedure, the database is filtered out to remove cases where the forces are lower than 5 eV/ Å. Such value was chosen large enough to keep some part of the repulsive interaction but not too large to avoid a fitting that is focused only on the repulsive part.

Within the fitting procedure, the database is randomly divided in a training (95%) and a test (5%) sets.

III. RESULTS

A. Influence of the descriptor space

First, five different types of descriptors were tested. In particular, we used three functions that are "peak" functions ie. Gaussian, Lorentzian and Log-normal peaks and two functions that are usually employed for orbital calculations and that diverge at short distance ie. Slattertype (STO) and Gaussian-type (GTO) orbital [Please see Table I]. Then, the LassoLars method is employed to obtain five different interaction potentials using α = 10 -7 . For the three-body and the N-body interactions, we used respectively l = [1, 2, 3, 4, 5] and m = [START_REF] Stone | The theory of intermolecular forces[END_REF][START_REF] Behler | [END_REF]6,7]. Altogether, we have around a thousand available descriptors which include Au-Au, Fe-Fe and Fe-Au interactions for each type of descriptors. In Fig. 2.a, we show the fitting error measured as the root mean square error (RMSE) for the five different types of descriptors. RMSE measured on training and test data sets are similar which means that at this stage, no overfitting is observed. In addition, it appears that the two orbital functions are not as good as the peak functions although the STO still give an RMSE equal to 0.17 eV/ Å. Gaussian and Lorentzian descriptors are able to reproduce the forces with an RMSE respectively equal to 0.13 and 0.14 eV/ Å. Such values are similar to what is obtained with the generally employed MLIP methods including neural networks, Gaussian approximation method and linearized potentials [61]. In Fig. 2.b, the number of non-zero coefficient is plotted for the five different types of descriptors. In the cases of GTO and STO functions, much fewer descriptors were selected in comparison to the peak functions. In overall, it remains that the LassoLars algorithm allows one to drastically decrease the number of employed descriptors with respect to the number of available descriptors.

In addition, before studying the transferability issues, we wish to illustrate a second advantage of using the Las-soLars algorithm which consists in having a penalizing parameter that controls both the accuracy and the complexity of the obtained potential. According to Fig. 3, by increasing α, the number of non-zero coefficients and the computational cost can be reduced at the expense of increasing the RMSE. As such, with the LassoLars algorithm, one can simply choose which degree of accuracy or complexity is required for their usage. Finally, the presence of a plateau for the smallest values of α shows that the LassoLars regression only selects relevant descriptors thus reducing the potential complexity. Similar to what was obtained previously using α = 10 -7 , we note that:

(1) the 3 peak functions are the most accurate and behave similarly and (2) the STO function gives slightly higher RMSE yet with much fewer non-zero coefficients. Fig. 3 evidences that the LassoLars algorithm gives the ability to finely control the complexity of the potential at the expense of the accuracy on the DFT database. In the following, we will test how this complexity can influence the MLIP transferability.

B. Complexity vs Transferability

For that purpose, three additional morphologies of gold-iron nanoparticles were designed: two Janus and one core-shell [see Fig. demixed structures. Fig. 4(b) shows the corresponding RMSE on the forces without having trained the potential on these structures and using only Gaussian functions. Fig. 4(c) shows that most of the errors are located at the gold/iron interface which was not included in the training database. Surprisingly, the RMSE behavior is nonmonotonic with a minimum located for the three structures around α = 8 × 10 -6 which corresponds to 90 nonzero coefficients. More specifically, when transferring the obtained potential to Fe100Au100 Janus nanoparticles, the RMSE obtained with α = 10 -7 is 30% higher than with the less complex potential that was obtained with α = 8 × 10 -6 . Therefore, while increasing the complexity of the potential leads to a better agreement within the training database, it does not necessarily lead to an improvement of the RMSE in unlearned structures using in our case different chemical ordering (ie. Janus and Core- shell instead of alloy nanoparticles). Our challenging test demonstrates that precautions should be made when using machine-learning approaches and that increasing the complexity does not automatically lead to a better overall potential.

According to Fig. 5, the non-monotonic behavior that was highlighted when using only Gaussian functions is also observed with the two other peak functions ie. Lorentzian and Log-normal. We also note that, GTO is again very inaccurate and should most probably be avoided for usage as descriptor. Finally, the STO functions while being less accurate on the training database seems to give in overall better results in terms of transferability.

Finally, a combination of three different descriptors (Gaussian, Gaussian-type orbital and Lorentzian peak) is tested in order to employ simultaneously two different types of peak functions along with a function that diverge at short repulsion. STO was chosen for its remarkable ability to decrease the number of non-zero coefficient while Lorentzian peaks also showed slightly less non-monotonocity. On the one hand, regarding its fitting performance, this combination gives an RMSE similar to that obtained previously yet with fewer non-zero coefficients [See Fig. 5]. Having these two additional functions gives more flexibility in the descriptor space and fewer functions can therefore be selected for the same accuracy. On the other hand, for the transferability towards the unlearned structures, the combination gives RMSE that are comparable to the peak functions and even better in two cases (ie. Fe110Au111 and Core-Shell). More importantly, it allows for a further reduction of the nonmonotonicity. We note that being able to combine three different types of descriptors at the same time is an additional advantage of using a constrained linear regression scheme such as LassoLars.

C. Quality of the potential

Even if the aim of the paper was not to obtain an allpurpose MLIP for gold-iron interactions, we still wish to asses the quality of the obtained potential. Regarding the force accuracy on the learned database, Fig. 6(a.b) shows the correlation plot obtained with two different values of α. For the RMSE, we obtained in both cases values of 0.14eV/ Å when using the three descriptors at the same time. Such value is on par with most of the currently employed MLIP methods [27] and by comparison, the EAM potential that was recently developed for Au-Fe nanoparticle [62] gives a value of 1.4 eV/ Å. Our MLIP is thus already a drastic improvement in force evaluation for the studied nanoparticles [See Fig. 6.c].

Furthermore, additional simulations with the MLIP were performed to check some properties in the bulk phase although it was not included in the training set. Simulations were carried out using the large-scale molecular dynamics software LAMMPS [63] in which the proposed MLIP was implemented. Periodic boundary conditions were employed and the different minimization runs were performed down to a net force of 10 -6 eV/ Å. We measured eight different lattice constants (pure iron bcc, pure gold fcc, alloys with 25%, 50% and 75% of gold in both bcc and fcc phases). In addition, our fitting did not include any energy. Therefore, to take into account the atomic energy, the MLIP energies are shifted in order to match the cohesive energy of pure iron and gold most stable states ie. bcc and fcc respectively. Then, we also measured the cohesive energy for each alloying structures. Results are compared in Fig. 6.(d,e) to DFT calculations that were previously obtained. [47] The errors are lower than 3% for lattice spacing and 12% for cohesive energy which is already satisfying considering that in the database, we employed alloying proportions that are much closer to 50% and only used nano-crystals. Therefore, being able to reach such small relative errors even for those extreme alloying proportions (25% and 75%) and in addition with bulk structures is an encouraging result for our MLIP.

Finally, we measured phonon dispersion curves using supercell approach implemented in PHONOPY [64] [See Fig. 6(e,f)]. In practice, the dynamical matrix was obtained by moving each symmetrically independent atoms by 0.01 Å. We used a supercell of size 4×4×4. The agreement between DFT and MLIP curves is not as good as what is usually obtained in MLIP works [17,65] but it remains qualitatively satisfying if one considers that our MLIP was not trained on any bulk structures.

Even if our potential has not been designed to reproduce bulk properties, these results are already very encouraging. Obtaining an accurate MLIP potential that is transferable to any phase and/or structure is a considerable challenge for multi-component systems and would require to carry out additional DFT calculations to build a bigger database including bulk structures but also interfaces. Besides, in order to target a specific application, one should also perform "on-the-fly" optimization of the potential as proposed by Jinnouchi et al. [38,39]. However this is not the main purpose of this paper, which focuses on shedding light onto the relationship between complexity and transferability in machine learning force fields. It remains that our current MLIP potential may be used as a first step when studying bimetallic Fe-Au nanoparticles.

IV. DISCUSSION AND CONCLUSION

To summarize, this work aims at measuring transferability issues that can occur when using MLIP in untrained structures. To begin, we presented the employed MLIP method that includes a linearized potential and a penalizing regression scheme. Then, we discussed the influence of the descriptor space and showed that although the three different peak functions behave similarly in terms of accuracy and number of non-zero coefficients, the repulsive functions lead to a worst accuracy but with a lower number of non-zero coefficients. Next, we demonstrated that by using the LassoLars algorithm instead of the previously employed linear regression scheme, one can finely tune the complexity of the potential along with its accuracy. This is because the penalizing parameter α allows for turning off most of the initially proposed descriptors thus controlling the overall complexity of the potential. With this ability in hand, we measured transferability issues using three unlearned structures that are qualitatively different from what was considered in the training database. We showed that while the accuracy on the trained structures decrease monotonically as the value of α is decreased, it is not the case for those untrained structures. Indeed, when using only one type of descriptor, it exists an optimal value of α that allows for the best transferability. Finally, we introduced a way to overcome this transferability issue which consists in using different types of descriptors simultaneously. This again shows an other advantage of using the LassoLars algorithm which is able to actively select the most appropriate descriptors. Finally, we computed some properties of the Fe-Au in bulk and showed that the obtained potential is already qualitatively satisfying. But, before being able to really our potential from practical applications, we plan to improve further it by adding bulk and interface DFT calculations within the database and by implementing "on the fly" learning. As a perspective, the obtained MLIP will be further improved and then employed to study nucleation in core-shell FeAu nanoparticles as observed in experiments. [41][42][43][44][45] Before closing, we wish to discuss two additional points.

First, the very intuitive MLIP expression that was employed here allows us to give some insights on the nature of the interactions. Here, we focus on the combined MLIP where Gaussian, STO and Lorentzian peak were used simultaneously. We can indeed distinguish between each descriptors (Gaussian, STO, Lorentzian). For that purpose, we compute the ratio between the absolute value of the forces given by each descriptor and the sum of the absolute values of the three descriptors and then perform an average over each force components (x,y,z) of all of the atoms within the training database (alloy, pure iron and pure gold nanoparticles). Fig. 7.a shows that the preponderant functions are different depending on the considered interactions thus highlting the advantage of using simultaneously the three types of descriptors. More interestingly, the same can be done in order to distinguish between two-body, three-body and N-body contributions [see Fig. 7.b]. In our case, the two-body and three-body contributions are more important than the N-body contributions. This may explain why the previously employed EAM potentials that do not possess any explicit angular contributions could not accurately compute the forces.

Moreover, we would like to raise an additional implication of our work in which the transferability issues were measured and connected to the complexity of the potential. Indeed, as previously discussed, when using the LassoLars regression scheme, the complexity of the potential can be adjusted using the penalizing parameter. In the alternative MLIP methods, the same is done by modifying (1) the number of neurons and hidden layers in the case of neural-network potential [START_REF] Behler | [END_REF] and (2) the number of selected configurations after sparsification in the case of gaussian approximation model [9,32]. For some users of these techniques, the rule of thumb may be to use these adjusting parameters in order to increase the complexity of the potential which necessarily improves the accuracy on the learning database. Yet, our work indicates that transferability issues should be expected by such operation.

APPENDIX A. Advantage of using LassoLars against other linear regression scheme

In this section, we wish to further motivate the choice of our regression scheme by working on a case study where we will compare LassoLars with two commonly employed linear regression schemes naming Ridge and the coordinate descent Lasso. For this test, the database is generated with an equimolar mixture of binary Lennard-Jones particles thus allowing to directly verify the obtained MLIP. With such a binary system, an additional challenge for the fitting algorithm is to distin-guish self-species and cross-species interactions. In practice, positions and forces were measured for 50 configurations of 64 atoms in the liquid regime. Regarding the basis of descriptors, only two-body interactions were considered and we used 17 Lennard-Jones functions with different distance parameters including those in the original simulations. All of the four employed methods manage to retrieve a linear combination of Lennard-Jones functions that matches the original interactions. However, it appears that only LassoLars can find the correct coefficients ω n setting all coefficients to 0 except those of the original interactions [see Fig. 8]. Such result shows the advantage of using LassoLars instead of the commonly employed linear regression methods.

B. Numerical implementation of the Φ-LassoLars method

In this section, we give some additional details on the described MLIP. First, obtaining an MLIP using the Φ-LassoLars method consists in two steps. In the first step, a homemade C++ parallelized code using OpenMP was developed to construct a matrix:

• Each columns of the matrix designates a specific descriptor which can be 2B, 3B or N B for all the functions and considered values for their parameters described in Table 1.

• Each rows of the matrix correspond to the force on a given direction, atom and structure.

In the second step, a python code concatenates the obtained matrices for each structures and read the associ-ated DFT forces. The same python code finally employs the LassoLars method as implemented in the sklearn package to obtain the linear coefficients associated to each columns of the matrix. Then, in order to use the obtained MLIP, the same C++ code is employed to read the obtained coefficients and generate input files for LAMMPS simulations. Those consists on three parts:

1. For the 2B interactions, we directly add all of the selected linear contributions and generate a table file that can be read by LAMMPS using pair style table.

2. For the 3B interactions, we build a homemade routine that is added to LAMMPS and use an personalized input file.

3. For the NB interactions, we use a python code based on atsim.potentials [START_REF] Rushton | atsim.potentials -potential model tabulation for atomic scale simulation[END_REF] to generate EAMlike files that can be directly read by LAMMPS.

Finally, the pair style hybrid/overlay is employed to combine all the contributions.

Regarding the computational timing for building the a matrix for a structure of 64 atoms containing all of the functions (2B, 3B, NB for all values in Table 1) for one type of descriptors is obtained in approximately 3 minutes. This process can be parallelized since each structure can be treated independently. Then, after reading all of the input matrices, a LassoLars fitting takes less than 20 seconds when using only one type of functions and less than 2 minutes when using simultaneously three types of functions. Results are obtained on one Intel E5-2650 processor.

ACKNOWLEDGEMENT JL acknowledges financial support of the Fonds de la Recherche Scientifique -FNRS. Computational resources have been provided by the Consortium des Equipements de Calcul Intensif (CECI), by the Fédération Lyonnaise de Modélisation et Sciences Numériques (FLMSN) and by the Regional Computer Center CALMIP in Toulouse (grant n • p1141). JL thanks David Tew for having introduced him to the machine-learning potential problem and C. Patrick Royall for supervising the early method development. JL is also grateful to Francesco Turci, Joshua F. Robinson and Atsuto Seko for fruitful discussions. Authors are finally grateful to Abdul-Rahman Allouche and Albert P. Bartók for proofreading the manuscript.

FIG. 1 .

 1 FIG. 1. Images of the initial structures employed in the database.

FIG. 2 .

 2 FIG. 2. (a) RMSE obtained for five different descriptors and measured on the training and the test sets. (b) Number of non-zero coefficients obtained with α = 10 -7

  FIG.3. Influence of the penalizing factor α on testing RMSE and on (b) the number of non-zero coefficients.

FIG. 4 .

 4 FIG. 4. (a) Images of untrained configurations with morphology Janus and Core-shell, (b) Fitting error as a function of fitting parameter α for each untrained configurations and for the training set and (c) Map of the force errors on the Fe100Au100 Janus structure for different values of α. The green line designates the iron/gold interface.

FIG. 5 .

 5 FIG. 5. (a) Number of non-zero coefficients and (b) RMSE as a function of the constrain parameter α when using simultaneously Gaussian, STO and GTO functions as descriptors. The dotted lines correspond to results when using only Gaussian functions ie. results from Fig. 4b.

FIG. 6 .

 6 FIG. 6. (a-c) Correlation plot between DFT forces and forces from MLIP with α = 10 -7 (b) and with α = 5 × 10 -6 (c) and from EAM (d). Blue and red points correspond to results within the training and the test sets. (d) Relative errors obtained on the lattice parameters. (e-f) Phonon dispersion of respectively Gold FCC and Iron BCC. The plain lines correspond to MLIP results while the dotted lines are obtained by DFT calculations.

FIG. 7 .

 7 FIG.7. Force contributions for each types of (a) descriptors and (b) multi-body components averaged over alloy, pure iron and pure gold nanoparticles.

FIG. 8 .

 8 FIG.8. Values of the obtained coefficients ωn using different linear regression scheme: (a) Ridge, (b) Lasso, (c) LassoLars. The penalty parameter was set to 10 -5 . The red points and the blue line correspond respectively to the original interactions and the fitting results.

TABLE I .

 I Summary of the tested descriptors.

	Function name	Equation	List of an	List of bn	N func
	Gaussian	fn(Rij) = exp(-an(Rij -bn) 2 )	0.5, 1.0, 1.5	1, 2, 3, 4, 5, 6	990
	Lorentzian peak	fn(Rij) = 1/((Rij -an) bn + 1).	1, 2, 3, 4, 5, 6	2, 4, 6, 8	1320
	Log-normal peak	fn(Rij) = exp(-ln[(Rij -an)/bn] 2 )	1, 2, 3, 4, 5, 6	1.0, 1.5, 2.0, 2.5	1320
	Slatter-type orbital	fn(Rij) = R an ij exp(-bnRij)	-2, -1, 0, 1, 2	2, 4, 6, 8, 10	1375
	Gaussian-type orbital	fn(Rij) = R an ij exp(-bnR 2 ij )	-2, -1, 0, 1, 2	2, 4, 6, 8, 10	1375


Andrea Díaz-Barriga Yáñez

Auriane Couderc

Léa Longo

Annabelle Merchie

Hanna Chesnokova

Emma Langlois

Catherine Thevenot
email: catherine.thevenot@unil.ch

Jérôme Prado
email: jerome.prado@univ-lyon1.fr

Learning to run the number line: The development of attentional shifts during single-digit arithmetic

Keywords: arithmetic, spatial attention, numerical cognition, development

Recent studies suggest that solving single-digit subtraction and addition problems is associated with left and right shifts of attention in adults. Here, we explored the development of these spatial shifts in children from 7 to 11. In two experiments, children were asked to solve single-digit addition, subtraction, and multiplication problems in which operands and arithmetic sign were shown sequentially. Although the first operand and the arithmetic sign were presented on the center of a screen, the second operand was presented either to the left or to the right visual field. We found that subtraction problems were increasingly associated with a leftward bias from 7 to 11, such that problem solving was facilitated when the second operand was in the left visual field around age 10-11 (i.e., grade 5). We also found evidence that at least some children may already associate addition problems with the right side of space by age 9-10 (i.e., grade 4). No developmental increase in either leftward or rightward bias was observed for multiplication problems. We argue that these attentional shifts reflect the progressive automatization of calculation procedures that involve specific movements to the left or right of a sequential representation of numbers during subtraction and addition.

Introduction

Increasing evidence indicates that our ability to process numbers is grounded in spatial representations [START_REF] Fischer | Spatial biases in mental arithmetic[END_REF] . For example, small numbers are associated with the left side of space and large numbers with the right side of space. [START_REF] Dehaene | The number sense. How the mind creates mathematics[END_REF][START_REF] Dehaene | The Mental Representation of Parity and Number Magnitude[END_REF][START_REF] Wood | On the cognitive link between space and number: a meta-analysis of the SNARC effect[END_REF] Numbers also automatically bias spatial attention, such that targets are detected faster in the left visual field (LVF) when they follow small numbers, and faster in the right visual field (RVF) when they follow large numbers. [START_REF] Dodd | Attentional SNARC: There's something special about numbers (let us count the ways)[END_REF][START_REF] Fischer | Perceiving numbers causes spatial shifts of attention[END_REF] Taken together, this suggests that numerical magnitudes may be represented on a mental number line (MNL) that is organized horizontally and in ascending order from left to right. [START_REF] Dehaene | The Mental Representation of Parity and Number Magnitude[END_REF][START_REF] Gevers | The mental representation of ordinal sequences is spatially organized Wim[END_REF] Recent studies suggest that such spatial associations may not be limited to number processing per se, but might also be observed during arithmetic calculation in adults. For example, when adults are asked to estimate the result of a subtraction or an addition problem, they tend to underestimate the result of a subtraction and overestimate the result of an addition. [START_REF] Knops | Dynamic representations underlying symbolic and nonsymbolic calculation: Evidence from the operational momentum effect[END_REF][START_REF] Mccrink | Moving along the number line: Operational momentum in nonsymbolic arithmetic[END_REF] Several explanations for this effect (termed operational momentum, OME) have been proposed. [START_REF] Chen | Spatial intuition in elementary arithmetic: A neurocomputational account[END_REF][START_REF] Mccrink | Operational momentum in large-number addition and subtraction by 9-month-olds[END_REF][START_REF] Pinhas | Mental movements without magnitude? A study of spatial biases in symbolic arithmetic[END_REF] Nevertheless, one major explanation is that participants rely on attentional shifts along the MNL to estimate results of problems. In a recent study, we argued that such attentional shifts might occur during exact symbolic arithmetic as well. [START_REF] Mathieu | Running the number line: Rapid shifts of attention in single-digit arithmetic[END_REF] Specifically, we asked adults to solve single-digit arithmetic problems presented on a computer screen. While the first operand and the arithmetic sign were presented sequentially at the center of the screen, the second operand was presented either in the LVF or the RVF (see Fig. 1). Results indicated that participants were faster to solve addition problems when the second operand was presented in the RVF than the LVF, while they were faster to solve subtraction problems when the second operand was presented in the LVF than the RVF. Thus, spatial shifts of attention are elicited during exact single-digit arithmetic in adults.

Importantly, no spatial bias was observed in multiplication problems (which are explicitly learned by rote in school). Thus, it is possible that these horizontal shifts of attention reflect calculation procedures relying on left-right movements along the MNL that have been automatized after years of practice with arithmetic calculation. [START_REF] Fayol | The use of procedural knowledge in simple addition and subtraction problems[END_REF] To date, few studies have investigated the associations between arithmetic and space in children. However, available evidence suggests that these associations may be relatively late developing. For example, studies have failed to find an OME in 6-and 7-year-olds [START_REF] Knops | Examining the presence and determinants of operational momentum in childhood[END_REF] as well as in 9-year-olds. [START_REF] Kucian | Operational momentum effect in children with and without developmental dyscalculia[END_REF] It has thus been suggested that associations between space and arithmetic might be the result of years of education and that an "unconscious shift of attention on the MNL becomes evident only with increasing expertise and automatization". [START_REF] Kucian | Operational momentum effect in children with and without developmental dyscalculia[END_REF] This proposal is consistent with the idea that automatic leftward and rightward movements along a sequential representation of numbers might stem from the repeated use of deliberate counting during arithmetic learning in elementary school. [START_REF] Fayol | The use of procedural knowledge in simple addition and subtraction problems[END_REF][START_REF] Barrouillet | On the problem-size effect in small additions: Can we really discard any counting-based account?[END_REF][START_REF] Thevenot | Ten-year-old children strategies in mental addition: A counting model account[END_REF][START_REF] Uittenhove | Fast automated counting procedures in addition problem solving: When are they used and why are they mistaken for retrieval?[END_REF] Yet, to our knowledge, there is no evidence that increasing practice with arithmetic is associated with the emergence of shifts of attention along the MNL during mental calculation in children.

To test this hypothesis, we explored in two experiments the developmental emergence of spatial biases during mental calculation in children from 7 to 11. Using the same paradigm as Mathieu and colleagues [START_REF] Mathieu | Running the number line: Rapid shifts of attention in single-digit arithmetic[END_REF] , one group of children was presented with single-digit addition and subtraction problems (Experiment 1) and another group with single-digit addition and multiplication problems (Experiment 2). We expected addition problems to be increasingly associated with a rightward bias with age. Conversely, we expected subtraction problems to be increasingly associated with a leftward bias with age. However, there should be no emergence of associations between multiplication problems and space, as these problems are directly retrieved from memory and may not involve calculation procedures relying on shifts along the MNL.

Experiment 1

Materials and Methods

Participants

Procedure

All children were tested in one single session that lasted approximately 30 minutes.

Children from Wave 1 were tested in a quiet area of their school, whereas children from Wave 2 were tested either in the lab or in a quiet area of their school. Each child was tested individually and gave verbal assent to participant. Parent's written consents were also obtained. The testing session started with the Math Fluency test from the Woodcock-Johnson battery [START_REF] Woodcock | Woodcock-Johnson III Tests of cognitive abilities[END_REF] to assess children's arithmetic skills (see below), followed by the experimental task. 13

Measures

Math fluency: Children's arithmetic skills were first assessed with the Woodcock-Johnson III Math Fluency subtest. The Math Fluency is a paper-and-pencil test that includes single-digit addition, subtraction, and multiplication problems. Children needed to solve as many problems as possible within 3 minutes. [START_REF] Woodcock | Woodcock-Johnson III Tests of cognitive abilities[END_REF] Experimental task: The experimental task was computer-based and adapted from Mathieu and colleagues. [START_REF] Mathieu | Running the number line: Rapid shifts of attention in single-digit arithmetic[END_REF] Small arithmetic problems included pairs of non-identical operands between 1 and 5 [(2, 1); (3, 1); (3, 2); (4,1); (4, 2); (4, 3); (5, 1); (5, 2); (5, 3); (5, 4)], and large arithmetic problems included pairs of non-identical operands between 5 and 9 [(6, 5); (7, 5); (7, 6); (8, 5); (8, 6); (8,7); (9, 5); (9, 6); (9, 7); (9, 8)]. Both small and large problems contained the number 5 in order to have the same number of problems in both categories. These pairs of non-identical operands were used to construct 20 addition problems (10 small problems and 10 large problems) and 20 subtraction problems (10 small problems and 10 large problems) with their second operand presented once to the right and once to the left, resulting in a total of 80 problems (i.e., 40 addition problems and 40 subtraction problems). For both addition and subtraction problems, the larger of the two operands was always presented first. This ensured that results from the subtraction problems were positive and that children could not anticipate the type of problem to come. The presentation of the trials was pseudo-randomized, so that no more than three problems of the same type (i.e., same size and spatial location of the second operand) would appear consecutively. Four scenarios were created by generating 4 random lists of trials. Each scenario was separated into 2 runs of 40 operations. The experiment systematically started with a practice run of 8 trials including tie problems (e.g. 7 -7), problems with 0 (e.g. 3 + 0), and problems with small and large operands (e.g. 7 + 3).

In each trial, operands and arithmetic sign were displayed in white Times New Roman 36-point font on a black background (see Fig. 1). Each trial started with the presentation of a white fixation dot for 500 ms followed by the presentation of the first operand that lasted an additional 500 ms on screen. After a first delay of 500 ms, either a + or a -sign appeared on screen for 150 ms at the center of the screen. A second delay of 300 ms separated the disappearance of the arithmetic sign from the second operand. This delay was chosen because it was the delay for which arithmetic shifts were maximal in Mathieu et al., 2016. The second operand was then displayed for 150 ms, either 5º to the left or 5º to the right of the center of the screen. Children needed to solve the calculation verbally in less than 5 seconds, otherwise the software automatically moved on to the next arithmetic problem.

All response times (RTs) were recorded through a headset microphone and corresponded to the period between the presentation of the second operand and the onset of the answer. The experiment was controlled by the DmDX software [START_REF] Forster | DMDX: A windows display program with millisecond accuracy[END_REF] and RTs were checked off-line and manually adjusted with CheckVocal 22 for each participant. Before the experiment started, children were given the instructions printed in a A4 white sheet and were able to ask any questions if needed. Children were located 44 cm from the 15-inch computer screen. A chin-rest was used to avoid head movements.

Results and discussion

Only correct responses were analyzed. Correct trials constituted 90.33% of trials for small problems, but only 62.78% of trials for large problems. To maximize power and ensure that our results would not be confounded by differences in accuracy between older and younger children, we exclusively focused our analyses on small problems. Small addition problems (91.5%) were responded more accurately than small subtraction problems (89.15%) (t (103) =2.09, p =. 04). For each participant and operation, we subtracted the mean RT of trials in which the second operand appeared on the right from the mean RT of trials in which the second operand appeared on the left. This difference in RT (dRT) served as dependent variable in the following analyses (see Ref. 13). Outliers (dRT smaller or greater than 2.5 SDs from the mean for each participant in either small addition or subtraction problems) were removed from the analyses (this corresponded to 3 children).

Descriptive statistics

Descriptive statistics for age, Math fluency scores, dRTs for small addition and subtraction problems, as well as RTs for small addition and subtraction problems are shown in Table 1.

Preliminary analysis

Independent t tests were conducted to test whether there was an effect of Wave (Wave 1 vs Wave 2) and gender on dRTs for small addition and subtraction problems. Results revealed no Wave (t(99)=1.09, p=.28 for dRT in small addition problems, and t(99)=.08, p=.94 for small subtraction problems) or gender (t(99)=.46, p=.64 for dRT for small addition problems, and t(99)=.67, p=.50 for small subtraction problems 1 ) effects in any of the dRTs.

dRTs analyses

As children get older (and calculation procedures would become more automatized), we predicted that dRTs would become more positive for addition problems (i.e., an increase in rightward bias) and more negative for subtraction problems (i.e., an increase in leftward bias). Thus, dRTs for addition and subtraction problems were analyzed in multiple regression analyses with age as the main predictor. Because age was positively correlated with math fluency scores (r=.30, p<.01, two-tailed) and negatively correlated with overall RTs for small addition problems (r=-.22, p<.05, two-tailed) and small subtraction problems (r=-.20, p<.05, two-tailed), these predictors were also included in the multiple regression analyses in order to investigate the amount of unique contributions of age to differences in dRTs (see Table 2). [START_REF] Fischer | Spatial biases in mental arithmetic[END_REF] Levene's test was significant, suggesting a violation of the equal variance assumption, thus a Mann-Whitney test was conducted. Results also revealed no gender effect, U=1238, p=.32, z=-.25 for small addition problems, and U=1190, p=.98, z=-.58 for small subtraction problems.

For both multiple regression analyses, the Durbin-Watson test was found to be within acceptable parameters (2.34 and 2.13 respectively), suggesting that the assumption of independent errors was met and that results may apply to a wider population. [START_REF] Field | Discovering Statistics Using SPSS[END_REF] Multicollinearity was also systematically checked using the variance inflation factor (VIF).

The largest VIF was well below 10 for both regressions. Similarly, the tolerance data were all within acceptable boundaries (all greater than 0.1). Therefore, it was concluded that there was no collinearity within the data. [START_REF] Field | Discovering Statistics Using SPSS[END_REF] Over and above differences in math fluency and overall RTs, dRTs for small addition problems did not increase with age (t (97) =-.17, p=.57, one-tailed 2 ) (see Table 2 and Fig. 2A). Such a lack of positive relationship between age and dRTs for small addition problems is not consistent with our hypotheses, as we predicted an increase in rightward bias for addition problems while children get older and therefore should automatize addition procedures. One possibility is that the bias might already be present in the youngest children investigated (such that the rightward bias might develop earlier than 3 rd grade). To investigate this possibility, we tested whether dRTs for small addition problems were larger than 0 in children from 3 rd grade (i.e., approximately age 8-9), 4 th grade (i.e., approximately age 9-10), and 5 th grade (i.e., approximately age 10-11) in our sample. One-sample t-tests indicated that dRTs for small addition problems were not larger than 0 in children from 3 rd (t (46) =.92, p=.18, one-tailed, d=.13), 4 th (t (32) =-1.46, p=.92, one-tailed, d=-.25), or 5 th grade (t (20) =.14, p=.45, one-tailed, d=.03). Therefore, contrary to our hypothesis, we did not find any evidence for either an increase or an early presence of a rightward bias in small addition problems in this sample of elementary school children. Although it is always difficult to interpret a null result, potential reasons for this lack of effect will be discussed in the general discussion.

In line with our hypotheses, however, we found a significant negative relationship between age and dRTs for small subtraction problems (t (97) =-1.89, p=0.03, one-tailed [START_REF] Dehaene | The Mental Representation of Parity and Number Magnitude[END_REF] , see Table 2). In other words, dRTs for small subtraction problems were more negative in older than younger children (see Fig. 2B). To further explore this effect, we tested whether dRTs [START_REF] Dehaene | The number sense. How the mind creates mathematics[END_REF] Our hypothesis was unidirectional because we only expected dRTs for small addition problems to increase (and not decrease) with age. [START_REF] Dehaene | The Mental Representation of Parity and Number Magnitude[END_REF] Our hypothesis was unidirectional because we only expected dRTs for small subtraction problems to decrease (and not increase) with age.

for small subtraction problems were smaller than 0 in children from grades 3, 4, and 5 in our sample. We found that dRTs for small subtraction problems were not smaller than 0 in 3 rd (t (46) =1.93, p=.97, one-tailed, d=.28) and 4 th grade (t (32) =.70, p=.76, one-tailed, d=.12).

However, dRTs for small subtraction problems were significantly smaller than 0 in 5 th grade (t (20) =-2.36, p=.01, one-tailed, d=-.51). Therefore, our results provide evidence for an increase in leftward bias when children solve subtraction problems over the course of elementary school, such that children exhibit shorter response times when the second operand is on the LVF (as compared to the RVF) around the ages of 10-11 (i.e., 5 th grade). This finding is in line with the proposal that, as children get older and practice arithmetic problem solving, they may automatize calculation procedures relying on attentional shifts along the MNL. [START_REF] Mathieu | Running the number line: Rapid shifts of attention in single-digit arithmetic[END_REF] In the case of subtraction, this would involve shifting attention to the left of the MNL, such that the " -" sign may automatically direct attention to the left side of space during problem solving.

These shifts of attention, however, should not be observed if answers are not calculated but directly retrieved from memory. [START_REF] Fayol | The use of procedural knowledge in simple addition and subtraction problems[END_REF] For example, Mathieu and colleagues [START_REF] Mathieu | Running the number line: Rapid shifts of attention in single-digit arithmetic[END_REF] found that the position of the second operand (LVF versus RVF) did not affect solving times during multiplication problem solving. This is consistent with the idea that solving multiplication problems does not involve calculation procedures but instead direct retrieval from memory (because single-digit multiplication problems are learned by rote in school; [START_REF] Mathieu | Running the number line: Rapid shifts of attention in single-digit arithmetic[END_REF] ). Thus, in Experiment 2, we presented another group of children with multiplication (as well as addition) problems using the same paradigm as in Experiment 1.

Experiment 2

Materials and Methods

Participants

As for Experiment 1, data collection was conducted in two different time waves.

Procedure

Experiment 2 followed the same procedure as Experiment 1, with all children being tested in one single session. However, the testing session lasted between 30 and 40 minutes because the experimental task included two more runs (see below). The testing session started with the Math Fluency test from the Woodcock-Johnson battery [START_REF] Woodcock | Woodcock-Johnson III Tests of cognitive abilities[END_REF] followed by the experimental task.

Measures

The Woodcock-Johnson Math Fluency subtest used in Experiment 2 was the same as in Experiment 1. Regarding the experimental task, the pairs of non-identical operands, the classification of problems (i.e. small and large), the apparatus, and the stimulus timing were the same as in Experiment 1. Hence, only full details for the stimuli and for the experimental procedure of Experiment 2 are provided next.

For Experiment 2, the pairs of non-identical operands were used to construct 20 addition problems (10 small problems and 10 large problems) and 20 multiplication problems (10 small problems and 10 large problems) with their second operand presented once to the right and once to the left, resulting in 40 addition problems and 40 multiplication problems.

The largest operand could be presented as either the first or the second operand, such that each problem was presented once for each order of presentation of operand (first vs second position). Therefore, there were a total of 160 problems, which were separated into 4 runs of 40 operations. The trials were pseudo-randomized and 4 different scenarios based on 4 different lists were created (with 4 runs in each scenario).

Results and discussion

As in Experiment 1, only correct responses in small arithmetic problems were analyzed. Correct trials constituted 89.4% of trials for small problems, but only 56.89% of trials for large problems. There was no accuracy difference between small addition problems (89.75%) and small multiplication problems (89.05%) (t (68) =.68, p = .50). dRTs were also calculated and served as dependent variable in the following analyses. Outliers (dRT smaller or greater than 2.5 SDs from the mean for each participant in either small addition or multiplication problems) were removed from the analyses (this corresponded to 1 child).

Additionally, data from one child on the WJ Math fluency task was considered as missing in further statistical analyses due to failure to understand task instructions. Analyses were therefore conducted using all available data using pairwise deletion to avoid bias and reduction in power.

Descriptive statistics

Descriptive statistics for age, Math fluency scores, dRTs for small addition and multiplication problems, as well as overall RTs for all small addition and multiplication problems are shown in Table 3.

Preliminary analysis

Independent t tests were conducted to test whether there was an effect of Wave (Wave 1 vs Wave 2) and gender on dRTs for small addition and multiplication problems. Results revealed no Wave (t(66)=1.7, p=.10 for dRTs in small addition problems, and t(66)=-.87, p=.39 for small multiplication problems) or gender (t(66)=-1.39, p=.17 for dRTs in small addition problems, and t(66)=-1.35, p=.18 for small multiplication problems) effects in any of the dRTs. dRTs analyses dRTs for addition and multiplication problems were analyzed in multiple regression analyses with age as the main predictor of interest. Just as in Experiment 1, we expected addition problems to be increasingly associated with a rightward bias in older than younger children. We also predicted a lack of developmental increase in either leftward or rightward bias for multiplication problems. As stated before, these problems are thought to involve direct retrieval from memory and may not involve calculation procedures relying on movements along the MNL. As in Experiment 1, math fluency scores and overall RTs for either small addition or small multiplication problems were included in the multiple regression analyses. Also like in Experiment 1, for both multiple regression analyses, the Durbin-Watson test was systematically checked and was found to be within acceptable parameters (1.85 and 2 respectively), suggesting that the assumption of independent errors was met and that results may apply to a wider population. [START_REF] Field | Discovering Statistics Using SPSS[END_REF] Multicollinearity was also checked using the VIF. The largest VIF was well below 10 for both regressions. Similarly, the tolerance data were all within acceptable boundaries (all greater than 0.1). Therefore, it was concluded that there was no collinearity within the data. [START_REF] Field | Discovering Statistics Using SPSS[END_REF] Results showed that, over and above differences in math fluency and overall RTs, dRTs for small addition problems increased with age (t (63) =1.67, p=0.0499, one-tailed [START_REF] Wood | On the cognitive link between space and number: a meta-analysis of the SNARC effect[END_REF] , see Table 4 and Fig. 3A). In other words, dRTs for small addition problems were more positive for older than younger children. Further analyses showed that dRTs for small addition problems were larger than 0 in children from both 4 th (t (41) =2.96, p=.003, one-tailed, d=.46) and 5 th grade (t (25) =3.35, p=.001, one-tailed, d=.66). Thus, in this sample, children already exhibited shorter response times when the second operand was on the RVF (as compared to the LVF) in 4 th grade. Furthermore, this rightward bias increased with age. This finding suggests that, much like subtraction problem solving (Experiment 1), addition problem solving may be characterized by an increase of attentional shifts along the MNL as children get older. In the case of addition, however, these shifts might involve moving to the right of the MNL, such that the " + " sign may automatically direct attention to the right side of space during problem solving. Note, however, that these results are inconsistent with our own findings in Experiment 1, since neither a rightward bias nor a developmental increase in rightward bias was found in addition problems in Experiment 1. Potential reasons for the inconsistencies between experiments are discussed in the general discussion.

In contrast to small addition problems, we did not find any relationship between age and dRTs for multiplication problems (t (63) =.29, p=.77, two-tailed) (see Table 4 and Fig. 3B). Follow up analyses also showed that dRTs in small multiplication problems were not different than 0 in either 4 th (t (41) =-.84, p=.41, two-tailed) or 5 th grade (t (25) =.16, p=.89, two-tailed). Thus, we found no evidence of a leftward or rightward bias during multiplication problem solving in elementary school children. This is consistent with the fact that answers of multiplication problems are mainly learned by rote in school and therefore never really [START_REF] Wood | On the cognitive link between space and number: a meta-analysis of the SNARC effect[END_REF] Our hypothesis was unidirectional because we only expected dRTs for small addition problems to increase (and not decrease) with age. calculated by children. The lack of spatial bias with multiplication problems also makes it unlikely that the rightward bias observed with addition problems (in Experiment 2) is due to the size of the answer. Indeed, addition outcomes are always larger than both of the operands involved. Thus, it is possible that the rightward bias observed may simply result from an association between relatively large numbers and the right space of the MNL (e.g. [START_REF] Dehaene | The number sense. How the mind creates mathematics[END_REF][START_REF] Dehaene | The Mental Representation of Parity and Number Magnitude[END_REF][START_REF] Wood | On the cognitive link between space and number: a meta-analysis of the SNARC effect[END_REF][START_REF] Hartmann | Spatial biases during mental arithmetic: Evidence from eye movements on a blank screen[END_REF] ), rather than a calculation procedure per se. However, outcomes of multiplication problems in Experiment 2 were larger, or as large as, the operands involved (because there were no multiplications involving 0). Thus, the lack of rightward bias with multiplication problems discards the hypothesis that the rightward bias observed with addition problems is due to the size of the answer.

General Discussion

Several studies suggest that arithmetic calculation is associated with attentional shifts in adults. [START_REF] Mathieu | Running the number line: Rapid shifts of attention in single-digit arithmetic[END_REF][START_REF] Di | The Time Course of Spatial Attention Shifts in Elementary Arithmetic[END_REF][START_REF] Pinheiro-Chagas | Finger Tracking Reveals the Covert Stages of Mental Arithmetic[END_REF][START_REF] Zhu | Spatial Bias Induced by Simple Addition and Subtraction: From Eye Movement Evidence[END_REF][START_REF] Zhu | Spatial Attention Shifts in Addition and Subtraction Arithmetic: Evidence of Eye Movement[END_REF] To our knowledge, however, no study has investigated how and when these shifts emerge in children. The aim of the two experiments presented in this paper was to investigate the developmental emergence of spatial biases during mental calculation in children from 7 to 11. Children were asked to verbally solve single-digit addition, subtraction, and multiplication problems. Operands and arithmetic signs were presented sequentially. Although the first operand and the arithmetic sign were presented on the center of the screen, the second operand was displayed either in the left or in the right side of the fixation.

In Experiment 1, we found that small subtraction problems were increasingly associated with a leftward bias from 7 to 11, such that only children at the end of elementary school (i.e., grade 5 or age 10-11) showed an association between subtraction problems and the left side of space. We also found in Experiment 2 that at least some children may show an association between small addition problems and the right side of space as early as grade 4 (i.e., age 9-10), and that this rightward bias tends to increase at least until the end of elementary school. We can see at least 3 potential explanations for the emergence of these associations.

First, because operands were kept constant across operations, results were overall smaller for subtraction than for addition problems. It could then be argued that number-space associations might have contributed to the difference in spatial bias between subtraction and addition. However, this is unlikely because number-space associations have been demonstrated before 3rd grade in previous studies [START_REF] Hoffmann | Developing number-space associations: SNARC effects using a color discrimination task in 5-year-olds[END_REF][START_REF] Yang | Development of spatial representation of numbers: A study of the SNARC effect in Chinese children[END_REF] and the association between subtraction and the left side of space was only observed in 5th grade in Experiment 1.

Furthermore, no spatial association was observed for multiplication problems in Experiment 2, despite the fact that multiplication problems lead to results that are even higher than addition problems.

Second, it is possible that with practice children progressively associate the "-" and the "+" signs with simple heuristics such as "the result of a subtraction should always be smaller than the first operand" and "the result of an addition should always be larger than the first operand". [START_REF] Mccrink | Operational momentum in large-number addition and subtraction by 9-month-olds[END_REF][START_REF] Hartmann | Spatial biases during mental arithmetic: Evidence from eye movements on a blank screen[END_REF][START_REF] Marghetis | Doing arithmetic by hand: Hand movements during exact arithmetic reveal systematic, dynamic spatial processing[END_REF] Simple mappings between arithmetic signs and space could then explain an operation-dependent spatial bias. This account, however, would also predict an association between the multiplication sign and the right side of space given that results of multiplication problems are larger than the first operand (at least when it is greater than 1). Yet, there was no such effect in Experiment 2.

Third, a number of recent studies have suggested that the repeated use of counting when young children solve basic addition and subtraction problems might lead to an automatization of these counting procedures, rather than to the construction of a network of arithmetic facts in memory as posited by prior literature. [START_REF] Fayol | The use of procedural knowledge in simple addition and subtraction problems[END_REF][START_REF] Barrouillet | On the problem-size effect in small additions: Can we really discard any counting-based account?[END_REF][START_REF] Thevenot | Ten-year-old children strategies in mental addition: A counting model account[END_REF][START_REF] Uittenhove | Fast automated counting procedures in addition problem solving: When are they used and why are they mistaken for retrieval?[END_REF] For example, studies in adults [START_REF] Barrouillet | On the problem-size effect in small additions: Can we really discard any counting-based account?[END_REF][START_REF] Uittenhove | Fast automated counting procedures in addition problem solving: When are they used and why are they mistaken for retrieval?[END_REF] and 10-year-old children [START_REF] Thevenot | Ten-year-old children strategies in mental addition: A counting model account[END_REF] show that the time participants take to solve very small addition problems is not constant but increases linearly as a function of the distance between the original operand and the sum. This suggests that, even if they might not be aware of it, adults and skilled children might solve these basic problems by rapidly "moving" from a source to a target number along the MNL. Given the left-to-right orientation of that MNL in children and adults [START_REF] Fischer | Spatial biases in mental arithmetic[END_REF][START_REF] Hoffmann | Developing number-space associations: SNARC effects using a color discrimination task in 5-year-olds[END_REF][START_REF] Opfer | Early development of spatial-numeric associations: evidence from spatial and quantitative performance of preschoolers[END_REF] , these forward and backward movements are likely to resemble rightward and leftward shifts of attention. Because the automaticity of these shifts during calculation is likely to rely on children's practice with counting and attentional skills (both of which arguably increase from 7 to 11), such an account may explain our findings. That is, spatial shifts of attention may emerge in elementary school because children automatize movements along the MNL that are specific to subtraction and addition.

It is important to note, however, that our results regarding addition problems are inconsistent. That is, we observed a rightward bias in Experiment 2 but not in Experiment 1.

There are a number of possible explanations for this. However, two (non-exclusive) speculations may come to mind. First, it is possible that, in contrast to children in Experiment 2, a majority of children from Experiment 1 may have retrieved answers of small addition problems from long-term memory. For instance, some children might have learned answers of some addition problems by rote, as the French math curriculum emphasizes that small addition problems can be memorized in tables [START_REF] Blanquer | Enseignement du calcul : un enjeu majeur pour la maîtrise des principaux éléments de mathématiques à l'école primaire[END_REF] . This would make small addition problems somewhat similar to multiplication problems and might explain a lack of spatial association.

Note that this explanation remains consistent with the observation that the same children show a leftward bias with subtraction because subtraction problems are seldom learned by rote in the French curriculum. [START_REF] Blanquer | Enseignement du calcul : un enjeu majeur pour la maîtrise des principaux éléments de mathématiques à l'école primaire[END_REF] Unfortunately, this explanation is speculative because we did not collect information on the way children were taught addition problems in our study. However, variability in the methods used to teach addition between schools and classrooms may account from variability in our results.

A second possibility is that automatizing calculation procedures relying on movements along the MNL is likely to require a great deal of practice with counting and arithmetic calculation. Because this may also critically differ between children from different schools and classrooms, it is possible that children from Experiment 1 may have been exposed to less extensive practice with arithmetic calculation than children from Experiment 2. An examination of arithmetic performance between children from Experiment 1 and 2 suggests that it might be the case. For example, 4th-and 5th-graders in Experiment 2 had a significantly higher score on the Math fluency test than 4th-and 5th-graders in Experiment 1 (t (119) =2.16, p = .017, one-tailed). This raises the possibility that we might have failed to observe a rightward bias in children from Experiment 1 because these children were not as fluent with addition problem solving as children from Experiment 2. Note, however, that this possibility is mitigated by the fact that Math fluency scores were not significantly related to dRTs for small addition problems in either Experiment 1 or Experiment 2. Therefore, children from Experiment 1 and Experiment 2 might also differ with respect to cognitive skills that have been shown to affect the automaticity of counting procedures, such as working memory and processing speed. [START_REF] Fayol | The use of procedural knowledge in simple addition and subtraction problems[END_REF][START_REF] Barrouillet | On the problem-size effect in small additions: Can we really discard any counting-based account?[END_REF][START_REF] Thevenot | Ten-year-old children strategies in mental addition: A counting model account[END_REF][START_REF] Uittenhove | Fast automated counting procedures in addition problem solving: When are they used and why are they mistaken for retrieval?[END_REF] Because these skills were not measured here, future studies are needed to better understand individual variability in the development of attentional shifts during addition problem solving. Nonetheless, our findings provide the first evidence for the emergence of associations between arithmetic processing and space in elementary school children. Specifically, we provide evidence that at least some children already associate addition problems with the right side of space by age 9-10 (i.e., grade 4), whereas subtraction problems are associated with the left side of space by age 10-11 (i.e., grade 5). However, these spatial associations may also depend on the level of arithmetic fluency in children, as well as on the way arithmetic is taught in school. Future studies are needed to investigate if these associations become even more salient in older children, and how they relate to mathematical skills.

  Data collection was conducted in two different time waves. The first wave of data collection (Wave 1) occurred between December 2014 and March 2015, and the second wave of data collection (Wave 2) occured between November and December 2019. Participants in Wave 1 were recruited from three private elementary schools in the area of Lyon in France. They consisted in 63 children aged between 7 and 11 years. Data from 6 children were excluded from the analysis for medical (dyscalculia, deafness) and behavioral (lack of answers for one of the two types of operations) reasons. Thus, data from Wave 1 came from the remaining 57 children (32 Female, Mean age=9.44 years, SD=.81). Participants in Wave 2 were recruited from one private primary school in the area of Lyon (N=29) and via a social media website (N=18). Thus, 47 children aged between 7 and 10 years were in Wave 2 (20 Female, Mean age=8.75 years, SD=.61). Data from three children were further removed because they were outliers (see below). As such, final data came from 101 children (50 female) with ages between 7 and 11 years (Mean age= 9.10, SD=.80). All children were native French speakers. Overall, 47 children were in 3rd grade, 33 children were in 4th grade, and 21 children were in 5th grade. The experiment was performed in accordance with the ethical standards established by the Declaration of Helsinki.

  Participants in Wave 1 were recruited from one private primary school in the area of Lyon, France. They consisted in 38 children aged between 8 and 11 years (21 Female, Mean age=10 years, SD=.65). Participants in Wave 2 were recruited from another private elementary school in the same area (N=22) and via a social media website (N=9). Thus, 31 children aged between 8 and 10 years were in Wave 2 (20 Female, Mean age=9.42 years, SD=.48). Data from one child were further removed because they were outliers (see below). As such, the overall sample came from 68 children (41 female) with ages between 8 and 11 years (Mean age=9.74 years, SD=.65). All children were native French speakers. Overall, 42 children were in 4 th grade and 26 children were in 5 th grade. The experiment was performed in accordance with the ethical standards established by the Declaration of Helsinki.

FiguresFigure 1 .

 1 Figures

  Figure 1. Sequence and timing of a sample trial.

Figure 2 .

 2 Figure 2. Relationship between age and dRTs for small arithmetic problems in Experiment 1. (A) Small addition problems. (B) Small subtraction problems.

Figure 3 .

 3 Figure 3. Relationship between age and dRTs for small arithmetic problems in Experiment 2. (A) Small addition problems. (B) Small multiplication problems.

Table 1 .

 1 Descriptive statistics across all children in Experiment 1.

	Tables
	Variable	Mean (SD)	Range min-max Skewness Kurtosis
	Age (years)	9.10 (.80)	7.5-11	.52	-.61
	WJ Math fluency	43.61 (12.69)	23-82	.94	.69
	dRT small addition problems	-3.46 (260.92)	-651.87-621.14	.03	.47
	dRT small subtraction problems	25.49 (241.25)	-679.18-621.14	.25	.09
	RT small addition problems	1710.19 (382.87)	929.45-2661.13	.13	-.40
	RT small subtraction problems	1839.16 (460.59)	794.59-3195.25	.51	.23

Table 2 .

 2 Regression analyses of dRTs for small addition and subtraction problems in Experiment 1.

	Estimate	SE Estimate	β	t	p (two-tailed)

Note.

DV=Dependent variable; SE= Standard Error.

Table 3 .

 3 Descriptive statistics across all children in Experiment 2.

	Variable	Mean (SD)	Range min-max Skewness Kurtosis
	Age (years)	9.74 (.65)	8.33-11	.10	-.97
	WJ Math fluency (N=67)	51.97 (13.13)	27-86	.72	.06
	dRT Small addition problems	91.47 (173.94)	-316.51-455.38	.29	-.33
	dRT Small multiplication problems	-14.61 (209.89)	-500.64-460.62	.02	.06
	RT small addition problems	1629.33 (388.27)	813.90-2600.42	.47	.10
	RT small multiplication problems	1727.22 (355.62)	900.99-2771.22	.30	.56

Table 4 .

 4 Regression analyses of dRTs for small addition and multiplication problems in Experiment 2.

		Estimate	SE Estimate	β	t	p (two-tailed)
	DV: Addition problems
	Age	4.54	2.72	.20	1.67	.099
	Math fluency	.38	2.26	.03	.17	.87
	RT small addition problems	.07	.08	.15	.90	.37
	F (3, 63) =1.29, p=.29, R 2 =. 06
	DV: Multiplication problems
	Age	1.02	3.51	.04	.29	.77
	Math fluency	-.92	2.62	-.06	-.35	.73
	RT small multiplication problems	.05	.10	.08	-.47	.64
	F (3, 63) =.32, p=.79, R 2 =. 02

Note.

DV=Dependent variable; SE= Standard Error.

Acknowledgments

The concept and design of the present research project was developed by Catherine Thevenot and Jérôme Prado. Data collection for Wave 1 was conducted by Auriane Couderc. Data collection for Wave 2 was conducted by Andrea Díaz-Barriga Yáñez, Léa Longo, Hanna Chesnokova, and Emma Langlois. Interpretation of the data, writing, and editing of the manuscript were conducted by Andrea Díaz-Barriga Yáñez, Catherine Thevenot, and Jérôme Prado. Statistical analyses were conducted by Andrea Díaz-Barriga Yáñez, Annabelle Merchie, and Jérôme Prado. Andrea Díaz -Barriga Yáñez takes responsibility for the integrity of the data analyzed. We would like to thank all of the schools and the children who took part in the experiments.

Competing interests

None