Nicolas Combe

Christophe Taillan

Joseph Morillo

Joseph Morillo Unidimensional

Unidimensional model of the ad-atom diffusion on a substrate submitted to a standing acoustic wave II. Solutions of the ad-atom motion equation

come

Unidimensional model of the ad-atom diffusion on a substrate submitted to a standing acoustic wave II.

Solutions of the ad-atom motion equation

I. INTRODUCTION

The self-organization of materials at the nanoscale is a promising way to avoid the expensive lithography stage in the conception of the semi-conductors devices. Common self-assembling techniques rely on the Stranski-Krastanov growth modes, 1,2 buried dislocations networks, 3,4 or nano-patterned substrates. [5][6][7][8] An alternative approach to self-assemble materials at the nanoscale, the dynamic substrate structuring effect has been recently proposed. 9 In this approach, a Standing Acoustic Wave (StAW) governs the diffusion of the ad-atoms on a substrate. In the first paper of this series, 10 we have established an unidimensional model of an ad-atom diffusing on a substrate submitted to a StAW and derived the ad-atom motion equation, a Generalized Langevin equation:

m d 2 x dt 2 + t t0 γ(x(t), x(t ′ ), t -t ′ ) dx dt ′ (t ′ )dt ′ = - dΦ eff dx (x) + ξ(t) + F SAW (x, t). (1)
The first left hand side (l.h.s.) term of Eq. ( 1) is the usual inertial term. The second l.h.s. is a retarded friction force with memory kernel γ(x(t), x(t ′ ), t -t ′ ) where t 0 corresponds to the time when the StAW production mechanism is switched on: since γ(x(t), x(t ′ ), t -t ′ ) is a decaying function of |t -t ′ |, we can fix t 0 at -∞ in the integral without loss of generality. The first right hand side (r.h.s.) term is an effective inter-atomic (substrateadatom) periodic potential force derived from the effective potential Φ eff (x), the second r.h.s. term ξ(t) is a stochastic force and the last one F SAW (x, t) is an effective force induced by the StAW. 10 The goal of this paper is to study the solutions of Eq. (1) and to evidence the structuring role of the StAW on the adatom diffusion, through the effective F SAW (x, t) force. In Ref. 10 we showed that F SAW (x, t) reads F SAW cos(kx + ϕ) cos(ωt) with k and ω the wave vector and angular frequency of the StAW. However, due to the precise nature of the ad-atom-substrate interactions, F SAW (x, t) experiences also some variation at the substrate lattice parameter scale: the proportionality factor F SAW and the phase factor ϕ vary as a function of the exact position of the ad-atom in between two successive atomic substrate potential wells. Because we focus here on the structuring role of the StAW, we will only consider its large scale variation, i.e. F SAW (x, t) = F SAW sin(kx) cos(ωt) with a constant F SAW value and a constant phase factor that we fix at -π/2 for convenience, eluding thus the possibility for this force to vary on the substrate lattice parameter scale: such dependence essentially affects the dynamics of the ad-atom and weakly the structuring effects of the StAW. 9 The detailed study of the dynamics of the ad-atom will be reported elsewhere.

Concerning the stochastic force ξ(t) and the memory kernel γ(x, x ′ , t-t ′ ), which are related by the fluctuationdissipation theorem: 10,11 ξ(t)ξ(t + τ ) = k B T γ(x(t), x(t + τ ), τ ), (2) with k B the Boltzmann constant and T the temperature, their properties are very interaction-potential model dependant and thus cannot be studied on a general ground.

Though we have established their analytical expressions in a very peculiar potential case in the preceeding paper, 10 we will here use the more standard and general model of a centered gaussian noise for the stochastic force, with correlation time τ c . [START_REF] Pottier | Physique statistique hors d'équilibre[END_REF][START_REF] Gordon | [END_REF][START_REF] Forster | Hydrodynamic fluctuations, broken symmetry, and correlation functions[END_REF] • The stochastic force is then fully characterized by its zero mean value and autocorrelation function:

ξ(t) = 0, (3a)
ξ(t)ξ(t + τ ) = D e -|τ |/τc τ c ; (3b)
• and the memory kernel is given by (Eq. ( 2)):

γ(x(t), x(t + τ ), τ ) = γ τ c e -|τ | τc , (4)
where D/m 2 is the adatom diffusion coefficient in the velocity space in the absence of the StAW and γ is the friction coefficient, given by γ = k B T /D. Eqs. (3) relies on the high number of atoms in the substrate and its thermalization by an external macroscopic system in experimental conditions. 10,[START_REF] Pottier | Physique statistique hors d'équilibre[END_REF] In this model, the correlation time τ c is of the order of the inverse of the Debye frequency of the substrate. The meaning of the friction coefficient appears clearly in the limit of vanishing correlation time :

ξ(t)ξ(t + τ ) τc→0 ---→ 2Dδ(τ ), ( 5a
)
where δ is the Dirac distribution, and thus

t -∞ γ(x(t), x(t), t -t ′ ) dx dt ′ (t ′ )dt ′ τc→0 ---→ γ dx dt (t), ( 5b
)
the usual friction force proportional to the speed. In this limit, m/γ is the relaxation time τ R of the ad-atom dynamics in the absence of r.h.s. in Eq. (1). Eq. ( 1) involves two different length-scales, the StAW wave-length λ = 2π/k and the periodicity a of the crystalline potential; and three different time-scales, the correlation time τ c , the StAW period T = 2π/ω and the relaxation time τ R of the ad-atom dynamics. In the following analytical calculations, we assume that all these scales are separable and have the following properties:

a ≪ λ (6a) τ c ≪ T ≪ τ R (6b)
Eq. (6a) derives from the fact that we consider StAW wavelengths varying from few to hundreds of nanometers compared to the substrate lattice parameter of the order of 10 -1 nm. Eq. (6b) is motivated by ours Molecular Dynamic (MD) simulations of ad-atom diffusion on a substrate submitted to a nanometer wavelength StAW. 9 The effective crystalline potential (Φ eff ) will be modelled by its fundamental Fourier component:

Φ eff (x) = - aϕ 0 2π cos(2πx/a), (7)
where ϕ 0 is the amplitude of the corresponding force.

Because of the difference of length scales (Eq. (6a)), we assume for simplicity that the wavelength of the StAW force is an integer multiple of the lattice parameter λ = na.

Φ eff (x) = - ϕ 0 nk cos(nkx). ( 8
)
Within the presented model, Eq. ( 1) writes now: (9) where ξ is the Gaussian noise defined by Eqs. (3).

m d 2 x dt 2 + γ τ c t -∞ e -|t-t ′ | τc dx dt ′ (t ′ )dt ′ = -ϕ 0 sin(nkx) + ξ(t) + F SAW sin(kx) cos(ωt),
In a reduced dimensionless formulation Eq. ( 9) writes:

ẍ + γ τc t -∞ e -|t-t ′ | τc ẋd t′ = -φ0 sin(nx) + ξ( t) + F sin(x) cos( t), (10)
with x = kx, t = ωt, γ = γ/(mω), τc = ωτ c , ξ( t) = kξ(t)/(mω 2 ) and F = kF SAW /(mω 2 ) and where ẋ and ẍ are the first and second derivatives of x relative to t respectively.

Getting an explicit expression of the general solution of the non-linear stochastic integro-differential equation ( 9) or ( 10) is unreachable at least for us and, anyways is not our goal. As stated before, what we are interested in, is to evidence the conditions for the appearance of a structuring effect on the ad-atom diffusion due to the effective StAW force and how the other forces can affect it.

In Sect. II, we evidence the structuring effect of the StAW force, by considering Eq. ( 10) in the long time (t ≫ τ c ) and length scales (x ≫ a) limit without thermal fluctuations. In this limit, the periodic potential and the stochastic forces can be neglected and the retarded effects in the friction force too (limit τ c → 0). Eq. ( 10) reduces then to:

ẍ + γ ẋ = F sin(x) cos( t). ( 11
)
Since Eq. ( 11) is non-linear, we apply, in Sect. II some standard tools of the non-linear physics to characterize its solutions. First, using the multiple time scale analysis, [START_REF] Kevorkian | Multiple Scale and Singular Perturbation Methods[END_REF] we will evidence the existence of an effective potential U eff governing the diffusion in the long time scale limit.

The general solutions will then be studied using the fixed points stability analysis. Then, it will be numerically integrated and its solutions will be analyzed with the use of the Poincaré sections of the phase diagram and the calculation of their Lyapounov exponents. The stochastic force ξ(t) will then be reintroduced in Sect. III in the same τ c → 0 limit. It will be shown that it mainly induces some fluctuations around the solutions of Eq. (11). In Sect. IV, the retarded effects (τ c = 0) will be added, whereas the effects of the substrate effective cristalline potential will be reintroduced in Sect. V in the absence of retarded effects (limit τ c → 0). In all the sections from III to V , the structuring effect of the StAW force and its sensitivity to the other forces will be evidenced through the study of the position probability density of the ad-atom. Finally, in Sect. VI the relevance of the analysis of Sect. II on the solutions of the complete equation (Eq. ( 10)) will be demonstrated.

