Joël Puibasset
email: puibasset@cnrs-orleans.fr

Adsorption/desorption hysteresis of simple fluids confined in realistic heterogeneous silica mesopores of micrometric length: a new analysis exploiting a Multiscale Monte Carlo approach

Adsorption/desorption isotherms in porous materials are commonly used for characterization.

In order to analyze the data, accurate calculations of fluid adsorption in various complex pore models are required. The reversible, low adsorption portion of the isotherm is generally well described by molecular simulation, since the relevant fluid/substrate interactions are described

at molecular level. This molecular approach is however ineffective in the hysteresis region because the large scale spatial distribution of heterogeneities in the pore network is beyond the computer capabilities. On the other hand, coarse grained approaches are more suited to take into account this porous network complexity at large scale and discuss the hysteresis nature, but the molecular description is lost. In this paper, a multiscale approach is introduced which allows both a molecular description of fluid/substrate interactions, and taking into account the connectivity between the various domains in a porous material. The case of argon confined in heterogeneous tubular silica mesopores (MCM-41 or oxidized porous silicon) is considered. Comparison with the simple independent domain theory shows the strong influence of quenched disorder. It is also shown that the independent pore model significantly overestimates the hysteresis width. The effect of pore ends open at only one or at both ends is addressed.

I. INTRODUCTION

In the last decade, many efforts have been focused on the synthesis and characterization of inorganic mesoporous materials of well-defined geometry, because of numerous potential applications in modern science like phase separation, catalysis, chemical sensors and nanotechnologies. 1,2 Among the numerous mesoporous materials reported so far, the MCM-41 silica molecular sieves 3,4 and porous silicon 5 are probably one of the most widely studied, due to their simple structure and ease in preparation. They are interesting for industrial applications as well as for fundamental studies on adsorption phenomena, hysteresis and criticality [6][7][8][9][10] owing to exceptional properties like cylinder-like pore shapes with a high aspect ratio (pore diameter of few nanometers for a pore length of few micrometers), narrow pore size distribution, negligible interconnections between pores, etc. For these reasons they are considered as model materials suitable for comparison between theory and experiment. See for instance Ref [ 11 ] for a review on theoretical and simulation investigations in cylindrical pores and comparison with experiments. Molecular simulations and nonlocal density functional theory methods developed in cylindrical geometries have been shown to significantly improve the predictions of simple fluids adsorption and position of hysteresis over macroscopic Kelvin equation. [12][13][14][15][16][17] However, there is still room for improvement. Several issues have been proposed. The first one concerns the effect of pore ends which have theoretically proven to have a great influence on hysteresis width. [18][19][20] Another point concerns the possible roughness of the silica walls. It is not yet clear whether the silica surface is rough or not. [21][22][START_REF] Berenguer-Murcia | Studies in Surface Science and Catalysis[END_REF][START_REF] Fenelonov | [END_REF][25] However, the development of Monte Carlo simulations and reconstruction methods allows taking into account this possible roughness in order to improve the MCM-41 and porous silicon models. See for instance the recent work by Coasne et al. 26,27 on pore morphology and realistic reconstruction of MCM-41 pores at atomistic level. The last point concerns the long wavelength characteristics of the pores. It has early been recognized that, despite a quite regular pore shape along micrometer scale, the MCM-41 or porous silicon pore walls are not perfectly straight 9,22,28 which results in modulation of pore size and wall thickness along its axis. As a consequence of these long wavelength heterogeneities, the fluid/substrate interactions are modulated along the pore axis. The importance of such nanometer scale physico-chemical heterogeneities was recognized early in literature, 29,30 but studies focused on slit pore geometry until a recent work [31][32][33] showed the importance of such nanometer scale heterogeneities in cylindrical geometry. However, these studies did not focus on silica pores, and did not take into account their atomic roughness. Furthermore, only a limited number of heterogeneities could be considered due to the limitation in computer capabilities, and extrapolation to a realistic distribution of heterogeneities has to take into account their spatial distribution along the pore, which has never been considered yet in a molecular simulation approach.

The aim of this paper is to propose a multiscale Grand Canonical Monte Carlo (GCMC) simulation study which takes into account these requirements. The substrate is modeled by a tubular pore 1 micrometer long, 3nm in diameter, described at atomic level, made of 140 segments of equal length (roughly 7 nm) with various fluid/wall affinities to model physico/chemical heterogeneities along its axis. The main difficulty in such a study is the four orders of magnitude between the atomistic level description (Å) and the pore length (1µm).

The saturated system is expected to contain roughly 1.5×10 5 adsorbed argons. A multiscale approach is then required, which allows to take into account the interconnections between the 140 domains. Comparison with the independent domain theory of Everett 34 where the 140 domains are treated as independent pores, shows that their spatial interconnection is of greatest importance to understand the hysteresis shape. The paper is organized as follows:

Simulation details concerning the atomic description of each domain is given in section 2, as well as intermolecular potentials and GCMC details. The multiscale approach to perform the adsorption/desorption calculations in the entire 1 micrometer long pore made of 140 domains is then introduced in section 3. It is shown that the calculations may be performed separately in each domain taking into account the influence of the two neighboring ones. The corresponding results are given in section 4. They are finally exploited to generate the full pore adsorption/desorption isotherm, with comparison to published experimental data.

II. SIMULATION DETAILS

In order to perform accurate calculations of the adsorption/desorption properties of argon confined in mesoporous cylindrical MCM-41 or oxidized porous silicon pores, the system is described at the atomistic level. The mesopores are essentially modeled by heterogeneous cylinders of various characteristics, drilled in an initial substrate made of pure silica. The atomic positions of silicon and oxygen species are taken from a silica crystal, the cristobalite, because its density is close to that of the amorphous silica which constitutes the pore walls of the real mesoporous substrate. 22 The fluid/wall interaction intensity is then well reproduced since the density of interacting sites is correct. However, amorphous silica presents natural local density fluctuations which necessarily modulate the effective fluid/wall interaction along the pore. We can also mention the amorphous structure of silica which produces a natural corrugation of the pore surface, and, most importantly, the wall thickness which is expected to have non negligible variations along the pore. These characteristics are taken into account as follows. The pore is drilled in a random direction not parallel to the crystal axis to avoid excessively regular atomic ordering at the pore surface. All silicon species with dangling bonds are removed, which introduces a moderate roughness (the cylindrical shape is preserved). In this context, simulation results show that such smooth corrugation does not introduce drastic changes in adsorption/desorption isotherms. On the other hand, the possible silica density fluctuations and variations of the pore wall thickness are expected to introduce stronger effects. It is not possible to modulate the silica density by simply removing or adding some atoms without producing unrealistic silica structures. Performing simulation runs in order to relax the silica structure was also out of the scope of this paper. It was then decided to adopt the simplest solution which consists of modulating randomly the fluid/substrate interaction intensity along the pore axis around its nominal value to mimic the physicochemical heterogeneities previously mentioned. The modulation is typically around 10 to 20 percent, and limited to 30 percent. This generates an effective heterogeneous pore which is exploited to calculate the adsorption properties of the confined fluid. The pore diameter is taken to be 3 nm, and pore length of roughly The Grand Canonical Monte Carlo (GCMC) algorithm [START_REF] Nicholson | Computer Simulation and the Statistical Mechanics of Adsorption[END_REF][START_REF] Allen | Computer Simulation of Liquids[END_REF] is used to calculate the adsorption/desorption isotherms of confined argon. The chemical potential of Argon is imposed and related to the gas pressure in the reservoir using the ideal gas relation since it is a good approximation at the temperature at which the simulations are performed ). During the simulation runs, the amount adsorbed in the substrate fluctuates. For each point of the isotherms (see below), the calculations are performed until the statistical uncertainty is lower than few percent, taking into account the correlations between successive configurations in a refined statistical analysis. Generally, a total of 10 8

Monte Carlo steps per particle are enough. Figure 1 shows a typical Monte Carlo run performed for a given value of the chemical potential, and starting from an initial configuration given by the previous Monte Carlo run at lower chemical potential (adsorption).

The graph represents the amount of adsorbed argon as a function of the number of blocks of 10 4 Monte Carlo trial steps per adsorbed argon. Such blocks are usually required to generate uncorrelated configurations. As can be seen, the amount adsorbed increases, and seems to "converge" after 500 blocks. The Monte Carlo procedure is then run beyond this point to acquire statistical data. It is generally admitted that another 500 blocks would give average values with uncertainty around 5 percent (500 -1/2 ), provided that the 500 new configurations are independent. In most cases, this hypothesis is however not fulfilled, and the uncertainty is actually much larger. This would not be important for qualitative results. However, as shown later, quantitative calculations are required in this study, and the uncertainty has to be smaller than 5 percent. Thus, the correlation between the data has to be evaluated. Let us consider the amount adsorbed N(n) as a function of the block label n. The relative auto-correlation function

is (<N(n+t) N(n)>-< N(n)> 2 )/(<N(n) 2 >-<N(n)> 2 )
where <> denotes the average, and t is the delay in the number of blocks. The auto-correlation is 1 for no delay (t = 0) and decreases to zero for large delay. It is assumed that two configurations separated by t are almost independent (de-correlated) if the auto-correlation function has decreased to 0.10. In the example shown in Fig. 1, the de-correlation delay is around t = 50 blocks. As a consequence, at least 10 4 blocks of 10 4 Monte Carlo trials per adsorbed argon are required to reach an uncertainty of 7 percent (shown in the inset of Fig. 1).

III. PORE MODEL AND MULTISCALE APPROACH

The pore considered in this study is a cylinder, of diameter 3 nm, and 1.03 micrometer long.