II. MAIN EFFECTS OF THE STANDING ACOUSTIC WAVE ON AD-ATOM DIFFUSION

We study in this section the non-linear equation Eq. (11). Note that this equation is invariant under a simultaneous space and time translation: x → x + π and t → t + π. As a consequence, all the results concerning the structuring effects will be invariant under a space translation x → x + π.

A. Multiple time scale analysis

Eq. ( 11) is a non-linear deterministic equation that implies two different typical time scales: the dimensionless relaxation time τR = ωτ R (depending on γ) of the adatom dynamics and the period of the StAW force. To manage these time scales, this equation is first analyzed using the multiple scale method. [START_REF] Kevorkian | Multiple Scale and Singular Perturbation Methods[END_REF] This method, due to the extra degrees of freedom it introduces, allows to remove the secular divergencies that can arise in a standard perturbation approach. Note that Eq. ( 11) without friction has already been studied [START_REF] Escande | [END_REF] in a different framework.

We consider the limit ǫ = 1/ω → 0 (keeping constant c s = ω/k, the sound speed in the substrate). Eq. ( 11) reads then:

ẍ + ǫγ ẋ = ǫF sin(x) cos( t), ( 12
)
where γ = γ/m and F = F SAW /(mc s ) are order 0 quantities (O(ǫ 0 )) . We seek an approximate solution of Eq. ( 12) using the following expansion:

x( t, ǫ) = x0 ( t0 , t1 ) + ǫx 1 ( t0 , t1 ) +ǫ 2 x2 ( t0 , t1 ) + O(ǫ 3 ), (13)
involving the two time scales: t0 = t and t1 = ǫ t. Substituting Eq. ( 13) into Eq. ( 12) and identifying terms of the same order in ǫ, we get the following equations:

D 2 0 x0 = 0, (14a)
D 2 0 x1 = -2D 0 D 1 x0 -γD 0 x0 + F cos( t0 ) sin(x 0 ), (14b)
D 2 0 x2 = -2D 0 D 1 x1 -D 2 1 x0 -γ(D 0 x1 + D 1 x0 ) +F x1 cos( t0 ) cos(x 0 ), ( 14c
)
where the operator D n designs ∂/∂ tn with n ∈ {0, 1}. The solution of Eq. (14a) reads:

x0 = A( t1 ) t0 + B( t1 ). ( 15
)
The first term of Eq. ( 15), a secular term that diverges with t0 , is removed by setting A( t1 ) = 0, leading to:

x0 = x0 ( t1 ). ( 16
)
Substituting Eq. ( 16) into Eq. (14b), the particular solution x1 writes:

x1 ( t0 , t1 ) = -F cos( t0 ) sin(x 0 ). ( 17
)
Using these expressions of x0 (Eq. ( 16)) and x1 (Eq. ( 17)) in Eq. (14c), we write the solubility condition of this equation by the elimination of the secular term:

D 2 1 x0 + γD 1 x0 = - F 2 2 cos(x 0 ) sin(x 0 ), (18)
a differential equation governing the solution x0 on the time scale t1 . To give a physical meaning to Eq. ( 18), we note that using Eqs. ( 16) and (17) in Eq. ( 13), the solution of Eq. ( 11) writes to O(ǫ 1 ) order:

x( t) = x0 ( t1 ) -F cos( t0 ) sin(x 0 ( t1 )).

Therefore, X( t1 ) the average value of x( t) over a StAW period, writes to O(ǫ 1 ) order:

X( t1 ) = x( t) = 1 2π t0+2π t0 x( t)d t ≈ x0 ( t1 ). ( 20
)
Substituting x0 ( t1 ) by X( t1 ) in Eq. ( 18) an going back to the t variable, we obtain the long time evolution equation of X( t) = x( t) , the mean value of x over a StAW period:

d 2 X( t) d t2 + γ dX( t) d t = - F 2 4 sin(2X( t)). ( 21
)
Note that the scheme provided by Landau and Lifshift 17 that develops the variable x( t) in Eq. ( 11) as a sum of a slowly varying function X( t) and a quickly varying function ζ( t) yields a similar result. 9

The multiple scale method has allowed us to transform the non-autonomous Eq. ( 11) into an autonomous equation (Eq. ( 21)) on a longer time scale. This equation describes the motion of the ad-atom on a mesoscopic time scale, long compared to the period of the StAW, but small compared to the relaxation time τR of the ad-atom dynamics. The StAW force acting at the mesoscopic timescale derives from the effective potential U eff :

U eff (X) = F 2 4 sin 2 (X). ( 22
)
This potential, is periodic with period π and minima at X = 0[π]. Whatever the initial conditions the long time evolution described by Eq. ( 21) will be a dampted evolution towards one of the minima of U eff . Within the described approximation (Eq. ( 21)), the StAW leads then to a self-organization of the ad-atoms diffusion into a periodic array with period half the period of the StAW. In the following, U eff will appear to be an essential tool to interpret the ad-atom trajectories and the structuring effect of the StAW. Note that the (stable) fixed points of Eq. ( 21) which are the minima of U eff , are also those of Eq. (11). These results are approximated results, we must now come back to Eq. (11) to test their relevance in the general case. We will start with the study of the stability of the fixed points ((x, ẋ) = (0[π], 0)) of Eq. ( 11) which, as we will see now, can differ from that of Eq. ( 21) for certain values of the parameters of the equation.

B. Fixed points stability

As already mentionned, Eq. ( 11) is invariant under simple spatial and time translations so that all the fixed FIG. 1. Stability (Blank) (SD) and instability (dashed) (ID) domains of the fixed point (0,0) of Eq. ( 11) in the (γ, F ) parameter plane (PP). Red crosses design the couples of parameters chosen in the numerical studies of Sects. II, III and IV.

points (x, ẋ) = (0[π], 0) are equivalent. We hence reduce our stability analysis to one fixed point: (x, ẋ) = (0, 0). The linerarized version of Eq. ( 11) around this fixed point reads:

V = A( t)V ( t) (23)
with V = x ẋ and A( t) = 0 1 F cos( t) -γ .

We define R t t0 , the propagator of Eq. ( 23):

V ( t0 ) R t t0 -→ V ( t) (24)
Since A( t) = A( t + 2π) in Eq. ( 23), the Floquet theory provides the stability of the fixed point (x, ẋ) = (0, 0) of Eq. ( 11) from the eigenvalues of R 2π 0 : the fixed point is stable if all the eigenvalues of R 2π 0 are inside the unit circle of the complex plane.

The propagator R 2π 0 is calculated by numerically integrating Eq. ( 23) using a fourth order Runge-Kutta method and is then diagonalized. Fig. 1 reports the stability diagram of the fixed point (x, ẋ) = (0, 0) in the (γ, F ) parameter plane (PP). Fig. 1 reveals the existence of unstable domains which questions the validity of the self-organization effect evidenced for any (γ, F ) parameters in the preceding section. This is not contradictory since the Floquet analysis is exact concerning the fixed points stability, while the approximated multiple scale analysis relies on the smallness of the amplitude of the fast variations (StAW period time scale) compared to the slow variations (τ R time scale) of x. A condition which will be fulfilled as far as the strength of the friction force (γ) remains large enough compared to the strength of the effective StAW force ( F ), which explains the global separation between the stability (upper left triangle) domain (SD) and the instability (lower right triangle) domain (ID) in the PP. In addition, in absence of friction γ = 0, Eq. ( 23) is equivalent to a Mathieu equation for a parametric oscillator with a null eigen frequency. The stability diagram of the Mathieu equation can be for instance found in Fig. 2 of Ref. 18: it evidences some unstable solution regions, corresponding to the parametric resonances, separated by very tiny stable regions for an oscillator with a null eigen frequency. The large separated unstable subdomains (tongue) in Fig. 1 hence correspond to the parametric resonances that disappear when increasing the dissipation γ.

C. Numerical study

In order to check the self-organisation ability of the effective StAW force in the ID of the PP (Fig. 1), we need to turn to a precise examination of the trajectories.