As previously explained, the MCM-41 or porous silicon pores are certainly not perfectly (discrete Gaussian-like distribution). The total number of segments or domains is taken to be 140, and their spatial extension is 7.36 nm. The step function used in this study to modulate the fluid/substrate interaction along the pore axis is shown in Fig. 2a as a function of the segment label. The simplest approach to treat the adsorption problem would be within the context of the independent domain theory of Everett. 34 However, it is shown in this paper that this is a crude approximation, due to obvious mutual influence of the adsorbed fluid between the domains. However, due to computer limitations, it is impossible to treat the adsorption problem in a single simulation run for such a long (one micrometer) pore. A multiscale approach is then proposed where the heterogeneities are treated separately, but taking into account the possible influence of their immediate neighbors. It is proposed to divide the long pore (micrometer scale) in several elemental units, called segments, of length compatible with computer capabilities. For instance, in our pore model, the segments correspond to the different domains (see Fig 2b). For segment length larger than the range of the interactions, the calculations in each segment only require the knowledge of the fluid state in the two neighboring segments. In this perfectly exact approach, the calculations could be performed on a multi-processor computer where each processor focuses on one segment and retrieve information concerning the molecular configurations of neighboring segments from other processors at each simulation step. However such an exact approach is quite demanding and can be greatly simplified by a reasonable approximation: the adsorption properties of an elemental pore are essentially determined by the density of the neighboring fluid, and its thermal fluctuations. The molecular configurations generated in the neighboring segments may then be replaced by any set of molecular configurations sampling the same statistical ensemble. The advantage is that the calculations in the elemental pores can be made separately, while taking into account the influence of the fluid adsorbed in the neighboring segments. In other words, the two neighboring segments are treated as refined boundary conditions for the studied segment (see Fig. 2c). Note that, due to hysteresis phenomena, the fluid state in the neighboring pores is not uniquely defined: several possibilities for the boundaries have then to be considered. However, a maximum of two fluid states are actually possible in our segments (see section The last point to be mentioned is that the first and the last segments of the 1 micrometer long pore merit special attention. Two cases are possible depending on whether the long pore is open at both ends to the gas reservoir or closed at one end. In the first case, the boundary condition at both ends of the long pore is simply the gas state, because the chemical potential in the reservoir is always below the saturation point. In the second case, the boundary condition for the pore is the gas state for the open side. The closed side has to be treated separately as a segment closed at one side.

Let us now suppose that we know the adsorption/desorption isotherms calculated at the molecular level (depending on the boundaries conditions) for the 7 elemental pores (the results are given in next section). The multiscale approach which allows to calculate the adsorption/desorption isotherm for the entire 1 µm long pore is now explained. The initial state is the empty pore at very low pressure. All segments are then almost empty, and the boundary conditions for all domains are G-G. The total amount adsorbed is given by the sum of the elementary adsorption in the domains with the given boundaries. The chemical potential of the fluid is then increased (adsorption), and the amounts adsorbed in the elemental pores increase. The boundaries conditions of the domains are still all gaslike (G-G) until fluid condenses to liquidlike state in one of the segments. The neighboring segments are then reconsidered with the new condition (G-L) since one of their boundaries is now liquidlike. If the fluid state in these segments jumps to the liquid state, the next neighbors are also reconsidered, and so on until the fluid state remains unchanged in all segments and the total amount adsorbed is converged. This corresponds to an avalanche where the condensation in a domain propagates to several domains. The pressure is then increased again, and the same procedure is applied to calculate the new adsorbed amount. Finally, for large chemical potential values, complete saturation is reached (at this point, all segment boundaries are liquidlike except for the segments open to the reservoir for which the gaslike boundary is maintained). The same convergence algorithm is then applied for desorption.

IV. ADSORPTION/DESORPTION ISOTHERMS FOR ELEMENTAL DOMAINS

As previously mentioned, the simplest calculation of the adsorption phenomenon in the heterogeneous pore is to consider the domains as independent. The results given by this simple approach are first given. In this case, the properties of each domain are calculated without taking into account the influence of the neighboring pores, and it is generally assumed that the molecular simulations in each domain or elemental segment should be performed with the usual periodic boundary conditions (PBC). To do so, the boundaries as previously defined are removed, and PBC are imposed in the axial direction. We performed the calculations for the nominal interaction intensity Let us now consider our multiscale approach. In this case, the adsorption/desorption isotherms have to be calculated in each elemental pore with the imposed G-G, G-L or L-L boundary conditions. Figure 5 gives the confined argon results for a domain with the nominal interaction parameter

0 O Ar O Ar - -= ε ε
. The amount adsorbed is calculated by integrating the fluid density within the central segment (between z* = -10.78 and + 10.78, see Fig. 3). As can be seen, the two isotherms obtained with symmetric boundary conditions (G-G and L-L) present a hysteresis, while the dissymmetric boundary conditions (G-L) produce reversible results. In this last case, when the chemical potential reaches the coexistence point between the gaslike and liquidlike states, 14,15,31 the liquid/gas meniscus simply goes through the pore from its initial position imposed by the dense fluid condition at one boundary to the opposite side. Conversely, when decreasing the chemical potential from complete saturation, the same meniscus recedes from the gaslike condition imposed to the fluid on the opposite boundary.

This reversible isotherm is also the "equilibrium" isotherm in the sense that it minimizes the grand potential. 31 On the other hand, for gaslike fluid imposed at both ends (G-G), no meniscus preexists: as a consequence, when the chemical potential is larger than the coexistence value, the new liquidlike state has to nucleate from the gaslike state in the pore (eventually at the wall), which is accompanied by a free energy barrier which delays the capillary condensation. After saturation has been reached, desorption occurs at coexistence without nucleation barrier since liquid/gas menisci are present at both ends due to the gaslike conditions imposed at the boundaries. In the last case where two liquidlike boundaries are imposed (L-L), the condensation occurs at coexistence since liquid/gas menisci preexist at both ends, while desorption is delayed since a nucleation barrier has to be overcome, associated to the nucleation of the new gaslike state in the pore in the absence of preexistent liquid/gas meniscus. This shows that the boundary conditions have large impact on hysteresis loop. However, it is emphasized that boundaries do not change the relative stability of the adsorbed fluid state because the grand potential minima are unchanged: above the coexistence chemical potential, the gaslike (respectively liquidlike) state is always metastable (respectively stable), whereas below, the gaslike (respectively liquidlike) state is always stable (respectively metastable). The boundaries only influence the path followed by the system. In other words, the boundaries affect the nucleation barriers, but not the grand potential values of the local minima. Note that, despite some similarities, these results do not deal specifically with the effect of pore ends. 19 The last point to be mentioned is that, far from the hysteresis loop, the simulation data obtained for the three boundary conditions are identical, as expected since in this case only one stable state exists. This actually shows that the procedure does not introduce artifacts. On the other hand, within the hysteresis loop, where two metastable states exist, the boundary conditions have a large influence.

As a validation of the method, it is proposed to recover the PBC results of given by the fluid state on the other side. The system being homogeneous, the fluid state, and as a consequence the boundary conditions, are identical at both ends. One starts with G-G for low chemical potential, up to capillary condensation (which is delayed beyond the coexistence). At this point, the boundary conditions switch to liquidlike (L-L), and the pore is completely saturated. The chemical potential may then be increased to bulk saturation (around 8 in reduced units), and the amount adsorbed is deduced from the isotherm in the elemental pore with the imposed L-L boundary condition. Upon desorption, the same liquidlike branch given by the L-L boundary condition is described in the reverse direction. Note that according to this, pore emptying is delayed below the coexistence point. When the pore empties, the boundary conditions switch back to (G-G), and the initial gaslike isotherm is recovered. In conclusion, the PBC result is expected to be given by the gaslike branch with G-G condition at boundaries, and the liquidlike branch with L-L at boundaries. The result is given in Fig 6

(circles), as well as the results previously given in Fig. 4 obtained with PBC implemented directly in the code (triangles). As can be seen the agreement is excellent within uncertainties.

Our approach slightly underestimates the adsorbed amount very close to the capillary condensation, and slightly underestimates the chemical potential at which pore emptying occurs. However this is not believed to be significant. For instance, the pore emptying occurs via a stochastic nucleation phenomenon, and, very close to the transition, the probability is such that it may or may not have time to occur within a simulation run. Hence, the positions of the two desorption branches might have been reversed. However, for a slightly lower chemical potential, the probability increases drastically, and transition has time to occur in both simulation runs.

Figure 7 presents the same results (adsorption/desorption isotherms in the elemental pores) for the seven values of the fluid-wall interactions considered in this study. For clarity, the gaslike and liquidlike branches are given in separate figures, and logarithmic scale is used for gaslike branches to enhance low adsorption region. As can be seen, the stronger the interactions the larger the amount adsorbed along the gaslike branch (formation of a monomolecular adsorbed layer). Conversely, for weak interactions, few atoms adsorb on the substrate. The liquidlike branch is also higher for larger fluid/wall intensity. However, the effect is less pronounced for this dense liquidlike state. The last point concerns the position of the vertical adsorption and desorption branches. For stronger fluid/substrate interaction, these positions shift to lower chemical potential, accompanied by a narrowing of the hysteresis.

V. ADSORPTION/DESORPTION ISOTHERMS IN DISORDERED SILICA

MESOPORES

The adsorption/desorption properties of argon confined in the 1 micrometer long silica pore made of 140 elemental segments are now evaluated. The calculations are first performed within the independent domain theory. The isotherms for each domain with PBC are then required. They could be obtained directly in a periodic simulation box. However, we showed previously that these data can also be deduced from the isotherms already calculated with imposed boundary conditions (shown in Fig. 7), which saves computing time. The differences are quite small. The advantage in doing so is that we are using the same set of Monte Carlo data to calculate the isotherms within the independent domain theory and our multiscale approach, and hence we are sure that the possible differences between the results are intrinsic to the models and not due to statistical differences between two Monte Carlo runs. The adsorption/desorption isotherm obtained within the independent domain theory is given in ). The second noticeable point is that the "equilibrium" curve never meets the steep portions of the adsorption or desorption branches. On the other hand, the adsorption occurs at equilibrium just above the lower closure point of the hysteresis (µ* = -9), and desorption occurs at equilibrium just below the highest closure point (µ* around -8). The last point to be mentioned in that in this model the effect of the 1 micrometer pore ends cannot be addressed.

Let us now calculate the same adsorption/desorption isotherms in the framework of our multiscale approach in the pore model previously presented (see Fig. 2a). The results are given in Fig. 8b in the two cases where the 1 micrometer long pore is open at only one (dashed line) or at both ends (dotted line). As can be seen, the adsorption/desorption curves present a hysteresis which is much smaller than for the independent pore model (Fig. 8a). As previously, small steps are visible due to the discrete number of fluid/wall interactions used in this preliminary approach. These steps are expected to disappear in a continuous model. The results show that the adsorption branch is not vertical due to the presence of the heterogeneities. The argon uptake is quite important and smooth between µ = -8.55 and µ = -8.22). On the other hand, the desorption branch is almost vertical (pore emptying occurs essentially between µ = -8.67 and -8.70). The value µ = -8.55 corresponds to the limit of stability of the G-G branch in the most attractive pore. This means that rapid adsorption starts to occur in the heterogeneous pore when the most attractive domains have switched to liquidlike state. This creates a meniscus. Upon chemical potential increase, this meniscus can propagate to neighboring pores according to the gaslike branch limit values obtained with G-L boundaries in the various domains (see Fig. 7). The rapid adsorption stops at µ = -8.22, which corresponds to the coexistence between the gaslike and liquidlike states in the less attractive domain. This confirms the picture of meniscus propagation in the pore. Note that this does not necessarily mean that adsorption occurs at equilibrium, since the meniscus propagation may be delayed by the less attractive regions of the pore. ), but remains trapped by the most attractive region of the pore. As a consequence, the less attractive domains do not desorb. However, when the chemical potential is slow enough, the meniscus propagates freely through all domains, which induces a rapid emptying of the pore. This is reminiscent of a pore blocking effect. Note that in the pore open at both ends, the initial desorption step is slightly enhanced because two menisci start to recede at both ends. The main features are however unchanged. Figure 8b also shows the " equilibrium" isotherm calculated in our interconnected pore model. The grand potential in the 1 micrometer long pore is the sum of the grand potential in the elemental pores (extensive quantity). As a consequence, its global minimum (or global equilibrium) is given by the minimum value for each domain (local equilibrium).

As a consequence, the equilibrium curve shown in Fig. 8b (thin solid line) is actually the same as for the independent domain theory (Fig. 8a). The main difference with the previous case, is that the equilibrium curve meets the steep adsorption branch before the highest closure point, and meets the rapid desorption curve at the lowest closure point. This means that adsorption never occurs at equilibrium, even just above the lowest closure point, and desorption never occurs at equilibrium even just below the highest closure point. This result shows that, despite the presence of menisci in the picture previously given to explain adsorption/desorption phenomena, the system is not at equilibrium. Furthermore, it is frequently admitted that at least one of the adsorption or desorption branches occurs at equilibrium; which can then be used to evaluate the pore size distribution. This study shows that none of the branches is at equilibrium, and hence the standard methods to obtain pore size distributions are poorly justified. It is however shown by Fig. 8b that the adsorption branch finally meets the equilibrium curve above µ = -8.45. At this point, the most attractive domains have switched to the stable liquidlike state, and essentially all remaining domains have at least one neighbor in the liquidlike state. These domains then switch to their equilibrium state, given by the G-L or L-L boundary condition isotherms.

Experimental data exhibit qualitatively quite different results between MCM-41 and porous silicon mesopores. Generally, the hysteresis in MCM-41 is made of parallel adsorption/desorption branches (see for instance the work of Morishige et al. 6,7 ), while porous silicon exhibit strongly dissymmetric branches. 8,9 Furthermore, experimental studies on porous silicon have shown that the hysteresis is quite similar in pores open at only one or at both ends, in contradiction with usual theoretical predictions. 8,9 Inspection of Figs. 8a and8b shows that the independent pore model is in poor agreement with porous silicon data, while the model taking into account the interdependence between domains is in good qualitative agreement. This shows that the H 2 type hysteresis observed in porous silicon is in no way the signature of complex pore interconnections similar to Vycor glass, but, as proposed by Wallacher et al., 9 the consequence of a strong disorder in the porous material (presenting or not interconnections). This picture was first proposed in the pioneering work by Kierlik et al. [START_REF] Kierlik | [END_REF][40][41] . Note that Wallacher et al. 9 propose that the quenched disorder is imposed by some variations of the pore diameter along the pores. This work shows that any variation in the fluid/wall interaction along the pore (variations in pore size or shape, in silica density, chemical nature of the surface, etc.) produces the same effect. On the other hand, the MCM-41 experimental data do not exhibit H 2 type hysteresis. This indicates that the MCM-41 pores are probably weakly heterogeneous along their axis, but present a non-negligible size distribution. Further investigations with a distribution of 1 micrometer pores of various sizes presenting weak heterogeneities are underway to confirm this hypothesis. The last point to be addressed concerns the effect of pore ends on adsorption/desorption hysteresis. It was early recognized that pore ends should have a great influence on hysteresis width. [18][19][20] Homogeneous pores open at both ends are expected to exhibit a hysteresis, while the same pores closed at one end are expected to exhibit reversible isotherms. Experimental data do not confirm these predictions since similar hysteresis are observed in both cases. 8 On the other hand, this experimental fact is well reproduced by our simulations taking into account the pore heterogeneities and their interdependence. This suggests that the quenched disorder is much more relevant than the exact nature of pore ends concerning the hysteresis. Our molecular approach then confirms a conclusion of Wallacher et al. 9 that " […] the large number of minor [diameter] variations along the [pore] is a more effective source for metastable states and hysteresis than the single bottle neck even though the change of the diameter is larger than the [minor] variations […]" .

VI. CONCLUSION

A multiscale approach is proposed to address the issue of adsorption/desorption phenomena in heterogeneous tubular pores. The example of argon adsorption in MCM-41 or oxidized porous silicon pores has been chosen because of the four orders of magnitude between the typical pore length (1 micrometer), the pore size (few nanometers) and the atomic scale. In such a situation, the classical macroscopic approaches are unsatisfactory because macroscopic quantities like surface tension are not valid at nanometer scale, and usual molecular simulations approaches are intractable due to the huge number of particles at full saturation (1.5×10 5 ). This multiscale approach consists of dividing the pore in elemental segments of then the liquidlike meniscus has propagated to the corresponding neighboring domains. The meniscus may propagate further, but finally stabilizes, which give the adsorbed amount. The chemical potential in then increased again, and a new adsorbed amount is calculated. After complete saturation is reached (adsorption branch), the desorption branch may be calculated using the same procedure. This multiscale approach has the great advantage to allow for a molecular description of the system while keeping into account the interconnection between the various domains. To evidence the strong effects of domain interdependence, the isotherms are also calculated in the framework of the independent domain theory. Both models exhibit hysteresis. However, the independent domain theory significantly overestimates the hysteresis width, because each elemental pore has its own nucleation barrier, whereas in a real heterogeneous tubular pore with interdependent domains, a liquid meniscus may appear in the most favorable regions of the pore (lowest nucleation barrier) and then propagate to neighboring pores as soon as the chemical potential allows it, without introducing larger delay. In other words, the presence of neighbors has a significant influence on the energetic barrier between the metastable states. Concerning the hysteresis shapes, it is observed that the interdependence between the domains produces a dissymmetry between the adsorption and desorption branches (type H 2 ), in contradiction with the independent domains prediction. This

shows that the H 2 hysteresis is not the signature of a complex Vycor-like network of interconnections but the consequence of a strong disorder in the porous material. This is in agreement with experimental work on porous silicon. It is also noted that for pores open at both ends, adsorption and desorption are slightly favored since two menisci can propagate in the pore. However, the effect is not so important because of the strong disorder which dominates the effect of pore ends, as confirmed by experimental observations. Further work is in progress for a quantitative comparison of predicted isotherms with experimental data.

  uniform along their axis, and present a distribution of heterogeneities or domains. To generate long wavelength heterogeneities, 230, p (1.1) = p (0.9) = 0.205, p (1.2) = p (0.8) = 0.140 and p (1.3) = p (0.7) = 0.040

  Fig. 2d) where an extra contribution is added to the external potential to maintain the fluid in

  are given in Fig.4. Note that in the transition regions, more points are calculated to reduce the uncertainty. However, for clarity, only the most important ones are shown. As can be seen, adsorption increases with the chemical potential. A hysteresis is observed, with vertical (within uncertainty) adsorption and desorption branches, as expected for homogeneous cylindrical pores. These well-known results also show that a simple fluid confined in cylindrical pores exhibits essentially two states: gaslike and liquidlike, corresponding to the two branches of the isotherm. This is why in the fully interconnected problem, only gaslike and liquidlike boundaries are considered.

  Fig 4 with the results of Fig 5 with imposed boundary conditions. In this case, the neighbor segments are the segment itself (by definition of PBC): as a consequence, the boundary condition on one side is