An analytical resolution of Eq. ( 11) is out of scope, so we will solve it numerically using the forth order Runge Kutta method for different values of the parameters (γ, F ). We will also look at the trajectories for parameters in the SD to give a complete scope of the different exact behaviors. In order to enlighten some specificities of the trajectories, the resolution will be performed for different initial conditions (x(0), ẋ(0)) (IC). From our analysis of the stability diagram in the preceding section, we expect to observe specific behaviors associated to increasing values of F at constant γ, going from one stability domain to the next one through the midway instability domain. This is what we observed in our rather extended exploration of the PP and IC spaces. From this exploration, we identified four categories of trajectories: one corresponding to the SD; the three others to the ID, one in its core and the two others in the vicinity of its frontiers with its two neighboring SDs. We did not examine the very peculiar case γ = 0, which corresponds to undamped trajectories, since it is not relevant for the ad-atom diffusion on a substrate, and since it has been studied previously. [START_REF] Escande | [END_REF] Of course, we cannot absolutely exclude the possibility to have missed some specific behaviors, even if we consider it as highly improbable. In all the rest of the manuscript, we will focus on the solutions for the 2.0 constant γ value, since it provides a representative sample of the behaviors we have exhibited. The results are presented for increasing values of F (1.5, 3.15, 4.5 and 4.96) through the first instability subdomain (red crosses in Fig. 1). In each case, the trajectory is related to the effective potential U eff (X) ( blue solid line in Fig. 2 and Eq. ( 22)) in order to evidence any self-organization behavior.

Converging Trajectories

The trajectories have typically an oscillating behavior whose amplitude and off-set are decreasing. They can be somewhat different at the very beginning, depending on the sign of the x(0) ẋ(0) product but they have the same character at a longer timescale. Fig. 2a reports the solution for initial conditions (1.0, 0.0) corresponding to an initial position in the U eff (X) potential valley associated to its 0.0 minimum. We choose an initial position (1.0) rather away from the minimum of U eff (X) to best evidence the amplitude and off-set decreases. The initial speed has been fixed to zero in order to only present the characteristic time behavior.

The trajectories converge towards the fixed point (0, 0). Depending on the initial conditions, the trajectory can eventually escape from the 0.0 minimum to an adjacent minima. In that case the trajectory will converge to the corresponding U eff (X) minima (fixed point). The several time scales mentioned in Sect. II A are clearly visible: fast oscillations at a 2π period (StAW period) whose amplitude slowly decays on the τR time scale. The period of the oscillations can be one of the harmonics of the StAW period. For example, in the vicinity of the stability/instability frontier (γ, F ) = (2.0, 3.028), we observe the second harmonic period (4π). It appears then that, in the SD, the solutions of Eq. ( 11) can be adequately estimated using a multiple time scale analysis, the long time behavior being well described by Eq. ( 21).

Periodic trajectories

In the ID, in the vicinity of the first frontier between the SD and the ID, the unstable character of the fixed point results in an amplification of the oscillations (Fig. 2b). This amplification occurs up to an upper limit fixed by the non-linearities of the sinus in Eq. (11). At this time the trajectory becomes periodic, revealing the existence of a limit cycle, the attractor of the system in the phase space. The IC of the presented trajectory (x(0), ẋ(0)) = (0.1, 0.0) correspond to a position very close to the unstable fixed point to evidence the amplification of the fast oscillations (transient regime) leading to the periodic trajectory. The period of oscillation is 4π, i.e. twice the excitation period. Like in the converging trajectories, this period depends on the peculiar choice of the parameters (see Sect. II D).

The trajectory of Fig. 2b oscillates in its initial effective potential basin. However, increasing F , results in an increase of the upper limit. This limit may overcome the maxima of U eff (X) so that the trajectory may extend on the two neighboring effective potential basins: however, the average position of the adatom always belong to the same potential basin. U eff (X), if relevant to characterize the trajectory of the adatom, does not capture all the information contained in the phase space diagram. With this kind of trajectories the self-organisation character of the effective StAW force is not apparent but will be demonstrated in Sect. III studying the distribution of the ad-atom x-position.

Chaotic trajectories

In order to evidence the chaotic character of some trajectories in the core of the ID region ( F = 4.5) we report on Fig. 2c the solutions of Eq. ( 11) for two sets of very closed I.C.s (0.10, 0.0) and (0.12, 0.0). Both trajectories visit different effective potential basins following an apparently erratic motion. The two solutions, diverge quickly from each other, despite their very close ICs, suggesting the chaotic character of these trajectories. To examine this chaotic behavior, we calculate the Lyapounov coefficients (λ i ) of the autonomous system associated to Eq. (11). For the first trajectory, using the algorithm provided by Wolf et al, 19 the highest Lyapounov coefficient λ L is found to be positive (0.278) (binary base), an unambiguous evidence of the chaotic character of this solution.

The transition to the chaotic behavior will be discussed in a separate section (Sect. II D). At first sight, it seems that the self-organization character is lost, but as we will see later (Sect. III) studying the position probability density of the ad-atom, it actually is preserved.

Unbounded trajectories

We look now, in the ID domain, at a characteristic trajectory in the vicinity of the second frontier between the ID and the SD ( F = 4.96). The IC of the trajectory reported in Fig. 2d are the same as those in the two preceding cases, (0.1, 0.0). After a transient period where the particle stays in its original potential well, it leaves it, without being captured by any other potential well: the trajectory is unbounded. As one can see, from this point the trajectory is roughly linear and thus does not present any visible chaotic character. This is confirmed by the calculation of its Lyapounov exponents that are all negative (or null). After the transient period, the trajectory is monotonously increasing with a staircase character. Changing the IC can lead to a monotonously decreasing trajectory with a symmetric staircase character. These trajectories are analogous to the rotations of a pendulum about its pivot in the clockwise or counterclockwise directions and have been already evidenced in the phase space (x(0), ẋ(0)) in the absence of dissipation γ = 0. [START_REF] Escande | [END_REF] The final remark of the preceding section still applies here: contrary to the appearances, the self-organization character is preserved. It infers from the staircase character of the trajectory which corresponds to longer residence times in the potential valleys than the transition times between valleys. The self-organization character will be evidenced for all the types of trajectories in Sect.III, where the probability distribution will be studied in the presence of the fluctuating force (vanishing fluctuating force case).

D. Transition to Chaos

We will now characterize the domains of existence, in the PP plane, of the four observed types of trajectories, through the study of the Poincaré section ( ẋ = 0) of the solutions of Eq. ( 11) at times larger than the transient initial period. As in the preceding section, we explore the PP at the constant γ = 2 value. We increase F starting with a value F = 1.5 in the first SD up to F = 6.84 in the heart of the second ID, going through the first ID in between (Fig. 1). The critical values of F separating different behaviors have been calculated using the shooting and continuation methods 20 and are given here with a precision of 0.001. Note that these F critical values depend on the γ value. The calculations have been performed for a wide range of ICs. Whatever the parameters, we found at most two types of asymptotic trajectories presented in Fig. 3a (black and red dots) depending on the ICs. The presented points correspond to (0.1, 0.0) (black dots) and (-0.1, 0.0) (red dots) ICs. λ L , the greatest non-null Lyapounov coefficient (binary base) of the autonomous system associated to Eq. ( 11), has been calculated also to characterize the chaotic or not character of the trajectories. 19 Since the results are in-FIG. 3. a) Poincaré section of the phase space omitting transient regime as a function of F with γ = 2.0 from trajectories calculated by numerical resolution of Eq. ( 11) with initial conditions (0.1, 0.0) (black) and (-0.1, 0.0) (red). Blue curve: U eff (X), the effective potential (abscissa) in arbitrary units versus position X (ordinate) b) lambdaL the greatest non-null Lyapounov exponent (binary base) of Eq. ( 11) as a function of F for initial conditions (0.1, 0.0). The solid horizontal line (lambdaL = 0) is a guide to the eyes.

dependent of the ICs, we present in Fig. 3b its evolution with F for the (0.1, 0.0) IC only.

From Fig. 3 we see that in the first SD i.e. F < F1 , with F1 = 3.028 the first SD/ID limit, the asymptotic solution of Eq. ( 11) is the stable fixed point (0, 0) whatever the ICs (sec. II C 1). This result is coherent with the values of λ L < 0.

At F1 , λ L goes to 0 for the first time and we enter the first ID, the system undergoes a Hopf bifurcation towards an unique limit cycle (with twice the StAW period), the attractor of the flow (identical black and red points) in agreement with the results of sec. II C 2. More precisely, due to the translation invariance of Eq. ( 11), there is actually an infinite number of limit cycles, attractors of the system, one per fixed point (0[π], 0). This domain of periodic asymptotic trajectories extends from F1 to F2 = 3.654. Beyond, for F > F2 , λ L becomes positive and thus the trajectories chaotic. The first periodic-chaotic transition is thus at F = F2 .

In this periodic domain F1 < F < F2 , at F = 3.552, λ L goes to 0 and the limit cycle splits into two limit cycles (separation of the black and red points) without any period change. The system has then two attractors per fixed point. This doubling of the number of attractors, precedes for each attractor a cascade of period doublings at increasingly close values of F , each of them being associated to the cancellation of λ L . We have been able to observe 4 period doublings (at 3.632,3.650,3.653,3.654), though we did not try to optimize this number. Each period doubling corresponds to the apparition of the corresponding peak in the Fourier spectrum. The transition to chaos proceeds when the spectrum becomes continuous at F2 . Between F2 and F3 = 4.890, there is an alternation between large domains of chaotic trajectories and many very small domains of periodic or unbounded asymptotic solutions (with λ L < 0) . The chaotic trajectories are characterized by positive values of λ L and Poincaré sections containing an infinite number of points in the limit of an infinite trajectory. In the very small domains where λ L becomes negative, the solutions are either asymptotically periodic with, in the Poincaré section, a reduced number of points, or unbounded with an absence of points in the x-interval chosen for Fig. 3a.