Fig. 8a (

 8a Fig.8a(solid thick curve). As can be seen, the isotherm presents a large hysteresis between µ* = -9 and -8. The lowest value corresponds to the desorption pressure of the most attractive

  Figure 8b also shows small differences during adsorption between open (dotted line) and closed (dashed line) pores. The adsorption in the closed pore is slightly larger because an additional meniscus appears in the first segment. However, since the next two domains are not so attractive ( meniscus rapidly stops, and does not play a role more important than the others. This is a consequence of the strong disorder in our pore model. Let us now focus on desorption isotherm in the pore closed at one end. Desorption starts (slowly) at µ = -8.57, and rapid pore emptying occurs µ = -8.67. The first value corresponds to the desorption branch in the domain -L boundary conditions, while the second value corresponds to the liquidlike branch limit in the most attractive domain with G-L boundary conditions. This shows that, upon desorption, the meniscus starts to recede from the open pore end (the segment open to the gas reservoir has label 140, with 1

  length compatible with computer capabilities (few nanometers long). For each segment, the adsorption/desorption properties are determined, taking into account the fluid state in the neighboring segments. Hopefully, in the case of cylindrical segments, the fluid may exist in essentially two states: gaslike (G) or liquidlike (L). Four boundary conditions are then possible for each segment, to characterize the neighboring fluid state: G-G, G-L, L-G and L-L. This set of isotherms is then used to calculate the adsorption/desorption isotherms in the 1 micrometer long pore, as follows: At low chemical potential, all pores are almost empty, and the total adsorption is given by the sum of elementary adsorption in each domain with gaslike boundary conditions. While increasing the chemical potential, the amount adsorbed increases, until argon condenses in one of the domains. As a consequence, the boundary conditions for the two neighboring domains change (one of the boundaries switches to liquidlike), which may affect the fluid state in this domain (which may switch to liquidlike). If this happens,

FIGURE CAPTIONS FIG. 1 .FIG. 2 .FIG. 3 .FIG. 4 .FIG. 6 .FIG. 7 .FIG. 8 .

 CAPTIONS1234678 FIGURE CAPTIONS

  Figure 1

  Figure 7

ACKNOWLEDGMENTS

Dr. A. Delville, Dr. E. Kierlik, Dr. P. Porion, and Dr. G. Tarjus are gratefully acknowledged for fruitful discussions on adsorption in porous materials. The simulations were performed on workstations purchased thanks to grants from Région Centre (France) and on IBM SP4 supercomputers at the Centre de Ressources Informatique de Haute Normandie (CRIHAN, St Etienne du Rouvray, France) and the Institut de Développement des Ressources en Informatique Scientifique (IDRIS-CNRS, Orsay, France).


Tangui Le Gléau
email: tangui.legleau@orange.com

Xavier Marjou
email: xavier.marjou@orange.com

Tayeb Lemlouma
email: tayeb.lemlouma@irisa.fr

Benoit Radier
email: benoit.radier@orange.com

Towards circular and asymmetric cooperation in a multi-player Graph-based Iterated Prisoner's Dilemma

Keywords: Game theory, non-cooperative games, iterated prisoner's dilemma, tit-for-tat

In collaborations involving multiple actors, it is well known that tensions between individual interest and global welfare can emerge: actors are personally incentivized to have selfish behavior whereas mutual cooperation may provide a better outcome for all. Known as social dilemmas, these cooperation paradigms have aroused renewed interest in solving social issues, particularly in environmental and energy issues. Hybrids methods with Reinforcement Learning (RL) policies and Tit-for-Tat (TFT) strategies have proven successful to identify fruitful collaboration in complex social dilemmas. However, there are also many situations, where cooperation cannot always be given back directly, and has instead to be carried out through one or more intermediary actor(s). This specificity ruins win-win approaches like TFT. To address this specificity, we introduce a Graph-based Iterated Prisoner's Dilemma: a N-player game in which the possible cooperation between players is modeled by a weighted directed graph. In addition to this new paradigm, we propose a graph-based TFT algorithm that we evaluate on multiple scenarios and compare to other algorithms. Our experiments show that leveraging a graph-based structure in the original TFT algorithm allows it to spread favor better collaboration synergies in most situations.

INTRODUCTION

Humanity and industry face numerous issues of cooperation where various self-interested actors intervene. Sometimes, while there is no direct personal incentive to cooperate, mutual cooperation can provide better outcomes for all actors. Such situations involving multiple non-cooperative intelligent agents exist or will exist in various areas where resources are scarce like energy among electricity, road traffic among independent autonomous vehicles, and decentralized multi-agent systems (such as machine learning on private devices). Faced with the proliferation of these intelligent multi-agent situations and faced with social and environmental issues such as a possible upcoming shortage of energy and natural resources, it is increasingly vital to focus more on the studies of learning agents in non-cooperative games [START_REF] Hager | Artificial intelligence for social good[END_REF].

This category of games where players have no incentive to cooperate despite optimal mutual cooperation is known as social dilemmas. For decades, matrix games have been very popular: with atomic actions (cooperation/defection), three kinds of canon-ical games have emerged. They all have in common point at least one incentive to choose defection. It can either be fear (the fear of being exploited by a defector while cooperating) like in trust problems known as the Stag Hunt, or greed (the temptation to exploit a cooperator) in the Chicken Game. At last, when both fear and greed incentivize players to defect, they face the famous Prisoner's Dilemma (PD). Social dilemmas have been well studied for many years. The notion of general-sum non-cooperative games is quite ancient [START_REF] Nash | Non-cooperative games[END_REF] as well as the formulation of the PD [START_REF] Flood | Some experimental games[END_REF]. Therefore, numerous real life situations have been modeled by such dilemmas: from climate change and environmental security [START_REF] Soroos | Global change, environmental security, and the prisoner's dilemma[END_REF] to institutions negotiation with N-person version of the PD as well as tax policies [START_REF] Zheng | The ai economist: Improving equality and productivity with ai-driven tax policies[END_REF].

If the formulation of this problem is rather simple, the issue of efficient behaviors facing such dilemmas has raised a huge interest for decades. In a seminal work, [START_REF] Axelrod | The evolution of cooperation[END_REF] tackles the question of the Iterated Prisoner's Dilemma (IPD) and proposes to oppose various agents in a tournament where they should maximize their rewards playing successively a PD. According to Axelrod, the key of success rests on four properties: a "good" agent should be nice, provokable, forgiving and clear to understand. The winner of IPD tournament was the algorithm Tit-for-Tat (TFT) introduced by Anatol Rapoport [START_REF] Rapoport | Prisoner's dilemma: A study in conflict and cooperation[END_REF]: this simple algorithm (essentially reproducing the behavior of its partner) fulfills the four properties, manages to encourage cooperation and to be robust to defection. Subsequently, variants of TFT [START_REF] Verhoeff | The trader's dilemma: A continuous version of the prisoner's dilemma[END_REF] or alternatives [START_REF] Nowak | A strategy of win-stay, lose-shift that outperforms tit-for-tat in the prisoner's dilemma game[END_REF] have also fueled interest in social dilemmas. Reinforcement Learning (RL) methods have also been used to study the IPD [START_REF] Izquierdo | Reinforcement learning dynamics in social dilemmas[END_REF].

Later, the emergence of Deep Reinforcement Learning [START_REF] Mnih | Human-level control through deep reinforcement learning[END_REF] brought a new interest since the training of complex policies has made possible the study of cooperation in more realistic social dilemmas [START_REF] Leibo | Multiagent reinforcement learning in sequential social dilemmas[END_REF][START_REF] Jaques | Social influence as intrinsic motivation for multi-agent deep reinforcement learning[END_REF]. However, [START_REF] Lerer | Maintaining cooperation in complex social dilemmas using deep reinforcement learning[END_REF] has shown empirically that a Tit-for-Tat strategy remains necessary to help RL policies not to fall in Nash equilibrium (mutual defection), which suggests that studying TFT keeps a huge interest, in particular to complement RL policies in more general non-cooperative games.

Tit-for-Tat has been well studied and performs well. However in some realistic situations, direct cooperation may be less optimal or even impossible. Indeed, sometimes cooperation has to go through one or more indirect cooperator(s). Whereas it exists an optimal multi-agent collaboration, the cooperation is now not necessary bilateral and then the basic TFT strategies are inefficient.

To address and study this asymmetry and possibly enable circular cooperation, we introduce the Graphbased Iterated Prisoner's Dilemma (GIPD), a novel kind of N-player Prisoner's Dilemma with a graphstructure to better generalize the complexity of possible cooperation between players. The key element of our formalism is a weighted directed graph which sets the maximal authorized cooperation for each pair of players. To tackle this new paradigm, we extend the classic TFT to create a novel algorithm: the Graphbased Tit-for-Tat (GTFT). To address the asymmetry and weighted cooperation, our extension adopts a flow network approach to detect and maintain cycles of cooperation within the players.

We simulate our GTFT on dilemmas involving various patterns of circularity, and compare it to some classic TFT. We also conduct some experiments on the choice of the TFT functions involved in GTFT. We study the relevance of the graph-based structure and the choice of parameters with five designed social metrics.

In section 2, we define social dilemmas in particular the continuous PD. Before proposing our ap-proach of GTFT in section 5, we discuss about other formalisms in section 3 for the N-player IPD and introduce our model of GIPD in section 4. At last, methods, results and discussions are presented in sections 6 and 7.

DEFINITIONS OF SOCIAL DILEMMAS

Social dilemmas are multi-player games where at least one of its Nash equilibria is not optimal. The literature describes them as situations where multiple players can choose between cooperation and defection in a game where mutual cooperation is optimal whereas they all are incentivized to defect. In this section, we recall the formalism of those social dilemmas in particular the so-called Prisoner's Dilemma (PD) and its continuous version.

Classic Social Dilemma

Matrix games with two players provide a basic framework to introduce social dilemmas [START_REF] Axelrod | The evolution of cooperation[END_REF]: Two players A and B choose between cooperation and defection, and they receive one of the four possible payoffs given by Table 1: they are denoted in the literature by R (Reward for mutual cooperation), S (Sucker for being exploited by a defector), P (Punishment for mutual defection) and T (Temptation to exploit a cooperator). This matrix game is defined as a social dilemma if the four payoffs verify the following inequalities [START_REF] Macy | Learning dynamics in social dilemmas[END_REF]:

1. R > P: Mutual cooperation is better than mutual defection 2. R > S: Mutual cooperation is better than exploitation 3. At least one of these two inequalities:

• T > R: Greed (3a)

• P > S: Fear (3b) 4. R > 1 2 (S + T ): Mutual cooperation is better than equiprobable different choice (4)

The last condition (4) is relevant in the iterated version of the game to avoid that the rotation Cooperate/Defect becomes optimal. In the literature, there are three kinds of social dilemma depending of which incentive of inequality (3) is verified. In this paper, we focus on the case of Prisoner's Dilemma (where greed (3a) and fear (3b) are verified). The often used payoffs S < P < R < T are given by S = 0, P = 1, R = 3 and T = 5 [START_REF] Axelrod | The evolution of cooperation[END_REF]. This particular game admits one Nash Equilibrium which is (Defection, Defection) whereas the optimal outcome is (Cooperation, Cooperation).

Continuous Prisoner's Dilemma

The classic PD can be extended to a continuous version by extrapolating a payoff with the four values of Table 1. Let us say that two players A and B play a dilemma. Instead of choosing a discrete action between Cooperate and Defect, they now can choose a continuous cooperation degree a ∈ [0, 1] and b ∈ [0, 1] (0 for total defection and 1 for total cooperation). [START_REF] Verhoeff | The trader's dilemma: A continuous version of the prisoner's dilemma[END_REF] proposes to extrapolate the values of S, P, R and T as by the gain function G follows:

G : [0, 1] × [0, 1] → R G(x, y) = xyR + (1 -x)(1 -y)P + x(1 -y)S + (1 -x)yT
(1) We then define the payoffs of a continuous Prisoner's Dilemma between two players A and B choosing respectively actions a and b by V A ,V B = V (a, b), with:

V : [0, 1] × [0, 1] → R × R V (a, b) = G(a, b), G(b, a) (2)
For clarity, we can check that with discrete actions (i.e. a, b ∈ {0, 1} 2 ), we find the four payoffs in the discrete case given in the Table 1:

• V (0, 0) = P, P • V (1, 0) = S, T • V (0, 1) = T, S • V (1, 1) = R, R

RELATED WORK

Social dilemmas, and in particular the Prisoner's Dilemma, have been studied for decades. They have their roots in the Nash notion of non-cooperative games [START_REF] Nash | Non-cooperative games[END_REF] and then have been properly introduced in the Axelrod's tournament of Iterated Prisoner's Dilemmas [START_REF] Axelrod | The evolution of cooperation[END_REF].

The winner of this seminal tournament was the Tit-for-Tat strategy [START_REF] Rapoport | Prisoner's dilemma: A study in conflict and cooperation[END_REF], a rather simple strategy basically consisting in reproducing the opponent choice. Then, numerous improvements have been proposed to the original TFT: [START_REF] Wu | How to cope with noise in the iterated prisoner's dilemma[END_REF] tried to cope with the noise of TFT in IPD, [START_REF] Beaufils | Adaptive behaviour in the classical iterated prisoner's dilemma[END_REF] proposed a memory-based TFT. Other novel strategies similar to the TFT have been introduced e.g. win-stay and lose-shift [START_REF] Nowak | A strategy of win-stay, lose-shift that outperforms tit-for-tat in the prisoner's dilemma game[END_REF].

Reinforcement Learning techniques have also been adopted to train agents to face the IPD: [START_REF] Matignon | Independent reinforcement learners in cooperative markov games: a survey regarding coordination problems[END_REF]] conducted a study on decentralized multi-agent Q-learning in cooperative and non-cooperative games, [START_REF] Matignon | Hysteretic q-learning: an algorithm for decentralized reinforcement learning in cooperative multi-agent teams[END_REF] proposed the hysteretic Q-learning to cope with non-stationarity of non-cooperative agents. More recently, [START_REF] Lin | Online learning in iterated prisoner's dilemma to mimic human behavior[END_REF] conducted a very complete tournament of bandits, RL and TFT agents.

Moreover, variants of Prisoner's Dilemma have emerged: [START_REF] Verhoeff | The trader's dilemma: A continuous version of the prisoner's dilemma[END_REF] introduced the continuous PD and some corresponding adapted strategies of TFT ; [START_REF] Hamburger | N-person prisoner's dilemma[END_REF] proposed a formalism of N-Iterated Prisoner Dilemmas (NIPD) in which payoffs of players are based on the number of cooperators. Some experiments have been conducted on NIPD [START_REF] Yao | Evolutionary stability in the n-person iterated prisoner's dilemma[END_REF], and in particular with Reinforcement Learning [START_REF] Agudo | Reinforcement learning for the n-persons iterated prisoners' dilemma[END_REF]. Some studies also tackled the problem of graphs with Prisoner's Dilemmas [START_REF] Ashlock | Cooperation in prisoner's dilemma on graphs[END_REF], [START_REF] Luo | Prisoner's dilemma in graphs with heterogeneous agents[END_REF]. Their models simulate situations of multiple players placed in a (undirected) graph in which they can play a PD only with their neighbors. Our approach differs greatly in a sense that maximal authorized cooperation between players is given by a directed and weighted graph. Thus, the cooperation between each pair of players is directed and weighted (in particular non-existent) and can enable non-symmetrical and even circular situations.

To finish, recently social dilemmas have known a new interest with the emergence of deep RL: [START_REF] Leibo | Multiagent reinforcement learning in sequential social dilemmas[END_REF] introduced sequential social dilemmas in the form of more realistic games in which RL agents are trained to cooperate between each other. [START_REF] Lerer | Maintaining cooperation in complex social dilemmas using deep reinforcement learning[END_REF] showed empirically that RL agents fall in non-optimal Nash Equilibrium, and that TFT was necessary to help RL policies to mutually cooperate and maintain cooperation. The solution mixing TFT strategies and RL policies they brought is very promising, however their solution doesn't cover asymmetric and circular situations.

GRAPH-BASED ITERATED PRISONER'S DILEMMA

In this section, we propose to extend the N-player Iterated Prisoner's Dilemma with a graph structure.

The main idea is that maximal authorized cooperation within each ordered pair of players is given by a weighted directed graph.

N-player Prisoner's Dilemma

In this section, we describe the model used for the N-player PD, that we define here without graph extension. We consider a decentralized model, i.e. the N players play one continuous PD between each other: each player i chooses a continuous cooperation degree c i j ∈ [0, 1] towards every other player j (c i j = 0 for total defection and c i j = 1 for total cooperation). In other words, each player i chooses a vector of cooperation -→ C i = (c i j ) j∈ 1,N and all decisions of players form the matrix C = (c i j ) ∈ [0, 1] N×N . Note that we adopt the convention that one player cannot cooperate with oneself: ∀i, c ii = 0. After the decision of cooperation degrees, the N(N -1)/2 continuous PD made simultaneously within each pair of players (i, j) give to them the payoffs given by 2.2. Therefore, each player i receives the total payoff:

V i = ∑ j̸ =i c i j c ji R + (1 -c i j )(1 -c ji )P + c i j (1 -c ji )S +(1 -c i j )c ji T
(3)

Graph-based Iterated Prisoner's Dilemma (GIPD)

We introduce our formalism of GIPD that we define by :

• a number N of players

• a weighted directed graph G max of maximal cooperation defined by a weighted adjacency matrix

C max ∈ [0, 1] N×N • a vector D max ∈ (R + ) N of maximal cooperation effort
• a number of steps (or rounds) T max

In this game, at each step t, each agent i can choose N -1 cooperation degrees c t i j with two constraints: first, the maximal value for c t i j is fixed by the graph G max , and secondly its outgoing flow cannot exceed

D i max : C t i + = ∑ j̸ = j c t i j ≤ D max [i].
In practice, the first constraint is verified by clipping the total choice matrix C t = (c t i j ) by the adjacency matrix of maximal cooper-

ation graph G max : C t ← min(C t , C max ). It means for example that if C max [a, b] = 0, player a cannot cooper-
ate with b, so even if a chooses a degree c ab = 1, the effective degree will in fact be equal to 0. To ensure the second constraint, each row -→ C t i is normalized by the factor min(1, D max [i] C t i + ). When the effective matrix C t = (c t i j ) is computed, the agents receive the payoffs determined by the formulation 3. In this game, we assume that all players can observe the full effective choice matrix of the previous step C t-1 , that they have the knowledge of C max and D max , but T max is unknown to them. This formalism of social dilemma allows to model the classic prisoner's dilemma (Figure 1a) as well as what we call circular dilemmas i.e. the presence of a cyclic flow in the weighted directed graph of maximal cooperation G max (e.g. Figures 2c or 2d). For better intuition, we provide in Figure 2 some real examples related to the cooperation graphs represented in Figure 1.

Examples

Although applications of this formalism are out of scope of the paper, we nevertheless provide below (Figure 2) some didactic examples of dilemmas related to the GIPD presented in Figure 1. The edge weights of the maximal cooperation graphs can be roughly viewed in this case as proportional to the payoffs sum of the helper and the "receiver".

The above example of dilemmas represents the numerous situations where multiple actors can share some items (or services) to homogenize their resources (or energy) by mutually giving their surplus to others.

GRAPH-BASED TIT-FOR-TAT

In this section, we update the classic Tit-for-Tat algorithm to adapt it to the GIPD (defined in 4.3). We start by recalling the principle of the classic Tit-for-Tat algorithm, in particular in a continuous version, and then we detail the novel Graph-based Tit-for-Tat algorithm we propose.

Definition of Tit-for-Tat

We assume that a TFT function f k T FT with k steps of memory is a function indicating at each time step t what a player A (facing a player B) should choose as cooperation degree a t ∈ [0, 1] according to the k previous cooperation degrees: 1c) is a 3-player circular PD and (1d) is a 4-player circular PD with an alternative with four 3-player cycles of cooperation (less optimal : flow of 0.5 instead of 1.0). The values of maximal cooperation effort of D max are indicated in gray. (d) Figure 2: A game where N players can decide to cooperate by sharing their items or defect by keeping them. The items have a decreasing marginal utility (e.g. Nk + 1 for the k th same item) thereby conferring the game a dilemma property: the optimal outcome is the mutual cooperation which is exchanging equitably items, but greed and fear lead to the non-optimal mutual defection. According to the topology of resources of agents, circular situations can emerge (2c and 2d).

f k T FT : N × [0, 1] 2k → [0, 1] f k T FT (t, a t-1 , b t-1 , ..., a t-k , b t-k ) = a t (4)
The majority of TFT algorithms needs only the previous step (k = 1) but some approaches needs more steps, like the gradual behavior with memory proposed by [START_REF] Beaufils | Adaptive behaviour in the classical iterated prisoner's dilemma[END_REF] and we can theoretically envisage the case k = +∞ with the use of models learning the history (Recurrent Neural Networks).

For example, the original Tit-for-Tat with discrete choice (cooperation/defection i.e. a cooperation degree in {0, 1}) is defined by cooperating on the first move (cooperation equal to 1) then doing whatever the other player did on the previous move. For better understanding of our notations, we detail the formalism of the classic TFT:

TFT : N × {0, 1} → {0, 1} TFT(t, b t-1 ) = 1, if t = 0, b t-1 , if t > 0.
(5)

Continuous Tit-for-Tat

To address the continuous version of an IPD (defined in 2.2), some continuous versions of Tit-for-Tat exist. We gather propositions of [START_REF] Verhoeff | The trader's dilemma: A continuous version of the prisoner's dilemma[END_REF] and [START_REF] Gléau | Game theory approach in multi-agent resources sharing[END_REF] as well as a contribution (the stochastic incentive) in the following expression:

TFT α,β,γ,r 0 ,c 0 : N × [0, 1] 2 → [0, 1] TFT α,β,γ,r 0 ,c 0 (t, a t-1 , b t-1 ), r t =      c 0 , r 0 if t = 0 αa t-1 + (1 -α)(r t-1 + (1 -r t-1 )b t-1 ), [r t-1 + β(b t-1 -a t-1 )] + + r 0 X γ if t > 0 (
6) In this expression, α is an inertia coefficient to smooth reaction. r t is a coefficient to incentivize cooperation (r 0 is the initial value of r t ). β is the adaptive coefficient to make r t dynamic: it increases (resp. decreases) when partner cooperation increases (resp. decreases). The objective of that adaptive coefficient is to gain safety. At last, X γ is a Bernoulli variable of probability γ, it means that r 0 X γ is equal to r 0 with probability γ and 0 otherwise. This stochastic parameter allows to stimulate incentive in case of null value of r t for all players. At last, for our simulations, we use an null initial degree c 0 = 0.0. To sum up the notations, we consider three kinds of TFT functions:

• TFT alpha: we omit adaptive and stochastic parameters (β = γ = 0), i.e. the incentive coefficient is constant r t = r 0

• TFT beta: we omit stochastic parameter (γ = 0)

• TFT gamma: we use all parameters

Graph-based Tit-for-Tat

We introduce the Graph-based Tit-for-Tat1 . The goal of this algorithm is to extend the classic TFT in particular to cope with asymmetry and the existence of cycles of cooperation. Our idea relies on viewing the cooperation graph as a flow network [START_REF] Ford | Network flow theory[END_REF] and finding the maximum flow which makes a cycle. We detail our algorithm denoted GRAPHTFT for the point of view of player k. The principle is that GRAPHTFT(k) has inner variables that it updates at each step: a cooperation graph C k , a graph which allows to penalize defectors and a source flow D k , an amount of cooperation that it is ready to provide. In the following, we consider that the algorithm uses a TFT function

f T FT : [0, 1] × [0, 1] → [0, 1] such that a t+1 = f T FT (a t , b t

Max flow problem

To find the maximum cyclic flow of cooperation, our algorithm converts its inner cooperation graph into a flow network (Figure 3), and finds the maximal flow.

We use two kinds of algorithms with polynomial complexity (O(∆N 2 ) where ∆ is the discretization number):

• Ford-Fulkerson algorithm [START_REF] Ford | Maximal flow through a network[END_REF]: find the maximal flow (but with shortest path)

Algorithm 1: GRAPHTFT (for agent k)

Input: Max cooperation graph C max and max source flow D max given by the game and a TFT function f T FT

Initialize: C k ← C max , D k ← D max [k] First step: Choose ∀ j ̸ = k, ⃗ C k [ j] ← f T FT (t = 0) for t ∈ [1, T max ] do
for each other agent j do From C t-1 , compute outgoing flow of cooperation of j : (C t-1 j ) + Execute a TFT on j :

c t k j = f T FT (c t-1 k j , (C t-1 j ) + ) Modify the inner cooperation graph: C k [k, j] ← c t k j C max end From C t-1
, compute the incoming flow of cooperation for k : (C t-1 k ) - Update by TFT the next source flow:

D k ← f T FT (D k , (C t-1 k ) -)
Generate a new flow network F from k to k with a source of capacity D k and capacities given by C k (see Figure 3) From F , extract the sub-graph R of maximum flow of cooperation Choose cooperation degrees from max flow: [START_REF] Orlin | A polynomial time primal network simplex algorithm for minimum cost flows[END_REF]: a variant of Ford-Fulkerson, with a minimization of a certain cost per chosen edge (for the cost, we choose the inverse of cooperation so that the flow search prefers longest cycles.)

-→ C t k ← R [k, :] end • Min-cost max flow

METHODS

In this section, we define and detail some methodological points such as our designed metrics as well as the agents and tournaments we use in our simulations.

Metrics

Let N agents denoted ⃗ π = (π i ) acting for the N players of a tournament T during T max steps. We denote V i (T ,⃗ π)(t) the T max payoffs received by the player i and V i (T ,⃗ π) the sum of its T max payoffs. Moreover, we define the social welfare SW (T ,⃗ π,t) as the sum of payoffs of the N agents at time t. In the following, we succinctly define the five metrics we design (adapted from [START_REF] Lerer | Maintaining cooperation in complex social dilemmas using deep reinforcement learning[END_REF]). First, the Utilitarian metric U measures how close the social welfare is from the optimum which is mutual cooperation (C) compared to the mutual defection (D), and the speed Sp measures how fast the Utilitarian metric reaches its maximal value. Then the incentivecompatibility IC is destined to measure the capacity to incentivize cooperation within other players. We define the incentive-compatibility of the agent π by the difference between what one agent (i = 0) receives by cooperating with all other agents π compared to what it would have received by defecting. The safety S f measures the risk an agent takes by preferring π rather than defection face to defectors ; let us note that since defection is dominant (by definition of dilemma), this metric is always non positive: the higher it is, the more the agent is safe. At last, the forgiveness Fg measures how the social welfare is impacted when a "repentant defector" R begins to cooperate only after τ 0 steps. For the sake of clarity, we denote x ⃗ 1 = (x) i the constant joint policy of N or N -1 agents choosing x.

U(T , π,t) = SW (T , π ⃗ 1,t) -SW (T , D ⃗ 1, 0) SW (T ,C ⃗ 1, 0) -SW (T , D ⃗ 1, 0) , Sp(T , π) = 1 τU max τ 0 U(T , π,t)dt IC(T , π) = V 0 (T , (C, π ⃗ 1)) -V 0 (T , (D, π ⃗ 1)) V 0 (T ,C ⃗ 1) -V 0 (T , D ⃗ 1) , S f (T , π) = V 0 (T , (π, D ⃗ 1)) -V 0 (T , D ⃗ 1) V 0 (T , D ⃗ 1) -V 0 (T , (C, D ⃗ 1)) Fg(T , π) = 1 T max -τ 0 T max -1 ∑ t=τ 0 SW (T , (R, π, ..., π),t) -SW (T , D ⃗ 1, τ 0 ) SW (T ,C ⃗ 1, τ 0 ) -SW (T , D ⃗ 1, τ 0 )

Agents and tournaments

In our simulations, we compare several agents of TFT. We use as baseline the vanilla continuous TFT and we compare the different graph algorithms of our novel GRAPHTFT. Moreover, we also compare the choice of the TFT function in our algorithm and evaluate the addition of our improvements (Alpha, Beta, Gamma). For our simulations, we created two kinds of tournaments of N > 2 players with some patterns of circularity. One is purely circular that we denote CIRC(N): the weight of edge (i, j) is equal to 1.0 if j = (i + 1) mod N and 0.0 otherwise. The other one is double circular DOUBLE(N), the same with an alternative co-operation edge (to study the case when one defector breaks the cooperation cycle): the weight of edge (i, j) is equal to 1.0 if j = (i + 1) mod N or j = (i + 2) mod N and 0.0 otherwise.

RESULTS AND DISCUSSION

In this section, we present the results of some simulations where we compare our extension of TFT suited for the GIPD compared to the classic TFT. We first study the relevance of the graph structure of our GTFT and then the impact of the choice of the TFT function. 6) with GRAPHTFT (with min cost approach). We compare the three kinds of TFT (Figure 5a), then the impact of the initial incentive r 0 (Figure 5b) and the adaptive coefficient β (Figure 5c) on the TFT beta algorithm.

Impact of the graph-processing algorithm

First, we focus on the necessity of adding a graphprocessing algorithm and study the impact of the chosen algorithm. The three cases we study are:

• noGraphTFT: classic continuous TFT with no graph processing. • grTFT Fulkerson: our GRAPHTFT (algo 1) with a Ford-Fulkerson algorithm to find the maximum flow. • grTFT minCost: our GRAPHTFT (algo 1) with a min-cost max flow approach: it searches the max flow but with constraint of maximal circulated flow (cooperation through maximum number of players).

We simulate these three kinds of agents in the games DOUBLECIRC with 3, 4 or 6 players and compare the results of the metrics in Figure 4.

1 A simulation of the cycle issue with 3 players: https: //youtu.be/VHhEZ8Wu XQ or with 5 players: https://youtu. be/J-weuvOkkBc 2 A simulation with 4 players: https://youtu.be/ s08j24LMr0U

3 A simulation with 6 players: https://youtu.be/ AL29LFbh3n8

The observations confirm the intuition that the lack of graph processing in circular situations totally defeats the classic TFT in the 6-player game since there is no possibility of symmetrical cooperation. In the 4-player game, the classic TFT can find a less optimal solution since it exists two less optimal 2-player cycles. At last, the classic TFT can address the 3 players DOUBLECIRC because a symmetrical alternative is possible but a bit less optimal. Regarding the choice of graph processing, the best one is the approach with min cost because it allows to select the most pro-social cycles of cooperation (unlike the Ford-Fulkerson algorithm). It is therefore more suited to incentivize cooperation and is safer to defection. However, when it exists several equivalent optimal cycles (e.g. with 3 or 5 players), both GTFT agents naturally hardly synchronize themselves.

Impact of choice of Tit-for-Tat function

We study here the choice of Tit-for-Tat function f T FT used in our GRAPHTFT (algo 1). We first evaluate the impact of the choice of TFT among the described algorithms in section 5.2 (TFT alpha with inertia, TFT beta with adaptive incentive and TFT gamma with stochastic incentive). We then study the impact of initial incentive coefficient r 0 and adaptive coefficient β with the TFT beta.

In Figure 5a, we can observe that the adaptive coefficient β is interesting since it makes TFT beta and TFT gamma safer and more incentive compatible. TFT gamma seems to be a bit more efficient in a sense that the outcome is more optimal. Regarding the impact of parameters of the TFT beta, we can remark that a higher initial incentive coefficient r 0 allows to reach a more optimal outcome more rapidly without reducing safety and incentive-compatibility. Regarding the impact of adaptive coefficient β, we can observe that a tiny non-null coefficient is sufficient to make the agent safer and more incentive compatible.

Results Summary

Although our GTFT is not perfectly optimal in ambiguous situations (with multiple optimal cycles), we can conclude that the key point of our experiments is the importance of a graph-processing algorithm. The min-cost max flow is the best approach since it is more incentive-compatible. Regarding the choice of TFT function, the TFT beta function is clearly safer than TFT alpha while TFT gamma is slightly more efficient.

CONCLUSION

In this paper, we introduced a novel paradigm for the N-player Prisoner's Dilemma where maximal cooperation between agents is induced by a weighted directed graph. This new model is particularly suited to address the asymmetry of cooperation and in particular the circular social dilemmas: a specific situation of dilemma where players can form a cycle of cooperation in which no player can cooperate with its "helper". We showed that classic solutions like Tit-for-Tat strategies cannot solve properly this specific issue, we therefore also proposed in the paper a Graph-based Tit-for-Tat which generalizes the classic TFT with a flow network approach. We evaluated this new algorithm in some scenarios and compare it to some baselines. As major conclusions, we can observe that adding a graph processing in the TFT is relevant since our GTFT outperforms the original TFT in most of situations. As further works, it could be very interesting to address the ambiguous cases with multiple equivalent optimal cycles. We recall our main contributions:

• We introduced and formalized a novel Graphbased Iterated Prisoner's Dilemma: a formalism able to generalize the N-player IPD involving asymmetrical or circular cooperation.

• We designed and formalized several social metrics adapted to this GIPD.

• We constructed a novel Graph-based Tit-for-Tat able to cope with circular cooperation, it is based on continuous TFT and max-flow algorithms.

We are convinced that this new GTFT paradigm which solves circular dilemmas should offer a lot of perspectives particularly in addition to the recent techniques mixing RL and TFT. Finally, in view of the expectations regarding the digital sobriety and the ethical stakes of artificial intelligence, we reiterated the importance of focusing urgently on non-cooperative games, and striving to include this kind of paradigm in the design of our future intelligent systems.

A Theoretical approach

Our conclusions were mainly experimental. In this section, we propose to study some properties of our policy TFT in a cycle of players (CIRC(K)). Proposition 1. (Efficiency) Let a cycle of K players (CIRC(K)) who all follow the policy TFT α,β,0,r 0 ,c 0 .

Then, ∀i ∈ [0, K -1], with j = i + 1 mod K, we have ∀t, c t i j = c t → where c t → is defined by :

c t → = 1 -Q t (1 -c 0 ) avec Q = α + (1 -α)(1 -r 0 ) (7)
Proof. For the K players, the policy of the incoming flow of cooperation is the same as its outgoing flow, so we have the same function of cooperation degree for all players : ∀i ∈ [0, K -1], with j = i + 1 mod K, we have ∀t, c t i j = c t → . First, let us focus on r t . Since β(c t-1 (i-1)i)c t-1 i(i+1) ) = β(c t-1 →c t-1 → ) = 0 and γ = 0.0, we have ∀t, r t = r 0 . We can rewrite the expression of the sequence c t → :