At F3 = 4.890, λ L goes again to 0 and there is a transition towards unbounded asymptotic trajectories (λ L < 0 and absence of points in the x-interval of the Poincaré section Fig. 3a). This unbounded domain extends up to the end of the first ID at F = F4 = 5.090 where again λ L = 0.

For F4 < F < F5 , where F5 = 6.183 corresponds to the second SD/ID limit, the asymptotic solution of Eq. ( 11) is again the stable fixed point (0, 0) with λ L < 0, whatever the ICs (sec. II C 1).

At F5 we enter the second ID and we observe a behavior very similar to the one in the first ID. The main difference is that the first periodic domain starts directly with two limit cycles, then there is also a cascade of period doublings, the first one at F = 6.522, leading also to a region of alternation of wide chaotic domains and small periodic or unbounded domains.

The system described by Eq. ( 11) is very rich and complex. The apparition of chaotic solutions is actually not surprising: imposing a standing acoustic wave involves the interaction of two harmonics in a parametric-like excitation. The equation (11) without dissipation has been studied in order to investigate regimes where resonances induced by both harmonics interact. [START_REF] Escande | [END_REF] The kicked rotor 21 is also very similar to our system excepted that an infinity of harmonics are involved in the excitation, so that their resonances can interact leading to chaotic solutions. Finally let's mention the work of Van Dooren 20 who studied the dynamics of a pendulum with a forced sinusoidal horizontal support motion: in the absence of gravity, this system reduces to our undamped system (γ = 0).

From this numerical study, it appears that the trajectory of an ad-atom on a substrate submitted to a StAW can be of four different types, depending on the parameters (γ, F ). While increasing values of F at constant γ, the domains corresponding to the different kinds of trajectories are successively i/ converging trajectories in the SD domains, followed in the ID domains by ii/ a periodic domain, then iii/ an alternation of mainly chaotic solutions together with small periodic or unbounded domains and ends up with iv/ an unbounded domain, and so on when entering the next SD. While the structuring effect of the StAW is obvious in the case of converging trajectories and to a less extend periodic trajectories, since it drives the particle into given regions in the configuration space, it is less obvious for the other cases, in particular for the chaotic domains.

III. STOCHASTIC FORCES

In this section, we reintroduce the Gaussian stochastic force ξ(t) in Eq. ( 11) in the τ c → 0 limit:

ẍ( t) + γ ẋ( t) = ξ( t) + F cos( t) sin(x), (25)
with

ξ( t) = 0, (26a)
ξ( t) ξ( t + τ ) = 2 Dδ(τ ), (26b)
where D = D/(ωm 2 c 2 s ) (Eqs. ( 3) and ( 5)). Note that D and γ are not independent due to the fluctuationdissipation theorem (Eq (2)):

D γ = k B T mc 2 s (27)
We will use in this section the same 2.0 fixed value of γ and the same four values of F as in the previous section, corresponding to the four exhibited types of solutions. Three values of D (0.0001, 0.01 and 0.1), or equivalently of temperatures for given substrate (c s ) and ad-atom (m), will be investigated, corresponding to the quasi-absence of the stochastic force (Sect. II, Eq. ( 11)), a medium and a strong stochastic force respectively.

Due to the stochastic character of Eq. ( 25) we adopt here a statistical point of view. It appears, as demonstrated in this section, that on a longtime scale compared to the relaxation time τR the memory of the IC is lost and the system is in a quasi-steady state: P (x), the distribution of the ad-atom x-position is mainly timeindependent (rigorously, it involves a tiny periodic contribution at the StAW frequency). On a timescale larger than the StAW period, the structuring effect of the StAW will be revealed through the correlation between the x variations of P (x) and those of the periodic effective potential U eff . Physically, P (x) will point out the preferential sites where the ad-atom spends most of its time.

Assuming the ergodicity of the system, P (x) can be obtained from a single long trajectory of one particle (after elimination of the initial transient period, with whatever IC) or from a set of trajectories. In addition, due to the translational invariance x → x + π, t → t + π of Eq. ( 11), P (x) is expected to be π-periodic. We calculate P (x) from a number of trajectories obtained for different ICs and realizations of ξ( t): 100 trajectories of t = 10000 time units each, The ICs are taken at random in the x = 0 effective potential valley with ẋ = 0. The points of the trajectory outside the initial valley are translated back to the x = 0 valley coherently with the translational invariance of Eq. (11). P (x) is then estimated from the histogram of the ad-atom position of these trajectories and the following normalisation condition:

π/2 -π/2 P (x)dx = 1. ( 28
)
The results are reported in Fig. 4 over half the StAW wavelength, i.e. a period of U eff together with the effective potential U eff .

1. Converging trajectories F = 1.5, Fig. 4a. At a very low diffusion coefficient value ( D = 0.0001), P (x) is strongly peaked at x = 0: in the quasi-absence of fluctuating force, the trajectory still converges to the minimum of the effective potential. The tiny stochastic force induces small fluctuations of the position in the vicinity of the minimum. This fluctuations are not strong enough to induce a transition to an adjacent valley (inter-valley transition) on the simulation time scale (P (±π/2) = 0 at the maxima of U eff ). At D = 0.01, the width of the peak centered on the fixed point of Eq. (11) has not sufficiently increased to induce a significant inter-valley transition, whereas it does at D = 0.1 (P (±π/2) ≃ 15%): the ad-atom diffuses from a basin of attraction to a neighboring one. The converging character of the trajectory is lost. However the selforganization is preserved through the peaked character of P (x) centered on the minima of U eff .

2.

Periodic trajectories F = 3.15 -Fig. 4b. Here also at very low diffusion coefficient D = 0.0001, the periodic character of the trajectory is roughly preserved with small fluctuations around the initial periodic trajectory in the absence of ξ; and the particle visits a wide region of a basin of attraction of the effective potential. The peaks and features observed on the plot of P (x) for D = 0.0001 are due to the specific shape of the limit cycle (or trajectory shown in Fig. 2b). Increasing the diffusion coefficient D induces some fluctuations around this limit cycle. They can even activate the crossing of the effective potential barriers (on the simulation time scale), clearly evidenced for D = 0.1 by the significant value of P (π/2). Increasing the diffusion coefficient smooths the structural role of the StAW: the stochastic fluctuations give rise to a wide Gaussian-like distribution centered on the effective potential minimum and whose width increases with the diffusion coefficient.

3. Chaotic trajectories F = 4.5 -Fig. 4c. In the quasi-absence of stochastic force ( D = 0.0001) the chaotic character of the trajectory leads to a P (x) distribution correlated to U eff : even if the ad-atom is not trapped in a given potential valley (P (π/2) = 0 and Fig. 2c), P (x) has a pronounced maxima at the minima of U eff . Increasing the diffusion coefficient yields the same qualitative observations as in the previous periodic trajectory case.

4. Unbounded trajectories F = 4.96 -Fig. 4d. At the very low D = 0.0001 value, from the staircase character of the trajectory presented in Fig. 2d one expects a non uniform P (x) with a marked peak inside each potential well. A precise examination of the trajectory in Fig. 2d reveals that the plateaux are at x positions slightly larger than the center of the potential wells. Such a trajectory contributes then to the x1 > 0 peak observed on Fig. 4d (top). The second peak at the symmetric -x 1 position results from mean trajectories towards x < 0 positions. The x1 > 0 non centered position of the distribution associated to trajectories towards increasing mean x values results from the definite direction x > 0 or x < 0 of the observed trajectories (x > 0 in Fig. 2d). Increasing the diffusion coefficient results in an increased width of each peak, leading to a unique central peak for D = 0.1. As in the three preceding cases, the structuring effect of the StAW is also evidenced in that case whatever the strength of the stochastic force in the studied range.

As a conclusion, we would like to emphasize that for all the trajectories types, and as soon as the stochastic force is significant, the distribution P (x) is a Gaussianlike distribution centered on the effective potential minimum. Therefore, the StAW has a structuring effect on the diffusion of the particle: it induces preferential sites in which the particle spends more time. These sites precisely correspond to the minima of the effective potential evidenced in Sect. II. The stochastic force essentially counterbalances the structurating role of the StAW by smoothing the distribution P (x).

IV. RETARDED EFFECTS

We have neglected the correlation time τ c in Sects. II and III. In surface diffusion problems, such an approximation is valid for heavy ad-atoms [START_REF] Forster | Hydrodynamic fluctuations, broken symmetry, and correlation functions[END_REF] which is a very peculiar case. In this section, we investigate the effect of a significant value of the correlation time τ c compared to the StAW period. The correlation time τ c is involved in the friction term through the memory kernel (Eq. ( 4)) and in the auto-correlation function of the stochastic force ξ( t) (Eq. (3b)). Reintroducing the retarded effects Eqs. (25), (26a) and (26b) write now:

ẍ( t) + γ τc t -∞ ẋ( t′ )e -( t-t′ )/τc d t′ = F cos( t) sin(x) + ξ( t), (29) ξ
( t) = 0, (30a) ξ
( t) ξ( t + τ ) = D e -|τ|/τc τc , (30b)
with the same Eq. ( 27) between D and γ. Let's first estimate a physical range for the time τc = ωτ c . On the one hand, in the case of the diffusion of an ad-atom, 1/τ c is roughly of the order of the Debye frequency f D , corresponding to the maximum frequency of atomic vibrations in the crystalline substrate, i.e. f D ≈ 10 13 Hz for common crystals, and consequently τ c ≈ 0.1 ps. On the other hand, as mentioned in the introduction, the StAW wavelengths of interest vary from few to hundreds of nanometers, i.e. typically from 5 nm to 1 µm. With typical sound speeds in solids around 3000 m.s -1 , the StAW frequency ω lies in the range 1.25 10 10 -2.5 10 12 rad.s -1 . Consequently, τc will be in the range 0.002 -0.40. 29), (γ, F ) = (2.0, 1.5), with a small diffusion coefficient (D = 0.0001) and initial conditions (x, ẋ) = (1.0, 0.0), for three different values of τc: 0.001 (black), 0.3 (red), 0.4 (green). Blue curve : U eff (X), the effective potential (abscissa) in arbitrary units versus position X (ordinate).

Three typical values of τc , 0.001, 0.3 and 0.4 will be used to investigate the effect of τ c in the two extreme cases: quasi absence ( D = 0.0001) and significant ( D = 0.01) fluctuations, two values of D already used in the preceding sections. Due to the separation of the different time scales, we expect the solutions of Eq. ( 29) with a low τc value (0.001) to be very similar to solutions of Eq. ( 25). Only values of τc non negligible compared to the StAW period are expected to produce solutions of Eq. (29) significantly different from those of Eq. ( 25). We reduce our study to the converging trajectory case previously studied (γ = 2, F = 1.5), with the same (1.0,0.0) IC as in Sect. II C 1, since this type of trajectory is the most favorable for self-organisation and thus will be usually preferred in any application.

Eq. ( 29) is numerically solved using a Leap-frog algorithm. 22 ξ( t) values satisfying Eqs. (30a) and (30b) are generated with the algorithm of Ref. 23 while the integral of the friction term including the memory kernel are calculated using the algorithm given in Ref. 13.

The solutions for D = 0.0001 are reported on Fig. relaxation time τR of the ad-atom dynamics. Such behavior has already been mentioned and explained in the literature: the velocity autocorrelation function for the ad-atom shows both a ballistic and a diffusive regime, the width, roughly τR of the balistic regime, is related to τc . [START_REF] Gordon | [END_REF] At odds, when the stochastic force is large ( D = 0.01), the possibility for the ad-atom to cross an effective poten-tial barrier on the simulation time scale increases with τc : we report in Fig. 6a the trajectories of the ad-atom as a function of time for the three different mentioned values of τc . Increasing τc , and hence the ballistic regime, gives more importance to the very unlikely extreme values of ξ (t) and thus results in a higher effective diffusion coefficient (not shown). Fig. 6b reports the P (x) distributions: τc slightly affects the histogram P in the vicinity of the minimum of U eff . The effect is more pronounced on its wings: they increase significantly leading to non zero P values at the U eff maxima for τc = 0.3 and 0.4, coherently with the possibility for the ad-atom to escape from its original potential valley evidenced in Fig. 6a). Nevertheless, no matter the value of τc , the shape of P (x) is still gaussian-like evidencing the structuring effect of the StAW

We can thus conclude that if the retarded effects quantitatively modify the trajectories, they weakly affect the structuring effect induced by the StAW.

V. EFFECTIVE CRYSTALLINE POTENTIAL

In this section, we consider the additional effect of the effective crystalline potential Φ eff (x) on the motion of the ad-atom, in the same negligible correlation time limit (τ c → 0). Eq. ( 10) writes then:

ẍ( t) + γ ẋ( t) = -φ0 sin(nx) + F cos( t) sin(x) + ξ( t) (31)
The calculations will be performed with n = 24, a value comparable to the ones we used in our MD simulations. 9 First, we study the modifications induced by Φ eff on the fixed-point stability diagram described in Sect. II -Fig. 1 in absence of stochastic forces. Note that in Eq. (31) both the crystalline and StAW forces cancel for x = 0. We could easily imagine that a dephasing of Φ eff compared to the StAW will shift the fixed points or yield to the absence of fixed points. However, since the StAW wavelength is large compared to the lattice parameter, we do not expect such a dephasing to qualitatively modify the ad-atom trajectories, especially in the presence of the stochastic force. Fig. 7 reports, as in Fig. 1, the stability diagram in the (γ, F ) PP for increasing values of φ0 = 0.0, 0.05, 0.15, 0.25 and 0.35. The main effect of Φ eff is to shift the instability domains towards higher F values at constant γ. The effective crystalline potential has then a stabilizing effect on the trajectories.

To investigate the structuring effect of the StAW, as in the preceding section, we will study the effect of Φ eff in the presence of the stochastic force only for one couple (γ = 2, F = 1.5), corresponding to converging trajectories in the absence of crystalline potential. Here we are mainly interested on what happens when Φ eff is substantially higher than the effective StAW potential U eff and is thus potentially able to challenge or to overcome the structuring effect of the StAW: actually, the StAW 8, the structuring effect seems to increase with the increasing of the effective crystalline potential. However, this effect is essentially due to the normalization: fitting the set of maxima of P (x) by a Gaussian curve leads to approximately the same Gaussian width for all the values of φ0 .

Hence, the structuring effect is weakly affect by the values (even for significant values) of φ0 . Note the very good qualitative agreement between the distribution P (x) in Fig. 8 and the results reported in our Molecular dynamics simulations 9 .

VI. EFFECTIVE POTENTIAL AND PROBABILITY DISTRIBUTION

To quantify the structuring effect of the StAW, we note that the structuration can be characterized by an effective energy difference ∆E eff SAW at the mesoscopic scale between the minima and maxima of the effective potential deduced from the P (x) curves of Fig. 8:

∆E eff SAW = k B T ln [max(P max )/min(P max )] (32)
with P max the ensemble of local maxima of P (x) and max (min) the ensemble maximum(minimum). Since both effective energy differences, ∆E eff SAW (Eq. ( 32)) and effective potential U eff (Eq. ( 22)), govern the diffusion of the ad-atom, and since U eff quadratically depends on the amplitude of the force F , we can reasonably expect the same quadratic dependence for ∆E eff SAW . However, note that while U eff does not take into account the stochastic force ξ(t) nor the effective crystalline potential, ∆E eff SAW implicitly takes into account all these contributions. Fig. 9 exhibits the dependence of ∆E eff SAW as a function of F 2 for γ = 2, φ0 = 3.0 and for D = 0.05, 0.1 and 0.2. ∆E eff SAW has been calculated from P (x) curves similar to the ones of Fig. 8. We did not succeed to calculate ∆E eff SAW as a function of F for D < 0.05: the average time needed by the particle to escape from an effective crystalline potential being too long to obtain good statistics, due to the activated character of this event with a 1/D exponential dependence. The linear dependence of ∆E eff SAW as a function of F 2 is observed in the small forces regions where the perturbative calculations of Sect. II A scientifically sound. We note that increasing the diffusion coefficient or the friction γ (not shown) yields a smaller StAW structuring effect i.e. as already mentioned, the thermal noise challenges the StAW structuring effect.

The results of this section demonstrates that the StAW amplitude and the temperature, acting on the friction γ and the diffusion coefficient D allow to tune the structuration induced by the StAW. These information are essential to identify the key parameters if one wishes to experimentally implement the dynamic substrate structuring effect described here.

VII. CONCLUSION

In this work, we have studied the solutions of the motion equation of an ad-atom diffusing on a substrate submitted to a StAW. To this aim, we have applied some standard tools of the non-linear physics to a simplified version of the motion equation keeping the most relevant terms. Noticeably, an effective potential governing the slow dynamics of the ad-atom has been derived. We have shown that this effective potential controls the distribution of the ad-atom x-position even when considering the other additional terms (noise, retarded effects, effective crystalline potential) in the ad-atom motion equation. We underline the relevance of our model when comparing calculated distribution probabilities to positions histograms collected from molecular dynamics simulations 9 .