c t → = TFT α,β,0,r 0 ,c 0 (t, c t-1 → , c t-1 → ) = αc t-1 → + (1α)(r 0 + (1r 0 )c t-1 → ) = [α + (1α)(1r 0 )]c t-1 → + (1α)r 0 = Qa t-1 + R with Q = α + (1α)(1r 0 ) and R = (1α)r 0 (8) It is an arithmetico-geometric sequence of common ratio Q and common difference R. If we note W = R 1-Q , the solution of c t → is given by:

c t → = Q t (c 0 -W ) +W (9)
Let us compute the value of W :

W = R 1 -Q = (1 -α)r 0 1 -(α + (1 -α)(1 -r 0 )) = (1 -α)r 0 (1 -α)(1 -(1 -r 0 )) = 1 (10)
Therefore, the final expression of c t → is given by:

c t → = 1 -Q t (1 -c 0 ) with Q = α + (1 -α)(1 -r 0 ) (11)
Proposition 2. (Safety) Let a player A, using the TFT function TFT α,β,0,r 0 ,c 0 to choose a t , he is facing a pure defector B (∀t, b t = 0.0). If β ≥ r 0 (1-α) a 0

, then it exists τ such that ∀t ≥ τ, r t = 0.0 and therefore lim t→+∞ a t = 0 Proof. First, we have :