Our present study opens different perspectives. In a recent publication 9 , we have discussed about the possibilities to experimentally implement the dynamic substrate structuring effect. We have noticeably mentioned the impossibility to produce StAW with a wavelength in the nanometer range on the substrate surface using the current available experimental setups: to our knowledge, and in the state of the art, the production of surface acoustic waves with wavelengths in the 100 nm range is possible using optical excitations 24 . However, the dynamic structurating effect does not directly depend on the wavelength. Actually, the dynamic effect is expected to exist as long as the ad-atom does not have the time to diffuse on a wavelength λ during a period of the StAW i.e.

k ≪ c s D eff ( 33
)
where D eff is the effective diffusion coefficient of the adatom in the effective crystalline potential. One may thus consider the possibility to use StAW with hundred nanometers wavelengths. Our present model and study will then be a very fruitful tool to evince the optimized parameters (especially the temperature and the StAW amplitude) leading to an efficient structuring effect. Theoretically, Eq. ( 11) exhibits a very rich and complex dynamics. We have studied the solutions of this equation, but we currently consider the possibility to derive the Fokker-Planck equations for the ad-atom probability density functions. Such equations would allow to directly derive the position distribution P (x), and even perhaps to derive the effective diffusion coefficient describing the diffusion of the ad-atom between the different basins of the effective potential.

Finally, in this paper, we have focus on the structuring effect induced by the StAW on the position distribution of the ad-atom and have eluded the study of the dynamics of the ad-atom. We have nevertheless exhibited the possibility for the ad-atom to follow different qualitative types of trajectories depending on the parameters γ, F and φ0 . The study of the ad-atom dynamics, including the dependence of the StAW force on the substrate lattice parameter scale is a natural perspective to this work.

  FIG.2. Black and red curves: x( t) solutions of Eq. (11) for γ = 2.0. Stability domain (SD): a) converging trajectories ( F = 1.5). Instability domain (ID): b) periodic trajectories ( F = 3.15), c) chaotic trajectories ( F = 4.5) and d) unbounded trajectories ( F = 4.96). Initial conditions are reported in the legend. Blue curves: U eff (X), the effective potential (abscissa) in arbitrary units versus position X (ordinate).

FIG. 4 .

 4 FIG. 4. Black, red and green curves : Histograms of P (x) for different values of the diffusion coefficient D: 0.0001, 0.01 and 0.1 respectively, with γ = 2.0 and different F values, i.e. types of trajectories as defined in sec. II C: (a) converging ( F = 1.5), (b) periodic ( F = 3.15), (c) chaotic ( F = 4.5) and (d) unbounded ( F = 4.96). Blue curve: U eff (X), the effective potential in arbitrary units.

4 FIG. 5 .

 45 FIG.5. Converging trajectories x( t) solutions of Eq. (29), (γ, F ) = (2.0, 1.5), with a small diffusion coefficient (D = 0.0001) and initial conditions (x, ẋ) = (1.0, 0.0), for three different values of τc: 0.001 (black), 0.3 (red), 0.4 (green). Blue curve : U eff (X), the effective potential (abscissa) in arbitrary units versus position X (ordinate).

5 .FIG. 6 .

 56 FIG.6. (a) : converging trajectories x( t) solutions of Eq. (29), (γ, F ) = (2.0, 1.5), with a diffusion coefficient D = 0.01 and initial conditions (x, ẋ) = (1.0, 0.0), for three different values of τc: 0.001 (black), 0.3 (red), 0.4 (green). Blue curve : U eff (X), the effective potential (abscissa) in arbitrary units versus position X (ordinate). (b): Corresponding histograms of P (x) from 100 trajectories of 10 000 time units with initial conditions evenly distributed between -π and π (same color code as in (a)).

FIG. 7 .

 7 FIG. 7. Stability (blank) and instability (dashed) domains of the fixed point (0,0) of Eq. (31) in the (γ, F ) parameter plane (PP) for increasing values of the effective cristalline force φ0. From top to bottom : φ0 = 0.00 (same as Fig. 1), 0.05, 0.15, 0.25 0.35.

FIG. 8 .

 8 FIG. 8. Histogram of P (x) calculated from converging trajectories ((γ, F ) = (2.0, 1.5)) solutions of Eq. (31) (100 trajectories of 200000 time units each), with D = 0.1 and increasing values of φ0 from top to bottom: φ0 = 0.7, 1.5, 3, 5.

2 FIG. 9 .

 29 FIG. 9. Effective energy difference ∆E effSAW as a function of ratio φ0/ F calculated from histograms P (x). Histogram P (x) are calculated from 100 trajectories of 200000 time units for each values of φ0/ F by solving Eq. (31) with γ = 2.0, φ0 = 3.0, D = 0.001, 0.01, 0.1.


Luc Touraille
email: touraille@isima.fr

Mamadou K Traoré
email: traore@isima.fr

David R C Hill

On the Interoperability of DEVS Components On-Line vs. Off-Line strategies

Keywords: DEVS, interoperability, web service, XML Résumé DEVS, interopérabilité, service web, XML

During the last years, the DEVS community provides many contributions towards the realization of a worldwide platform for collaborative Modeling & Simulation. The goal of such a platform would be to enable the sharing and reuse of models between scientists, as well as the seamless simulation of distributed and heterogeneous models. Therefore, one of the major research fields is the definition of architectures for integrating heterogeneous DEVS components, meaning simulators and/or models written in different frameworks and programming languages. In this work, we present three different strategies for providing such interoperability between DEVS components. The first focuses on standardizing exchanges between simulators, and has been explored in previous works. The two others strategies are more prospective; in keeping with Model-Driven Engineering, they place the model at the center of their architecture and make extensive use of model transformations. To make this possible, we defined a platform and language-independent format for describing and sharing DEVS models, called DEVS Markup Language.

Introduction

One of the current trends in software engineering is to develop applications not as monolithic entities like they were a few years back, but as a collection of loosely-coupled/highly-cohesive components that can be easily combined, substituted and reused, named services. In addition to the notion of "Software As A Service" introduced a decade ago [Keith at al. 2000] and now knowing business success, the Service Oriented Architecture (SOA) is now commonly implemented through web services, a widely used technology based on several accepted standards [START_REF] Newcomer | Understanding Soa With Web Services[END_REF]. Among the advantages of web services, we retain the fact that they facilitate the integration of heterogeneous components, notably through the standardization of interfaces description (Web Service Description Language), messages definition (XML Schema) and communication protocol (Service Object Access Protocol). Therefore, we can leverage this technology to provide straight-forward interoperability of DEVS components [START_REF] Zeigler | [END_REF]].

Our goal is to provide a simple and mostly automated way of executing simulations that involve distant and/or heterogeneous4 DEVS models. This can be achieved by taking two different approaches: simulator-based interoperability and model-based interoperability. Both these approaches need a standardized, platformindependent representation of DEVS models. We defined the DEVS Markup Language, an XML-based language for describing and sharing DEVS models.

Simulator-based interoperability

The main idea of this approach, as used in DEVS/SOA [Seo 2009], is to have a collection of simulation services distributed over the internet. These services provides several operations for simulating atomic or coupled DEVS models in a unified manner, by using the DEVS simulation protocol and the closure under coupling property of coupled models to make them appear as atomic model. The overall simulation is coordinated by a main service, which acts like an entry point for the user. This architecture is summarized in the Figure 1. In order to communicate, the DEVS services expose a standard interface that accepts the usual simulation protocol messages: initialize, get time of next event, run transition and so on.

From a user perspective, the simulation process consists in:

1. writing a coupled DEVS model, using its favorite framework. 2. providing a list of DEVS servers on which he wishes to distribute the simulation. 3. deploying DEVS models to the servers (either by downloading them from his location or activating existing remote models) 4. running the simulation

Model-based interoperability

This second type of interoperability is based on the fact that most of the time modelers already have or can easily install a DEVS simulator on their machine. However, it is much more difficult for them to retrieve and reuse models already existing on other scientists' computers. Therefore, we need some way to provide access to and share these numerous models.

In this perspective, we distinguish between on-line (or dynamic) and off-line (or static) interoperability, which are both centered on the models, and bring into play several methods related to Model Driven Engineering (MDE) [Schmidt 2006].

On-line (dynamic) model-based interoperability

In this solution, models themselves are deployed as web services instead of simulators (see Erreur ! Source du renvoi introuvable.). Using this model driven approach, the operations invoked through the network are no longer simulation mechanisms, but model functions, such as δ int , δ ext , the time advance function, etc. Given a model written for a specific framework, we need to generate its representation in a platform-independent language, such as the DEVS Markup Language (DML) [START_REF] Touraille | A Mark-up Language for the Storage, Retrieval, Sharing and Interoperability of DEVS Models[END_REF], described below, or the DEVS Modeling Language (DEVSML) [START_REF] Mittal | [END_REF]]. This description will then be used to generate:

-the web service description, as a WSDL file -the wrapper that will adapt the model interface to the service interface, by forwarding the operation calls.

Once this is done, the client can retrieve the description of the model, and generate a stub that will conforms to its framework and behave exactly as the original model, by invoking the wrapper service. The automatic generations cited above can be seen as automatic transformations in the model driven engineering terminology [Schmidt 2006]. Regarding distant coupled models, two solutions are possible. The first is to use some kind of adapter, as in the previous section, which will make the coupled model appear as an atomic one. The second possibility is to take into account that coupled model don't have any operational semantics. All they do is describing their ports and their couplings, in a static way. As a consequence, we can use their DML representation to generate a copy of the model in the local framework, without losing any useful information.

In a nutshell, the model provider:

1. writes a model in his favorite DEVS framework.

2. automatically generates the model's DML description, from which he generates a WSDL file and a web service encapsulating the model. 3. deploys the web service, providing access to the underlying model.

On his side, all the model consumer has to do when writing a coupled model in a given framework is:

1. find the location of the models he needs, using some kind of model directory or more classical discovery procedures such as web searches or acquaintances. 2. give the addresses of the models' web services to a generation tool that will download the WSDL files and generate client stubs adhering to the modeler framework. 3. run the simulation. The framework's simulator deals with local and distant models without even knowing it.

Off-line (static) model-based interoperability

On the previous approach, we used DML description of models to generate client and server stubs; however, DML includes not only the interface of models, but also their dynamics, described in a platform and language independent format (see Section 3). Consequently, we can go one step further in automatic model transformations and generate entire local copies of distant models.

Using transformations based on code parsing, we can generate the DML description of existing models written for a specific framework. These DML files can be stored in model repositories, accessible through Internet, so that anyone (with the proper authorizations) can access them. When facing a problem, the modeler starts by searching an existing model that could fulfill his requirements. If there is none, he can write his own but still take advantage of existing works by reusing them in his coupled models. To do so, the user download the model description from its repository then use generic transformations to generate the code corresponding to the model in his framework. He ends up simulating his models and those of other scientists on his own simulating platform.

Compared to on-line interoperability, this approach has the benefit that it drastically reduces the communications with remote servers. In fact, after the model description has been downloaded, there is no more need for accessing the network (hence the "off-line" appellation). However, the drawback of this is that a given model can have several unsynchronized versions all over the world. A solution to this issue could be to use some sort of update feed, but this is still more complicated than simply invoking a service without caring about the underlying implementation, which can continuously improve and evolve.

To sum up both interoperability solutions we exposed, we can say that the simulator-based approach is more aimed at coupling distributed or local models over distant servers, whereas the model-based approach aimed at integrating distributed models in a local simulation.

DML, a markup language for the sharing and interoperability of DEVS models

As stated in the previous sections, many solutions for providing seamless interoperability of DEVS models involve a standardized language for describing these models. The DEVS Markup Language (DML) [START_REF] Touraille | A Mark-up Language for the Storage, Retrieval, Sharing and Interoperability of DEVS Models[END_REF]] is an attempt to unify previous propositions for such a language [Janoušek et al. 2006] [Mittal et al. 2007] [Risco-Martín et al. 2007], while at the same time making it as generic and flexible as possible.

Introduction

An urgent necessity in the Modeling and Simulation domain is to define standards to make for interoperability and reusability of models and simulators. Thanks to these standards, a modeler should be able to seamlessly execute the following operations:

-store model specifications and data trajectories for further use.

-perform reasoning on these data, e.g. select a model in a library based on precise criteria, or execute analysis on a model structure or a generated output trajectory. -share models among the community, to enhance productivity and quality by reusing well-tested and validated models. -visualize simulation results through different interfaces (charts, animation), and graphically create/edit models. -transform model specifications into executable source code for several target simulators, and perform simulation.

Due to the lack of standard, these use cases are currently rather difficult to fulfill. Most actual models are described either in natural language, problematicif not impossibleto process with a computer, or in a specific programming language, using types and functions peculiar to one simulator. This heterogeneity obstructs the collaboration between modelers, and makes impossible the development of generic tools that could be used by the entire community.

To solve this problem, we propose a simulator-independent language to describe models formalized in the Discrete Event System Specification (DEVS) [START_REF] Zeigler | [END_REF]]. This language, built on eXtensible Markup Language (XML) [Abiteboul et al. 1999], must allow the specification of all the aspects needed to realize the previously exposed use cases.

The next section presents the works related to the definition of an XML representation of DEVS models. Section 3 exposes several limits of these previous propositions, and propounds some ideas that are further developed in section 4. Finally, we synthesize our results and present future work in section 5.

Related works

The DEVS Standardization group pinned down four areas of the DEVS framework needing standardization [START_REF] Vangheluwe | [END_REF]]:

-the DEVS formalism in itself, and its variants -model representation, which should describe the model structure and dynamics in a platformindependent manner -minimum requirements for a simulator to be labeled "DEVS-compliant" -model libraries, aimed at providing a collection of models usable out of the box This paper focuses on the second point, namely the standardization of model representation. To achieve this goal, the choice of XML seems obvious; indeed, a consensus emerged among the M&S community about XML being the best tool to describe models in a standardized way [START_REF] Brutzman | [END_REF][START_REF] Fishwick | [END_REF]]. This general agreement is due to numerous reasons, such as the wide use of XML in Service Oriented Architecture [START_REF] Newcomer | Understanding Soa With Web Services[END_REF], the great number of existing XML tools, its platform independence and many others.

Several works proposed XML Schemas or Document Type Definitions (DTDs) describing languages to represent DEVS models. [START_REF] Wang | [END_REF] develop a tool allowing the modeler to graphically build System Entity Structures [START_REF] Zeigler | [END_REF]], which can be used to capture the DEVS models structure. The graphical representation is then transformed into an XML document containing the model hierarchy.

Based on this premise, [START_REF] Rhöl | [END_REF] define XML Schemas for the DEVS formalism, including both coupled and atomic models. Model structure specification includes ports, couplings and state variables, and the emphasis is put on the use of XML data binding to automatically generate classes corresponding to XML models. However, the representation of the model dynamics is rather limited since each function is defined by a raw string, written in a specific programming language. This problem of behavior representation is tackled in several works. One approach [START_REF] Risco-Martín | [END_REF]] is to restrict the different functions to simple constructs. For example, the transition functions are described by a set of conditional expressions, each being associated with new values for the state variables. A more general and powerful solution is to use some kind of pseudo-language, as is proposed in [Janoušek et al. 2006]. The advantage of using a simple language is that it can be easily represented in XML and transformed into source code in any programming language. Finally, [START_REF] Mittal | [END_REF]] use existing XML-based source code representation. In this case, the favored language is Java, but the authors evoke the possibility to use a more generic XML language to describe the model logic.

Limits and perspectives

The major problem that emerges from these previous works is the difficulty to standardize the representation of the model dynamics. Several solutions are possible, each with their own advantages and drawbacks.

The first idea is to use an XML representation of a specific programming language, e.g. JAVA Markup Language [Badros 2000]. This approach has the great advantage that many tools are available to manipulate and transform the source code and the XML representation. However, it induces a strong coupling between the model and a particular language, therefore restricting its deployment scope to a small set of simulators.