a t = αa t-1 + (1 -α)(r t + (1 -r t )b t-1 ) = αa t-1 + (1 -α)r t ≥ αa t-1 (12)
With this geometric sequence, we can find a lower bound for a t : ∀t, a t ≥ a 0 α t (13)

Let us find an upper bound for r t :

r t = r t-1 + β(b t-1 -a t-1 ) = r t-1 -βa t-1 ≤ r t-1 -βa 0 α t-1 (14)
Hence, ∀u ≥ 1, r ur u-1 ≤ -βa 0 α u-1 (15)

We can add up the terms:

∀t ≥ 1, t ∑ u=1 r u -r u-1 ≤ -βa 0 t ∑ u=1 α u-1 ⇐⇒ r t -r 0 ≤ -βa 0 1 -α t 1 -α ⇐⇒ r t ≤ r 0 -βa 0 1 -α t 1 -α (16)
Then, the condition β ≥ r 0 (1-α) a 0 is a sufficient condition such that it exists τ for which ∀t ≥ τ, r t = 0.

Finally, ∀t ≥ τ, a t = a τ α (t-τ) and therefore lim t→+∞ a t = 0

  Some examples of Graph-based Prisoner's Dilemmas: (1a) is the classic 2-player PD. (1b) is a 3-player full version of PD. (

2.

  Update the source flow D k with f T FT according to the difference between what k (oneself) "has received" and what k (oneself) "had given" at previous step. 3. Create a flow network F , whose capacities are given by C k with a source vertex directed towards the vertex of k (oneself) with capacity D k and all edges of C k initially directed towards vertex k artificially redirected towards a sink vertex. Thus, we have a flow network allowing to find the maximum cyclic flow (i.e. from k to k). For better understanding, see Figure 3. 4. Compute the maximum flow R on F , i.e. a subgraph in F and extract the next choice of cooperation -→ C k ← R [k, :]

Figure 3 :

 3 Figure3: Steps to find a cyclic sub-graph of max flow from player 1 in the graph (3a). The graph is transformed into a flow network regarding player 1 with a source s and a sink t in step (3b) and the max flow is extracted in (3c).

  Impact of the choice of graph processing in the double circular tournament with 3, 4 and 6 players. The Tit-for-Tat algorithm used for all agents of that simulation is the algorithm TFT beta with parameters (α = 0.6, β = 0.6, r 0 = 0.7) (a) Impact of type (b) Impact of incentive r 0 (c) Impact of adaptive β Figure 5: Impact of the choice of TFT function on DOUBLECIRC(

Table 1

 1

	: Payoffs in a 2-player social dilemma
		Cooperate Defect
	Cooperate	(R, R)	(S, T )
	Defect	(T, S)	(P, P)

For reproducibility, the code has been made available: https://github.com/tlgleo/graph based TFT


Maider Ormaza

Nicolas Bachellier

Marisa N Faraggi

Benjamin Verlhac

Paula Abufager

Philippe Ohresser

Loïc Joly

Michelangelo Romeo

Fabrice Scheurer

Marie-Laure Bocquet

Nicolás Lorente

Laurent Limot

Single-molecule enhanced spin-flip detection

We studied the spin-flip excitations of a double-decker nickelocene molecule (Nc) adsorbed on Cu(100) by means of inelastic tunneling spectroscopy (IETS), X-ray magnetic circular dichroism (XMCD) and density functional theory calculations (DFT). The results show that the molecule preserves its magnetic moment and magnetic anisotropy not only on Cu(100), but also in different metallic environments including the tip apex. Taking advantage of the efficient spin-flip excitation of this molecule, we show how such a molecular functionalized tip boosts the inelastic signal of a surface supported Nc by almost one order of magnitude thanks to a double spin-excitation process.

1

Recent advances in addressing and controlling the spin states of a surface-supported object -atom or molecule-have further accredited the prospect of quantum computing [START_REF] Gatteschi | Molecular Nanomagnets[END_REF][START_REF] Leuenberger | [END_REF] and of an ultimate data-storage capacity [3]. Information encoding requires that the object must possess stable magnetic states, in particular magnetic anisotropy to yield distinct spin-dependent states in the absence of a magnetic field together with long magnetic relaxation times [4][5][6][7]. Scanning probe techniques have shown that inelastic electron tunneling spectroscopy (IETS) within the junction of a scanning tunneling microscope (STM) is a good starting point to study the stability of these spin states. STM-IETS allows for an all-electrical characterization of these states by promoting and detecting spin-flip excitations within the object of interest [8][9][10][11][12][13]. It can also provide an electrical control over these states [14], simplifying the information readout process. As spin excitations need however to be preserved from scattering events with itinerant electrons, single objects are usually placed on non-metallic surfaces such as thin-insulating layers [6,7,15,16] or superconductors [17].

In this sense, new approaches to improve the detection of spin-flip excitations are desirable. With this purpose we present here a novel strategy based on the molecular functionalization of a STM tip. Atomic point-like probe-particles, such as CO, H 2 , D 2 , Xe or CH 4 , decorating scanning probe tips have been proved to enhance the imaging resolution to the extent of resolving the internal structure of large organic molecules, to act as force sensors and signal transducers, to image the local potential energy landscape of an adsorbed molecule, or to improve the detection of molecular vibrations [18][19][20][21][22][23][24][25][26][27]. Tip functionalization with larger organic molecules could introduce additional mechanical, electronical or/and spin degrees of freedom very valuable for further development of the scanning probe techniques capabilities [28].

In this work, we are able to controllably attach a double-decker molecule, nickelocene (Nc hereafter), to the tip apex and exploit the spin-flip excitations exhibited by Nc to boost the signal of an IETS measurement. The chemical structure of nickelocene consists in a Ni atom sandwiched between two cyclopentadienyl rings (C 5 H 5 , Cp hereafter), endowing nickelocene with an electronic structure (a 1g ) 2 (e 2g ) 4 (e 1g ) 2 of spin S = 1 in the gas phase (Fig. 1a) [29].

By attaching a nickelocene to the tip apex of a STM, we show that two spin-flip excitations can be sequentially triggered, resulting in a large enhancement of the inelastic signal as illustrated by the IETS acquired above a surface-supported Nc. This unique situation can be achieved due to the chemical stability, the presence of magnetic anisotropy, the full preservation of the Nc spin on the tip apex, or more generally in a metallic environment, and the extremely efficient spin inelastic effect exhibited by this molecule.

To benchmark the electronic and magnetic properties of a Nc-terminated tip, we investigated first the nickelocene adsorption on a Cu(100) surface using a metal tip. The ring-shaped pattern visible in Fig. 1b indicates that individual nickelocene molecules bind to the surface through a Cp ring, while the other Cp ring is exposed to vacuum (see experimental details in supplementary material [START_REF]Materials and methods are available as supplementary materials[END_REF]). In particular, the Cp ring is centered on a hollow position of Cu(100) most likely adopting the D 5h symmetry (eclipsed configuration) [START_REF] Bachellier | [END_REF]. Figure 1c presents the differential conductance (dI/dV ) versus bias (V ) spectrum The dashed lines correspond to simulations for a S = 1 system with a magnetic anisotropy of D = 3.2 meV; the simulations are scaled to the data by a factor 0.8. Note that the magnetic field is parallel to the X-ray beam.

of an isolated Nc molecule by positioning the tip above the Cp ring, together with its numerical derivative. As shown, the dI/dV spectrum is reduced near zero bias (Fig. 1c) with symmetric steps up to 600% times higher in conductance at an energy of |3.2 ± 0.1| meV (the uncertainty translates the tip-dependency observed). Since the vibrational excitations of nickelocene adsorbed on the surface can be a priori discarded in this energy range [32],

the differential conductance steps reveal the occurrence of spin-flip excitations within the molecule. In view of the spin 1 character predicted upon adsorption on Cu(100) [START_REF] Bachellier | [END_REF], these excitations occur between the ground state (S=1, M=0) and the doubly-degenerated (S=1, M=±1) excited states of the molecule. The threshold energy observed in IETS corresponds then to a direct measurement of the longitudinal magnetic anisotropy energy, D [START_REF]Materials and methods are available as supplementary materials[END_REF]. The large step height observed in the conductance reflects a highly efficient spin-flip inelastic effect which we attribute to the interference of two spin-channels with different symmetries, an effect not seen in previous experiments where only one channel prevails [START_REF]Materials and methods are available as supplementary materials[END_REF].

To further characterize the magnetic origin of the excitation, we performed X-ray magnetic circular dichroism (XMCD) measurements in a magnetic field of 6.5 T on a collection of individual Nc molecules on Cu(100) [START_REF]Materials and methods are available as supplementary materials[END_REF]. The top panel of Fig. 2a shows X-ray adsorption spectra (XAS) recorded at the L 2,3 edges of Ni for left and right circularly polarized light.

The panel below presents the difference between them, the so-called XMCD spectrum. Having a clear dichroic signal is already an indication of the magnetic character of nickelocene on Cu(100). Figure 2b shows the evolution of the signal maximum at the Ni L 3 edge as a function of the temperature for grazing (red dots, θ = 60 • ) and normal (black dots, θ = 0 • ) incidence of the applied magnetic field. At low temperature, the signal shows a clear angular dependence revealing the existence of an easy magnetization direction, more specifically of an easy plane perpendicular to the molecular axis. For the sake of comparison, we have simulated [START_REF]Materials and methods are available as supplementary materials[END_REF] in Fig. 2b the magnetization of a S = 1 system under an external field of 6.5 T and possessing a magnetic anisotropy of D = 3.2 meV. Based on the XMCD results, we can assert that D is positive, implying that M = 0 corresponds to the ground state of nickelocene and the M = ±1 to the excited state.

The magnetic properties of isolated nickelocene turn out to be robust to external perturbations, possibly owing this behavior to the cyclic π* orbital of the Cp rings sandwiching the Ni atom. This differentiates Nc from other systems reported so far [3,7,9,35] in which the local environment is a crucial factor determining the magnetic anisotropy and/or the spin.

At higher coverages, for example, where Nc molecules self-assemble following a T-shaped geometry (Fig. 3b) typical for metallocenes [START_REF] Bachellier | [END_REF]36], no difference was observed in the IETS spectra of upstanding molecules within the layer (Fig. 3a). We have checked the impact of the surface symmetry by changing to Cu(111), also resulting in no difference for isolated or self-assembled nickelocene (Fig. 3a). Supporting these findings, we note that the spin 1 character and magnetic anisotropy are consistent with existing macroscopic measurements of Nc powder samples and crystallites [37,38] or matrix embedded samples [39] where D ranges from 3.17 meV to 4.16 meV; recent theoretical calculations for gas-phase Nc also predict D to be close to 5 meV [40]. In the following, we exploit such robustness to produce a portable source for spin-flip excitations by attaching a single nickelocene to the tip apex of the STM.