To enlarge this scope, a more generic representation could be employed, such as Metal [START_REF] Lemos | [END_REF]] or Object Oriented XML (O2XML) [START_REF] Wiharto | An Architecture for Retargeting Application Logic to Multiple Component Types in Multiple Languages[END_REF]. These markup languages enable a representation of most object-oriented programming languages, by specifying a common denominator embracing all base features of these languages, in a .Net Common Language Specifications fashion. Thereby, an XML document can be transformed into source code in any object-oriented target language. Even though this solution might seems very attractive, it prevents the modeler from using certain language features or libraries that could be crucial. For example, one might want to generate pseudo-random numbers with the Stochastic Simulation in Java library [L' Ecuyer et al. 2002], or perform matrices operations using Boost.uBLAS. As a result, such a high-level language with a small set of features would reduce the universality of the standard, and make very hard the representation of legacy models, which are usually written in a specific programming language and take advantage of its functionalities and libraries.

To resolve this issue and be as permissive as possible, we could just let the modeler write the model functions in his favorite language, and delegate the task of interpreting them to simulators or proxies. However, that would imply developing an architecture where a model could be composed of several modules in multiple languages, this model being itself simulated by an application possibly written in yet another language. Even though this architecture is conceivable, using technologies such as JNI, RMI, CORBA or Web services, the complexity and the loss of performance would be huge.

To reach a compromise between genericity and flexibility, we define in this work our own XML representation of source code, by drawing inspiration from all these approaches. In addition, we propose some directions about several aspects that were left unspecified in previous works, such as parameterization and initialization of models, representation of their graphical attributes and standardization of trajectories.

Towards a unified representation

In this section, we propose a basis for a unified and standardized representation for each characteristic of DEVS models, as described above.

Model structure

Describing the structure of DEVS models is rather straight-forward, and is simply a translation from the set formalism into an XML representation. Figure 3 depicts our XML Schema arborescence for atomic models in classic DEVS.

Input and output ports, parameters and state variables are characterized by a name and a type. To allow the deployment of the model on any simulator, this type should be as generic as possible; however, in real-world applications, the modeler might need to use certain types defined in particular language libraries. To satisfy these apparently contradictory requirements, we define three different kinds of types:

intrisic types, which are natively supported by all object-oriented languages (e.g. integer, string, float, etc.) custom types, i.e. types that are defined by the modeler in a language independent manner language dependent types, which have to be bind to a particular type, for each targeted languages. <state> <variable> <name>phase</name> <type kind="custom">procPhase</type> </variable> <variable> <name>buffer</name> <type kind="langDepend">Queue</type> </variable> <variable> <name>sigma</name> <type>double</type> </variable> </state> ... <type id="procPhase" xsi:type="enum"> <value>busy</value> <value>idle</value> </type> <type id="Queue" xsi:type="langDepend"> <impl lang="java"> java.util.LinkedList<String> </impl> <impl lang="c++"> std::queue<std::string> </impl> </state> 4 shows an example use of these different types. In this sample, the state variable "sigma" is a simple floating-point number, whereas "phase" has a more complex type: it is an enumeration, described below as having two possible values, busy and idle. Finally, the variable "queue" has an abstract type, which we bind to specific types in Java and C++. Therefore, to make this model deployable on a simulator written in another language, the modeler simply has to specify what type can be used as a queue in this language. We can also imagine having repositories of abstract types along with their bindings in numerous programming languages, to allow their reusability in several models.

The XML Schema for coupled DEVS models has no major additional subtleties. It includes ports, components and coupling declarations, as well as a tie-breaking function in the case of classic DEVS. Components are described by their name and the type of model they instantiate. The couplings are pairs of ports, if necessary qualified by their component's name. Finally, the Select function is described by ordered sets indicating priorities among particular group of components.

Model dynamics

Section 3.3 identifies several issues concerning the specification of atomic DEVS models behavior. These problems boil down to the following question: How can we describe application logic in a language-independent manner? Several works have been done to provide an XML representation of source code. Some of them are oriented towards a sole programming language, e.g. JAVAML [Badros 2000], PascalML [McArthur et al. 2002] or CppML [START_REF] Mamas | [END_REF], while others propose vocabularies that contain only the common denominator between all object-oriented languages [START_REF] Mamas | [END_REF] None of these works provide a framework being at the same time as generic as possible, and open to extension with language dependent features. Therefore, we propose a compromise between both approaches, by drawing inspiration from the concept of weaving used in Aspect Oriented Programming [START_REF] Kiczales | [END_REF], or more simply from the C/C++ macros.

The main idea is to use a pseudo-language such as O2XML so that the major part of the code can be reused across different simulators, while letting the possibility to include language specific snippets. This concept is illustrated in Figure 5, where generic programming constructs are represented in XML, and library uses are abstracted and bind to specific implementations in several languages. This semi-generic language can also be used to create custom types (i.e. classes, enumerations), which become usable in the variables and ports declarations. This approach has two major advantages:

allowing the use of programming language features and libraries.

permitting the representation of legacy models written for a specific simulator.

By limiting the language dependent code snippets to a minimum, we facilitate the migration of models between platforms and thereby enhance their portability and reusability.

Parameterization and initialization

The description of DEVS models as described above is not sufficient to produce an executable specification. Indeed, models can have parameters, which need to be set prior to simulation, and their initial state must be specified in some way.

Figure 6 depicts how a model, regarded as a template, can be instantiated into several parameterized model depending on the values assigned to each parameters. After that, the state must be initialized to obtain a ready-tosimulate model. A simulator can then be fed with this executable model, along with an input trajectory if need be, and produce the output of the simulation. To perform these operations, two approaches are possible depending on the degree of freedom required. The simplest way is to restrict parameterization and initialization to pairs of variable identifiers and values. However, these operations can sometimes be more complex than mere assignments, e.g. when dealing with file manipulation or random number generation. Therefore, our proposition let these two activities be described in a similar fashion as the model dynamics, using the not-so-generic source code representation outlined in the previous section.

Consequently, the modeler needs to define two functions for each atomic model in the hierarchy: one to set the parameters values, and the other one to initialize the state. Ideally, these functions should be kept separate from the model definition, so that it can be easily reused with different scenarios.

In order to simulate a DEVS model, four specifications are needed: structure, dynamics, parameterization and initialization. Given these specifications and if required an input trajectory, a simulator can produce results representing the model behavior.

Model behavior

The behavior of a model is normally fully represented by its output trajectory. However, it is common to desire its state trajectory and, in the case of coupled models, the behavior of its components as well. These two aspects can be described by set of values annotated with a timestamp. Figure 7 shows a sample state trajectory of a coupled model, namely the pipeSimple defined in [START_REF] Zeigler | [END_REF]].

The XML schema for this representation of trajectory can be deduced from the model structure: each variable or port become an XML tag accepting a value in accordance with its type. Another solution would be to have a tag variable with an attribute containing the name of the variable.

In the case of variables of complex type, the value at a given time can be obtained by serializing the object from memory, as done when exchanging data between web services. This supposes providing an XML Schema describing the structure of the object, for later use (e.g. in a visualization tool). For custom types, this schema can be deduced from the type definition made in the model specification; for language dependent types, the schema can often be automatically generated from the source code.

<behavior id="pipeSimple"> <state> <snapshot time="0"> <component id="p0"> <phase>idle</phase> <queue /> <sigma>infinity</sigma> </component> <!--other components --> </snapshot> <snapshot time="3"> <component id="p0"> <phase>busy</phase> <queue> <elem>job1</elem> </queue> <sigma>10</sigma> </component> <!--and so on -->

Figure 7: Sample state trajectory

This XML representation of trajectories has the great advantage that it can be easily processed by a machine while remaining human readable; however, due to the potentially huge size of data generated, it might be necessary to switch to a more compact representation, such the eXternal Data Representation standard [START_REF] Vangheluwe | [END_REF].

Presentation

Finally, a modeler should be able to save and load its models to a repository without losing their graphical attributes (position on the layout, size and so on). However, these pieces of information are very dependent on the graphical user interface employed, and therefore should be uncoupled from the model specification, and stored in separate files, in the style of Cascading Style Sheets used in web development.

Conclusion

In this paper we proposed three possible approaches for providing interoperability between DEVS components. The first one, simulator-based interoperability is to standardize exchanges between simulators and have been successfully experimented in DEVS/SOA. The two others focus on models instead of simulators. In the on-line approach, we deploy models as web services so that they can be used remotely by any simulator, using model transformations for generating these services and the code for invoking them. In the off-line approach, we go one step further in the code generation process by generating a complete operational copy of the models on the local simulator. To make possible these two model-driven strategies, we use a platform and language-independent description of DEVS models done with the DEVS Markup Language.

The propositions made here are rather prospective; however, some works have been done towards the realization of the on-line model-based architecture, which promise interesting follow-ups. The off-line solution, involving multiple complex model transformations, need to be looked into more thoroughly to demonstrate its feasibility, and is a work in progress.

Figure 1 :

 1 Figure 1: DEVS components interoperability through simulators communication

Figure 2 :

 2 Figure 2: Local simulation of distant and local models

Figure 4 :

 4 Figure 4: Sample XML representation of state variables

Figure 3 :

 3 Figure 3: XML arborescence for classic atomic DEVS model

Figure 5 :

 5 Figure 5: Representation of model dynamics using a semi-generic language

Figure 6 :

 6 Figure 6: Parameterization and initialization of a model

<codeSnippet id="enqueue" queue="buffer" element="job" />

  [START_REF] Wiharto | An Architecture for Retargeting Application Logic to Multiple Component Types in Multiple Languages[END_REF].

	<if>
	<test>
	<binaryExpr op="equality">
	<lhs>
	<var>phase</var>
	</lhs>
	<rhs>
	<enumValue>
	<enum>procPhase</enum>
	<value>idle</value>
	</enumValue>
	</rhs>
	</binaryExpr>
	</test>
	<then>
	</then>
	</if>
	...

<codeSnippet id="enqueue"> <parameter name="queue" /> <parameter name="element" /> <impl lang="java"> queue

	.offer(element);
	</impl>
	<impl lang="c++">
	queue.push(element);
	</impl>
	</codeSnippet>
	Java
	if (phase == procPhase.idle) {	C++
	queue.offer(element);
	}
	if (phase == procPhase::idle)
	{
	queue.push(element);
	}

By heterogeneous we mean written in different programming languages and/or for different DEVS frameworks.