Similarly to Xe atoms [41], or physisorbed molecules [42,43], we can transfer a Nc to our monoatomic tip apex by scanning the target molecule at biases close to -1 mV (Figs. 3c and3d). The tip functionalization is reversible, i.e., in the same way the molecule goes to the tip it can be released back to an atom on the surface (+1 mV), indicating that the molecule remains intact during the transfer process. Counter-images of the molecular tip acquired above a Cu atom (inset of Fig. 3d) confirm the presence of a nickelocene molecule. The half-moon shape that is observed can be easily assigned to a tilted Nc [44]. The nickelocene on the tip exhibits similar magnetic properties as previously seen for nickelocene on the surface (Figs. 3a and4a). IETS acquired above the bare surface reveals a larger magnetic anisotropy of (D = 3.7±0.3) meV and conductance steps up to 900% times higher compared to the elastic conductance alone.

The magnetic properties of the molecular tip are computed by fully-relaxing a single nickelocene atop a Cu atom adsorbed on Cu(100). In Figures 3e-g anisotropy energy (MAE) [START_REF]Materials and methods are available as supplementary materials[END_REF]. The MAE is smallest when the molecule is on the surface (Fig. 3f, -1.27 meV) and recovers the gas phase value (Fig. 3e, -1.54 meV) when adsorbed on the Cu atom (Fig. 3g, -1.49 meV); the 17% difference in the MAE is consistent with the experimental changes found for D for Nc on the surface and on the tip. Note that the MAE is essentially affected by the local ligand field as remarked in other systems [45].

The presence of Nc at the tip apex offers the thrilling perspective of an enhanced detection of on-surface spin-flip excitations. To exemplify this effect we acquired IETS spectra above of a Nc molecule adsorbed on the surface using a molecular tip (see Fig. 4a-c). In this case, along with the conductance steps at ∆ = |3.4 ± 0.2| meV, a new set of prominent steps appears at |6.9 ± 0.4| meV, which are the signature of a double spin-excitation process sequentially involving the nickelocene on the tip and on the surface, or vice versa. These inelastic processes can be directly detected in the I(V ) curve shown in Fig. 4d. To explain the origin of these steps, we suppose first, for simplicity, that the molecules on the tip and on the surface have the same threshold energy ∆. When a tunneling electron has an energy below ∆ (region 1) only elastic tunneling processes can take place (Fig. 4e). Above that threshold but below 2∆ (region 2), one of the molecular spins, either the one on the tip or the one on the surface, can be excited (Fig. 4f). Above 2∆ (region 3) a tunneling electron has sufficient energy to excite inelastically both molecules (Fig. 4g). Taking now into account the slightly different threshold energies as evidenced above for Nc on the tip (3.7 meV)

and on the surface (3.2 meV), we find that their sum nicely matches the threshold energy found for a double excitation. In principle, in region 2 two additional steps (at 3.2 meV and 3.7 meV) rather than one (3.4 meV) should be observed, but their detection is beyond our energy resolution.

The most prominent effect associated to the double spin excitation is the drastic enhancement of the conduction step height observed at 2∆ (a factor 43 with respect the elastic signal in the spectrum of Fig. 4b) that we explain as follows. Let us suppose that the conductance is proportional to the probability of taking one electron from one electrode (tip) into the other electrode (surface), in other words that the conductance is proportional to the T -matrix modulus squared [START_REF]Materials and methods are available as supplementary materials[END_REF]. Since the T -matrix has now spin degrees of freedom in the tip and on the sample, we denote the new matrix elements as T n,m where n is the number of excitations of the tip and m the number of excitations of the sample. Then, the conductance will be proportional to n,m |T n,m | 2 , with a sum on the excitations since the channels are supposed independent. In this case, T 0,0 corresponds to the elastic transmission of an electron from the tip to the surface or vice versa (region 1). T 1,0 and T 0,1 correspond to the case in which there is only one inelastic event during the process, i.e., in the molecule on the tip, or on the surface. Finally, T 1,1 is related to the probability of an electron being transmitted inelastically between the electrodes exciting both molecules. Hence, if the applied bias is high enough to open the first excitation channel for the molecules on the tip and on the sample (region 3), we obtain that the conductance is proportional to |T 0,0 is 56 times higher than the elastic conductance alone (region 1).

| 2 + |T 1,0 | 2 + |T 0,1 | 2 + |T 1,1 | 2 .
To summarize, we have shown that a double molecular spin-flip excitation leading to an enhancement of the inelastic signal in spin-flip spectroscopy can be obtained using a Ncterminated tip. The tip-functionalization scheme described can be generalized to study the spin-flip excitations of other systems [10,11], provided that the inelastic energy thresholds are sufficiently large as to be resolved by the given experimental resolution. The enhanced detection will be in any case favored when the inelastic channel transmission is larger than the elastic one. On a more fundamental note, the portability of the spin-flip excitations, whether of nickelocene or eventually of other molecules, could serve as a spin-sensitive probe (with or without magnetic field) and thereby provide quantitative spin-dependent STM data.

This work has been supported by the Agence Nationale de la Recherche (Grant No.

ANR-13-BS10-0016, ANR-11-LABX-0058 NIE, ANR-10-LABX-0026 CSc). Experiments were performed on the DEIMOS beamline [46] at SOLEIL Synchrotron, France (proposal number 20150309). We are grateful to the SOLEIL staff for smoothly running the facility. We thank the national computational center IDRIS, CINES and TGCC (Grant -

[x2015087364]) for CPU time. * ormaza@ipcms.unistra.fr † limot@ipcms.unistra.fr

FIG. 1 .

 1 FIG. 1. Spin-flip excitations of a single Nc molecule. (a) Energy diagram of the zero-field splitting for Nc. The ground state is S=1, M =0 and the excited states correspond to S=1, M =±1. (b) STM topography (14.5 × 14.5 nm 2 , -15 mV, 20 pa) of Nc molecules adsorbed on a Cu(100) surface. Inset: Nc structure calculated by DFT. White, grey and green atoms correspond to H, C and Ni atoms. (c) Differential conductance spectrum and its derivative measured at 2.4 K acquired at 50 pA and -15 mV. The spectrum was measured by means of a lock-in amplifier using a modulation of 150 µV rms at 722 Hz. The d 2 I/dV 2 versus V spectrum was obtained by numerically deriving the dI/dV spectrum.

FIG. 2 .

 2 FIG. 2. XMCD measurements. (a) Absorption spectra for left and right circularly polarized light and XMCD signal at the Ni L 2,3 edges in a 6.5 T magnetic field (60 • off-normal) at 4.7 K for a collection of isolated molecules. (b) Maximum value of the L 3 XMCD signal for normal (grey circle) and 60 • (red circle) off-normal incidence under a 6.5 T field as a function of temperature.

FIG. 3 .

 3 Fig.3f, there is a very small spin polarization that accounts for 2% of the total polarization.From both adsorbed cases, our calculations show that the surface polarizes with a minimum impact on the molecule and the molecular spin is preserved in the three systems studied, obtaining a negligible value for the orbital moment of 0.1µ B . An inspection of the molecular geometry shows a small reduction of the Nc height (3.54 Å) upon adsorption on Cu(100) due to the depopulation of the anti-bonding frontier molecular orbital. This corresponds to a very small (0.1 e -) electron transfer from the molecule to the surface. Interestingly, the charge transfer is reversed when Nc is on the Cu atom; the molecule charge increases by 0.14 e -leading to a larger molecular height (3.62 Å). This translates into a different magnetic

FIG. 4 . 2 and 3 .

 423 FIG. 4. Double spin-flip excitation. (a) dI/dV spectrum of a Nc-terminated tip above Cu(100). Panel (b) and (c) show the dI/dV spectrum of a Nc molecule measured with a Nc-terminated tip and its derivative. Three tunneling regimes are highlighted (see text) indicated by regions 1, 2 and 3. To record the spectra the feedback loop was opened at 100 pA and -40 mV. Note that the different elastic background observed in (a) and (b) is due to a larger tip-surface distance in (b). Panel (d) presents the I(V ) curve showing the opening of the inelastic channels. Panels (e) to (g) present a sketch of the different tunneling mechanisms: (e) there is only elastic tunneling of electrons; (f) the spin of just one of the molecules is excited. In this case the excitation of the molecule on the tip from M = 0 to M = +1 is shown, (g) both molecules undergo spin-flip excitations.

  The observed inelastic versus elastic ratios show tip dependency. Considering the maximum step height values we observed, we have for the Nc on the tip that (|T 1,0 | 2 + |T 0,0 | 2 )/|T 0,0 | 2 ≈ 9 and for the Nc on the surface (|T 0,1 | 2 + |T 0,0 | 2 )/|T 0,0 | 2 ≈ 6. Since the excitations of the molecule on the tip and on the sample are independent, the transmission probability of a double spin-flipexcitation to occur is approximately |T 1,0 | 2 |T 0,1 | 2 /|T 0,0 | 4 = |T 1,1 | 2 /|T 0,0 | 2 ≈ 8 • 5 = 40. conductance is then in region 3, (|T 0,0 | 2 + |T 1,0 | 2 + |T 0,1 | 2 + |T 1,1 | 2 )/|T 0,0 | 2 =1+6+9+40, this