D Levesque

J. -J Weis

Stability of solid phases in the dipolar hard sphere system

Keywords: solid phases, dipolar hard spheres, Monte Carlo simulation, nonequilibrium work method solid phases, dipolar hard spheres, Monte-Carlo simulations, nonequilibrium work method

Free energy differences between solid phases of dipolar hard spheres are estimated by Monte Carlo simulation using a nonequilibrium work method. These calculations allow to determine which of the considered phases has the minimum free energy. The phase diagram which we obtain is confirmed by simulations in the isothermal-isobaric ensemble over a wide region of the density and temperature domain.

Introduction

We are very pleased to contribute to this special issue in honour of Professor Luciano Reatto and, on this occasion, to remind us of our fruitful collaboration.

The dipolar hard sphere (DHS) system presents a variety of fluid and solid states which have been actively studied in the literature [1]. The low density states, at small dipole moments are similar to those in simple monatomic fluids. At moderate densities and large dipole moments, in the fluid state, the DHS associate into linear chains of variable lengths [2][3][4]. Much endeavour has been devoted as to whether there could be liquid-vapour type coexistence between the low density fluid phase and the one of linear aggregates without, though, being able to provide an unambiguous conclusion about the existence of such a coexistence and the location of its critical point [5][6][7]. At densities typical of the liquid state and large dipole moments simulation results established the existence of a polarized phase for finite size systems with periodic boundary conditions [8][9][10]. The existence of such a phase in the thermodynamic limit has, however, been questioned on the basis of theoretical arguments [11,12]. Nonetheless, the dense fluid phases, polarized or not, coexist with solid phases whose symmetry has, to a large extent, still to be established. The aim of the present work is to contribute to the determination of the diagram of the solid phases and to estimate their stability as a function of density and dipole moment, i.e. temperature in the present case.

Available studies at zero or finite temperatures demonstrated the existence of a variety of solid phases for dipolar particles of spherical shape [13][14][15][16][17]. At low dipole moments there will be no polarization and the stable lattices should, a priori, be similar to those found for spherically symmetric particles interacting by a van der Waals type pair potential. These lattices are of a body-centred cubic (bcc) type at high temperature and face-centred cubic (fcc) or hexagonal close-packed (hcp) at low temperature [18]. For sufficiently high dipole moments such that the solid phase is polarized, the unit cell of the lattice can be modified. If, for instance, there is a polarization along the z-axis the cell edge of the bcc lattice can decrease in this direction and the bcc lattice transforms into a body-centred tetragonal (bct) lattice. If, in addition, the lattice parameters change along the x and y axes but the unit cell remains tetragonal the bcc lattice turns into a base-centred orthorhombic (bco) lattice. Similarly, fcc or hcp lattices can modify under the influence of polarization.

The number of possible lattices in the solid phase, taking into account the dipolar interaction and excluded volume effects, can obviously be quite large. In the present study only tetragonal lattice cells, bcc, fcc, bct, bco and also hcp will be considered. The relative stabilities of the corresponding phases will be established by Monte Carlo (MC) simulation.

A frequently used method to estimate free energy differences between phases is thermodynamic integration [START_REF] Frenkel | Understanding Molecular Simulation[END_REF]. This method assumes that the free energy difference can be calculated along a sequence of equilibrium states connecting a reference state of known free energy to the equilibrium state of the solid phase under consideration. This procedure requires that no phase transition is crossed when going from the initial reference state to the final equilibrium state. The possibility of the existence of first order phase transitions between solid phases as well as the absence of a convenient reference system for systems with rotational degrees of freedom determined us to discard this method in the present study.

Monte Carlo simulations in the isobaric ensemble [START_REF] Frenkel | Understanding Molecular Simulation[END_REF] are in principle suitable for determining solid phase diagrams. With an adequate choice of particle number, MC sampling of the size and the shape of the simulation volume possibly can drive the system to evolve e.g., from a bcc to a fcc phase or conversely, depending on their relative stability at the density and temperature considered. The principal limitation of this method is that despite a favourable free energy difference, MC sampling over an accessible amount of computation time may succeed the system to evolve only from an unstable phase to a mechanically more stable phase but not to the phase of minimum free energy. Due to this limitation we supplemented the simulations in the isobaric ensemble with an estimate of the free energy difference of the solid phases based on the Jarzynsky equality [START_REF] Jarzynski | [END_REF], also referred to as the nonequilibrium work (NEW) method in the literature. The latter allows to calculate the difference between two equilibrium states from the work which is spent when going from one state to the other. This calculation, contrary to those based on thermodynamic integration or determination of the mean work, does not require reversible transformations. However, for given temperature (or dipole moment) and density, the Jarzynski equality does not allow to access the minimum free energy of the solid phase. Indeed, the free energy difference is evaluated between a given initial equilibrium state and final states which are specified transformations of this initial state. The calculated free energies allow only to estimate the relative stability of the chosen final states. When these final states are similar to the solid phases spontaneously obtained in the isobaric ensemble simulations, it is possible to distinguish which among them is the state of lowest free energy.

The article is organized as follows. Section 2 gives details of the simulation methods and the use of the Jarzynski equality. Section 3 presents the simulation results in the isobaric ensemble and for the free energy differences. In Section 4 a tentative solid phase diagram of the DHS system is given in the density range ρ * = ρσ 3 = 1.06 -1.26 and dipole moments µ * = µ/ k B T σ 3 = 1, 2 and 3 where σ and µ are the diameter and dipole moment of the spheres and T the temperature. A summary of the main results is given in the conclusion.

2. Simulation method and free energy calculation.

The simulations in the isobaric ensemble are performed at fixed temperature and pressure, with the standard sampling algorithm alternating trial displacements of the particles and trial variations of the size and shape of the volume containing the system of N hard spheres with periodic boundary conditions. Shape variations are such that the simulation cell remains tetragonal. Having fixed N , the pressure and dipole moment and chosen an initial lattice and volume, the MC sampling will drive this state to one with lower free energy. From the analysis of the pair distribution functions g(r) and the average dimensions of the simulation volume in the final equilibrium state one can determine the symmetry of the lattice. Isobaric-isothermal (NpT) MC were performed for dipole moments µ * = 1, 2 and 3 by applying a uniform pressure p to an initial system having bcc, bct, fcc or hcp structure. By varying p/k B T from 2 to 35 a density range ρ * = 1.06-1.27 is covered. The number of particles was 1024, 1280, 864 and 1440 for the bcc, bct, fcc and hcp lattices, respectively. For each state the initial orientations of the dipole moments were chosen to be random and the total energy of the system was calculated by means of an Ewald sum [START_REF] Allen | Computer Simulation of Liquids[END_REF]. An attempted volume change is achieved by varying the three box dimensions separately according to L ′ α = L α 1+∆ α (0.5-λ α ) where the λ α (α = x, y, z) are random numbers in the interval [0, 1], maintaining the simulation cell orthogonal. The values of ∆ α were adjusted to have an acceptance ratio for the volume changes of about 0.3-0.4.

Ideally one would expect that whatever the initial lattice structure of the system, the MC procedure would drive it to the configuration of lowest free energy (within the constraint of an orthogonal cell). It turned out that for most states considered the initial compact lattices, hcp and fcc, remained stable and did not transform to the lattice structure identified by the NEW method described below as being of lowest free energy. As will be stated below the transformation of an initial hcp structure to e. g. a bcc or bco structure requires relative displacements of alternate basal planes but such moves had not been implemented in the MC sampling. Starting from an initial fcc lattice, which remains quite stable and becomes polarized for µ * > 1, the dipole arrangement most oftenly obtained consists of domains of lattice planes in which the dipole moments all point in the same direction (on average). The direction is tilted by 45 • with respect to the average total polarization vector (e.g. along x axis). While the tilt angle is the same for all domains, the average orientations of the dipole moment projections onto the yz plane (perpendicular to the polarization) are located along the plane diagonals, and vary from one domain to the other by a rotation of 90 • . Figure 1 gives a view of this dipole arrangement projected on the xy plane. Such a structure had already been evidenced in previous canonical (NVT) MC simulations [13]. On the other hand, when the initial state was chosen to be a bcc or bct lattice the system generally evolved to a state compatible with that of lowest free energy as obtained by the NEW method. Examples will be given below.

The Jarzynski equality [START_REF] Jarzynski | [END_REF] < e -W/kBT >= e -(FF -FI )/kBT = e -∆F/kBT (1)

enables to calculate the free energy difference ∆F between equilibrium states I and F of a system by evaluating the average of e -W/kBT along a set of paths, i.e. a set of configurations relating the states I and F . W is the work done on the system in the transition I to F along one such path. As repeatedly proven and discussed in the literature, paths from I to F can be configurations obtained from hamiltonian dynamics [START_REF] Jarzynski | [END_REF][START_REF] Jarzynski | [END_REF]23] or from a Markov process satisfying detailed balance [24,25]. It is the last category of paths which has been used in the present work. More specifically, the initial state of a path is a system of N hard spheres at density ρ * I = ρ I σ 3 and dipole moment µ * located on a bcc lattice of volume V with periodic boundary conditions defined by the three vectors a, b and c corresponding to the columns of the matrix h I αβ with det(h I αβ ) = V . The system in this initial state is transformed into a state of same volume V specified by the matrix h F αβ . The matrix h F αβ is such that the shape of the volume V corresponds to N hard spheres located on a fcc, bct, bco or hcp lattice at density ρ * F = ρ * I . This transformation is achieved by a succession of matrices

h i αβ (i = 0, 1, ... , n) with h 0 αβ = h I αβ , h n αβ = h F αβ and h i αβ = h i-1 αβ + δ h i αβ .
The small random variations δ h i αβ are given by

δ h i αβ = (λ i -λ i-1 )(h F αβ -h I αβ ) (2)
where the variables λ i (0 < i < n) are a sequence of n increasing random numbers in the interval [0, 1] with λ 0 = 0 and λ n = 1.0. A sequence of δ h i αβ specifies the path of a transformation from I to F and satisfies the constraint i=1,..,n δ h

i αβ = h F αβ -h I αβ .
The work W done along a path is calculated as the sum of the works

w i corresponding to each variation δ h i αβ W = i=1,..,n w i = i=1,..,n V dr α,β τ i αβ (r) η i αβ (3)
where τ i αβ (r) is the tensor of constraint at position r of the system when the configuration of the N spheres changes volume according to the variation from h i αβ to h i+1 αβ and η i αβ is the deformation induced by this variation [26]. In the DHS system, the contribution of the dipolar interaction τ i,dd αβ (r) is

τ i,dd αβ (r) = 1 2 l=1,N δ(r -r i l ) k( =l)=1,N r i l,α F i l,k,β (4)
where F i l,k,β is the β component of the dipole-dipole force acting between the spheres l et k. The contribution of the hard core potential is obtained from the expression of the tensor of constaints derived by Allen [27,28]. The deformation η i αβ corresponds to the change in positions of the spheres entailed by the variation of the volume, i.e. for sphere l :

r i+1 l,α -r i l,α = γ,β h i+1 αγ hi γβ r i l,β -r i l,α , with hi αβ the inverse matrix of h i αβ .
It is identical at each point of the volume and given by

η i αβ = 1 2 γ δh i αγ hi γβ + γ δh i βγ hi γα + ν ( γ δh i αγ hi γν γ δh i νγ hi γβ ) . ( 5
)
The Jarzynski equality and Jensen's inequality [START_REF] Chandler | Introduction to Modern Statistical Mechanics[END_REF] imply the inequality < W > ≥ ∆F where < W > is the average W calculated for the ensemble of transformations from I to F . The average value, < W ′ >, of the work W ′ necessary to make the inverse transformation F to I, satisfies the inequality -< W ′ > ≤ ∆F . These two inequalities giving upper and lower bounds for ∆F , can be used to estimate ∆F in a solid-solid phase transition [START_REF] Hoy | [END_REF].

The matrices h I αβ and h F αβ specifying the initial and final states of the transformations from the bcc lattice to the bct, bco or fcc lattices are such that if r I 1 , r I 2 , ..., r I N are the lattice positions of the spheres on the bcc lattice, r F 1 , r F 2 , ..., r F N obtained from r F l,α = β,γ h F αγ hI γβ r I l,β , will be those of the spheres on the bct, bco or fcc lattices. In these cases the transformations from h I αβ to h F αβ can be made continuously, the volume V maintaining its parallelepipedic shape. The transformation from the bcc to the hcp lattice is also continuous but a relative shift of neighbouring planes, e. g. those perpendicular to the z-direction has to be supplemented with the volume change. The initial volume of the bcc lattice is chosen monoclinic and this shape is conserved in the evolution towards the hcp lattice. The lattice positions of the final hcp lattice are related to those of the bcc lattice by r F l,α = β,γ h F αγ hI γβ r I l,β + g δ α2 . The displacement g along the y direction generates the succession ...ABAB... of [001] planes characteristic of the hcp lattice. The value of g is alternatively equal to 0 (plane A) or b/3 (plane B) where b is the smallest distance along y between two spheres of a [001] plane of the final hcp lattice [31]. At each stage of the evolution of the volume a small random displacement δ g i is applied to the atoms in B planes subject to the constraint i=1,2,...,n δ g i = b/3. These displacements contribute to the values of w i taking into account, in Eq. ( 3), that the η i αβ differ depending on whether the spheres are located in A or B planes.

The calculation of W along the path from I to F comprises three stages. The first one consists in a simulation in the canonical ensemble to equilibrate the N spheres on the bcc lattice at the chosen values of ρ * et µ * . The second stage starts with the change of the matrix from h 0 to h 1 . This trial move is either accepted or rejected if there is overlap of the spheres. In the latter case canonical sampling is continued, a new trial transition from 0 to 1 is attempted and the process repeated until acceptance of the transition; w 1 is then calculated. This procedure is continued until the transition n -1 to n. The third stage is a simulation in the canonical ensemble at the volume defined by h n αβ = h F αβ to calculate the pair correlation functions in the final state F . This allows to check that the correlation functions are in agreement with those expected in the bct, bco, fcc or hcp phase for the DHS system at chosen ρ * and µ * .

The random variations δh i αβ are of the order (h F αβ -h I αβ )/n. To obtain an acceptance ratio of about 50% for the transitions i to i + 1, n must be of the order of 10 4 at the considered solid phase densities, canonical sampling between transition attempts being ∼ 5 MC trial moves per sphere. The numerical estimate of the hard core contribution to τ αβ (r) has been verified for µ * = 0 by calculating ∆F between dense fluid states of hard spheres for which the free energy is known. The Ewald summation method has been applied for the calculation of the energy and τ dd αβ (r) [32]. The values of N are 1024 for the calculation of the transitions from bcc to bct or bco, 1200 for those from bcc to hcp and 1372 for those from bcc to fcc.

Results

In all calculations of ∆F , the initial configuration was a bcc lattice of spheres completely polarized along the z axis (polarization 1). In the density region 1.06 ≤ ρ * ≤ 1.26, at µ * = 1, the equilibrium state associated with this initial configuration has zero polarization. At µ * = 2 and 3, it has a polarization between 0.75 and 0.87. These equilibrium states are the reference states used to estimate ∆F relative to the other solid phases. For each chosen thermodynamic state the average value of e -W/kBT has been calculated for 16 paths between the reference bcc state and the final state of symmetry bct, bco, fcc or hcp. Each path is defined by an initial equilibrium bcc configuration and a succession of different random variations δh i αβ . This number of different paths turned out to be sufficient because the values of W have a low statistical variance if amplitude and orientation of the polarizations are similar in the final states. The statement has been verified by calculating, at µ * = 3 and ρ * = 1.20, the values of W for 64 paths corresponding to a transition from bcc to bct. To test the method a comparison has been made with results for the phase diagram of the totally polarized hard sphere system [16,17]. In agreement with this study, the value of ∆F obtained from the Jarzynski equality shows that among the phases bct, hcp and hcp compressed along z by about 10%, the phase of bct symmetry is the most stable at µ * = 2 and ρ * = 1.20 and the compressed phase of hcp symmetry the most stable at µ * = 3 and ρ * = 1.06.

In spite of the constraints imposed by the excluded volume effects and the homogeneity of the deformations η αβ it is possible that the transition from the bcc phase to another phase is accompanied by creation of defects. The contribution of these possible defects to ∆F seems marginal. Indeed, as mentioned above, the variance of the values of ∆F is small for states with close final polarizations. Moreover, the pair correlation functions calculated in the final states are in excellent agreement with those obtained in canonical simulations at fixed density at identical thermodynamic states. An example of the agreement is given in Fig. 2 for the bcc-bct transition at µ * = 3 and ρ * = 1.16. At given density the bct lattice parameters are obtained from the bcc parameters by multiplying them by (3/2) 1/6 , (3/2) 1/6 and (2/3) 1/3 along x, y and z. For the bco lattices these factors are λ(3/2) 1/6 , (3/2) 1/6 /λ and (2/3) 1/3 along x, y and z with λ = 1.05, 1, 10, 1.15 and 1.20 in our simulations.

For the states µ * = 1 and ρ * ≥ 1.20 at which the polarizations vanish both in the initial and final states, the difference ∆F * = ∆F/k B T is always positive with a statistical error of the order of 0.1∆F * , indicating that bcc is the stable phase. This result is confirmed by the NpT simulations as shown in Fig. 3a which compares the pair distribution functions (p.d.f.) at ρ = 1.20 obtained in the isobaric ensemble, starting from an initial bcc lattice, with that for the final state of the NEW method. At densities ρ * ≤ 1.10 the bcc phase seems also to be the most stable phase but the values of ∆F * are small ∼ 0.1 with a statistical error of ≃ 0.03. In the domain 1.10 ≤ ρ * ≤ 1.20, the bco phase with λ ≃ 1.10 is marginally more stable than bcc, ∆F * ≃ -0.1. Again the NpT calculation transforms the initial bcc lattice into a state with p.d.f. in good agreement with that of the bco lattice with L y /L x = 1.10 identified by the NEW method as having the lowest free energy (cf. Fig. 3b).

When µ * increases the solid phases get polarized. At µ * = 2 and increasing densities 1.06 ≤ ρ * < 1.20 the most stable phase changes from bct to bco with values of λ between 1.05 and 1.15, the polarization being parallel to the z direction. For these bct and bco phases ∆F * is of the order -1.3 with a statistical error of about 0.02∆F * . At ρ * ≥ 1.20, the variance of ∆F * increases. As a matter of fact, for approximately 15% of the bco final states, i. e. 2 to 4 out of 16, the polarization is not along z but either along a diagonal of the parallelepipedic volume or a diagonal of the xz or yz planes and ∆F * is equal to ≃ -0.8. However, the states polarized along z have the lowest free energy ∆F * ∼ -1.1 and their equilibrium p.d.f's are in close agreement with those obtained from NpT simulations (cf. Fig. 3c).

If the final states are fcc or hcp, at µ * = 2, ∆F * varies from ∼ 0.3 for ρ * = 1.06 to ∼ -0.5 for ρ * = 1.20 then gets positive again ∼ 0.5 for ρ * = 1.26. As for the bcc and bco phases, the final states are generally polarized along z but other polarization states exist as well where the orientation lies in a lattice plane. Moreover, for the fcc lattice, there exists polarized states, where, as described above in Section II, the average orientation of dipole moments can differ from one domain of lattice planes to the other, while having the same tilt angle with respect to z. However the different polarization states of the fcc and hcp phases all have a free energy which is higher by ∼ 1 compared to the bct or bco phases.

The bct and bco phases with λ < 1.15 have nearly identical stability at µ * = 3 and ρ * ≤ 1.20, differences in ∆F * being close to the statistical error. The bco phase is marginally more stable at ρ * = 1.26 with λ ≃ 1.10. The polarization of the bct or bco phases of lowest free energy is parallel to z and equals ∼ 0.90 -0.95. A typical value of ∆F * between bcc and bct or bcc and bco is ∼ -4.3. For the hcp or fcc phases ∆F * is larger by 1.5 for all polarizations of the final state. Comparison of the p.d.f. at µ = 3 for ρ = 1.16 and ρ = 1.26 obtained either by NpT simulation or NEW method are shown in Figs. 4a and4b, respectively.

The lattices associated with the most stable phases of the final states have been called, somewhat improperly, with reference to the perfect lattices of hard spheres with volume compatible with the periodic boundary conditions defined by the matrices h F αβ . It is noteworthy that such an appellation is approximate. Upon increase of µ * the local order of the hard spheres is strongly modified by the dipolar interaction. Figures 5 and6 are examples of the reorganization of the local order at ρ * = 1.06 and 1.26. The correlation functions show that if at µ * = 1 or 2 the order is close to that of hard spheres (µ * = 0), it differs significantly at µ * = 3. This reorganization is made explicit in Fig. 6 for ρ * = 1.26. Although the positions of the peaks corresponding to the most numerous neighbours are preserved, peaks of other layers can overlap with these main peaks, or overlap between themselves; occasionally new peaks can appear. These rearrangements result from the open character of the local order in the bct and bco lattices; they are less important in the compact lattices like fcc where they vanish for densities ρ * ≥ 1.20 (cf. Figs. 7 and8). The partial delocalization of shells of neighbours thus does not allow a complete characterization of the cristalline order from the sole pair correlation functions. This is particularly clear for perfect bco lattices in which the number of nearest neighbors in successive shells are 2, 8, 2, 2, 4, 4, 8, 4, ..., respectively. At densities ρ * ∼ 1.20, the locations of the first and second shells and those of the seventh and eighth shells differ only by ∼ 0.002σ so that these shells merge due to the molecular motion.

When the final states are divided into states of different polarizations they will have different local hard sphere order. These differences in structural arrangement explain the difference in free energies of these final states. Thus, at µ * = 3.0 and ρ * = 1.26, the final states bct and bco, with polarization parallel to z have ∆F * ∼ -4.0 whereas those with polarization along a diagonal of the volume have a higher free energy ∆F * ∼ -1.0. Fig. 9 shows the change of local order entailed by these variations of free energy and polarization.

In the case where the probability distribution W is gaussian of width σ W , it is easy to show from the Jarzynski equality that [START_REF] Jarzynski | [END_REF]

∆F =< W > - σ 2 W 2k B T . (6)
Our results are compatible with this equality provided the polarizations of the final states are unique. When several states of polarization are possible, the distribution of values W has several peaks, is not gaussian and the equality ( 6) is not valid anymore.

Discussion

The results for the free energy differences calculated from Eq.( 1) are summarized in Fig. 10 which shows the domains of stability of the bcc, bct and bco phases with an approximate location of their limits as a function of µ * and ρ * . Only these phases are plotted since the NEW method results, as dicussed in Sect. 3, indicate that for this domain of µ * and ρ * values, the fcc and hcp phases are always unstable. They are in relative agreement with those given by the simulations in the isobaric ensemble. The solid phases identified by the NEW method as having the lowest free energy (among those considered) and those obtained by the simulations in the NpT ensemble are in good agreement for solid densities ρ * ≥ 1.16. In particular the local order obtained from the p.d.f. are quasi identical. Additional NpT calculations in this density region at µ * = 1.5 and µ * = 2.5 also predict structures compatible with those obtained by the NEW method simulations at µ * = 2 and µ * = 3. At the lower densities 1.06 ≤ ρ * ≤ 1.10 such an agreement could not be found. For µ * = 2 or 3, although the free energy calculations show that the bct phases or the bco phases with ratio L y /L x close to 1.0 are the most stable, the NpT simulations do not evolve spontaneously towards solid arrangements of this type. With different initial choices of the initial configuration the NpT simulations stabilize in a variety of states manifesting the possible existence of several metastable states close to the state of minimum free energy. By comparison of the pair distribution functions, it appears that these states obtained by the NpT simulations do not belong to the set of final states considered in the NEW computations; therefore the relative values of their free energy could not be estimated.

Conclusion

Use of the Jarzynski equality provides an alternative to the method of thermodynamic integration for the estimate of free energy differences between equilibrium thermodynamic states. The reliability of the calculation relies on the choice of the paths selected to calculate the average in Eq. ( 1) [33][34][35]. As mentioned above, by associating to the variation of h α,β the relaxation of the sphere configurations, it is possible to prevent overlap of the hard cores to interrupt the transition from one phase to the other. Moreover, if by construction of the transition path one has V F = V I , the volumes V i differ from V F only by 0.05% thus avoiding creation of defects. The choice of the paths seems justified by the facts that Eq. ( 6) is satisfied and that, taking into account the relaxation of the intermediate configurations and the number of values of n of order 10000 to 30000, the paths nearly satisfy the equality valid only for reversible equilibrium paths

F F -F I = i=1,n dr α,β < τ i αβ (r) > i η i α,β (r), (7)
where < ... > i denotes a canonical average at step i of a transition path.

As stated in the introduction the crossing through a first order phase transition precludes that Eq. ( 7) can be a strict equality, contrary to the Jarzynski equality. This illustrates again the relative advantages of the NEW and thermodynamic integration methods. Eq. ( 1) is always valid and needs only an initial reference state. It requires though to perform an average on a set of adequately chosen paths and, in the present work, the rather time comsuming computation of the hard sphere tensor of constraint. The thermodynamic integration method is simple to implement and uses only one path. However, to circumvent the first order phase transition problem, especially for solid phases, it needs to use multiple initial reference states the free energy of which must be known analytically. Such initial states, well adapted to hard sphere systems with rotational degrees of freedom, seem difficult to construct.

In conclusion, the solid phases of dipolar hard spheres estimated by free energy calculations to be the most stable ones are of the bcc, bct or bco type. For dipole moments µ * ≥ 2 these phases are polarized. For the bct and bco phases the direction of polarization is the direction along which the bcc unit cell is compressed (in our calculations the z direction). The pair distribution functions show, however, that the local arrangements of a particle and its nearest neighbours can differ notably according to the values of ρ * and µ * . A detailed and complete characterization of the unit cell of these phases remains difficult, in particular, due to the proximity of nearest neighbour shells and the small number of particles in these shells (cf. Sect. 3).

The NpT simulation results are in accord with the free energy results for thermodynamic states of fixed µ * and densities larger than roughly 1.10 and confirm in this domain the phase diagram given in Fig. 10. For the states µ * ≥ 2 and densities lower than approximately 1.10, although the stable phases are identified without ambiguity from the free energy calculations, this domain of the phase diagram has to be confirmed by a method different from that based on the NEW method.

Fig. 1 :Fig. 2 :Fig. 3 :Fig. 4 :

 1234 Fig. 1: Projection on the xy plane of the dipole moments for a fcc lattice with polarization along the x direction at ρ * = 1.26 and µ * = 3.0 (NVT simulation with N=1372).

Fig. 5 :Fig. 6 :Fig. 7 :Fig. 8 :

 5678 Fig.5: Pair distribution function g(r) for equilibrium states of a bct phase at ρ * = 1.06 and increasing values of µ * . For µ * = 2 and 3 the polarizations are parallel to z and equal to ∼ 0.9.

Fig. 9 :

 9 Fig. 9: Pair distribution functions g(r) for final equilibrium states of a bcc to bct phase transition corresponding to different orientations of the polarization at ρ * = 1.26 and µ * = 3.0. Solid line: polarization parallel to z axis; red line: polarization parallel to a diagonal of the parallelepipedic volume.

Fig. 10 :

 10 Fig. 10: Phase diagram obtained from free energy differences computed by the Jarzynski equality (NEW method). Two superimposed symbols for a state indicates close stability of the two phases.

  Fig. 2 279x361mm (300 x 300 DPI)

  Fig. 3 231x335mm (300 x 300 DPI)

  Fig. 4 231x331mm (300 x 300 DPI)

  Fig. 5 279x361mm (300 x 300 DPI)

  Fig. 6 279x361mm (300 x 300 DPI)

  Fig. 7 279x361mm (300 x 300 DPI)

  Fig. 8 279x361mm (300 x 300 DPI)

  Fig. 9 279x361mm (300 x 300 DPI)

  Fig. 10 279x361mm (300 x 300 DPI)


Olivier Beaumont

Hubert Larchevêque

Loris Marchal

Non Linear Divisible Loads: There is No Free Lunch

Keywords: Divisible Load Theory, MapReduce, Scheduling, Resource Allocation, Matrix Multiplication, Sorting

Divisible Load Theory (DLT) has received a lot of attention in the past decade. A divisible load is a perfect parallel task, that can be split arbitrarily and executed in parallel on a set of possibly heterogeneous resources. The success of DLT is strongly related to the existence of many optimal resource allocation and scheduling algorithms, what strongly differs from general scheduling theory. Moreover, recently, close relationships have been underlined between DLT, that provides a fruitful theoretical framework for scheduling jobs on heterogeneous platforms, and MapReduce, that provides a simple and efficient programming framework to deploy applications on large scale distributed platforms.

The success of both have suggested to extend their framework to non-linear complexity tasks. In this paper, we show that both DLT and MapReduce are better suited to workloads with linear complexity. In particular, we prove that divisible load theory cannot directly be applied to quadratic workloads, such as it has been proposed recently. We precisely state the limits for classical DLT studies and we review and propose solutions based on a careful preparation of the dataset and clever data partitioning algorithms. In particular, through simulations, we show the possible impact of this approach on the volume of communications generated by MapReduce, in the context of Matrix Multiplication and Outer Product algorithms.

I. INTRODUCTION

A. Divisible Load Theory and MapReduce

Scheduling the tasks of a parallel application on the resources of a distributed computing platform is critical for achieving high performance, and thus it is the target of many research projects, both in theoretical studies and in software developments.

As far as theory is concerned, the scheduling problem has been studied for a variety of application models, depending on the nature of the dependencies between tasks. In the case of general dependencies, the wellknown directed acyclic task graph model is used, and many scheduling heuristics have been developed [START_REF] Shirazi | Scheduling and load balancing in parallel and distributed systems[END_REF]. A slightly less general model corresponds to the parallel execution of several (usually many, so that scheduling concentrates on the steady state phase) pipelined applications that execute concurrently on a given platform, in particular for applications whose structure is a linear chain of tasks. Such applications are ubiquitous in streaming environments, as for instance video and audio encoding and decoding, DSP applications, image processing, and so on [START_REF] Beynon | Datacutter: Middleware for filtering very large scientific datasets on archival storage systems[END_REF], [START_REF] Hary | Precedence-constrained task allocation onto point-to-point networks for pipelined execution[END_REF], [START_REF] Wu | Supporting distributed application workflows in heterogeneous computing environments[END_REF], [START_REF] Wu | On optimal resource utilization for distributed remote visualization[END_REF]. An extreme (in terms of dependencies) application model is that of independent tasks with no task synchronizations and no inter-task communications. Applications conforming to this admittedly simple model arise in most fields of science and engineering. The number of tasks and the task sizes (i.e. their computational costs) may or may not be set in advance. In this case, the scheduling problems are akin to off-line and on-line bin-packing and a number of heuristics have been proposed in the literature (see [START_REF] Chrétienne | Scheduling Theory and its Applications[END_REF], [START_REF] Hochbaum | Approximation Algorithms for NP-hard Problems[END_REF], [START_REF] Epstein | Online bin packing with resource augmentation[END_REF] for surveys in off-line and on-line contexts).

Another flavor of the independent tasks model is one in which the number of tasks and the task sizes can be chosen arbitrarily. This corresponds to the case when the application consists of an amount (a large amount in general) of computations, or load, that can be divided into any number of independent pieces. This corresponds to a perfectly parallel job, whose sub-tasks can themselves be processed in parallel on any number of resources. This divisible load model has been widely studied, once Divisible Load Theory (DLT) has been popularized by the landmark book [START_REF] Bharadwaj | Scheduling Divisible Loads in Parallel and Distributed Systems[END_REF]. DLT provides a practical framework for the mapping of independent tasks onto heterogeneous platforms, and has been applied to a large spectrum of scientific problems, including Kalman filtering [START_REF] Sohn | Optimizing computing costs using divisible load analysis[END_REF], image processing [START_REF] Hamdi | Dynamic load balancing of data parallel applications on a distributed network[END_REF], video and multimedia broadcasting [START_REF] Altilar | An optimal scheduling algorithm for parallel video processing[END_REF], [START_REF]Optimal scheduling algorithms for communication constrained parallel processing[END_REF] database searching [START_REF] Drozdowski | Selected Problems of Scheduling Tasks in Multiprocessor Computer Systems[END_REF], [START_REF] Blazewicz | Divisible task scheduling -concept and verification[END_REF], and the processing of large distributed files [START_REF] Wang | Modeling communication pipeline latency[END_REF]. These applications are amenable to the simple master-worker programming model and thus can be easily implemented and deployed on computing platforms ranging from small commodity clusters to computational grids [START_REF]Computational Grids: Blueprint for a New Computing Infrastructure[END_REF] and Clouds [START_REF] Armbrust | A view of cloud computing[END_REF], [START_REF] Shih | Performance study of parallel programming on cloud computing environments using mapreduce[END_REF], [START_REF] Iyer | On handling large-scale polynomial multiplications in compute cloud environments using divisible load paradigm[END_REF]. From a theoretical standpoint, the success of the divisible load model is mostly due to its analytical tractability. Optimal algorithms and closed-form formulas exist for the simplest instances of the divisible load problem. This is in sharp contrast with the theory of task graph scheduling and streaming applications scheduling, which abounds in NP-completeness [START_REF] Garey | Computers and Intractability, a Guide to the Theory of NP-Completeness[END_REF] and in inapproximability results [START_REF] Ausiello | Complexity and Approximation[END_REF].

On the side of software development for scheduling parallel applications, the MapReduce framework [START_REF] Dean | Mapreduce: Simplified data processing on large clusters[END_REF], [START_REF] Zaharia | Improving mapreduce performance in heterogeneous environments[END_REF] has recently received at lot of attention. Indeed, the MapReduce model, which has been popularized by Google, is particularly well-suited to parallel processing of arbitrary data. Just as in the case of divisible load, a large computation is broken into small tasks that run in parallel on multiple machines (the case of a sequence of Map and Reduce operations has been studied in [START_REF] Berlińska | Scheduling divisible mapreduce computations[END_REF], and scales easily to very large clusters of inexpensive commodity computers. Hadoop [START_REF] White | Hadoop: The definitive guide[END_REF] is the most popular open-source implementation of the MapReduce framework, originally developed by Yahoo to manage jobs that produce hundreds of terabytes of data on thousands of cores. Examples of applications implemented with Hadoop can be found at http://wiki.apache.org/hadoop/ PoweredBy. A crucial feature of MapReduce is its inherent capability of handling hardware failures and processing capabilities heterogeneity, thus hiding this complexity to the programmer, by relying on on-demand allocations and a detection of nodes that perform poorly (in order to re-assign tasks that slow down the process).

The MapReduce programming model has first been introduced to deal with linear complexity tasks such as standard text processing operations [START_REF] Dean | Mapreduce: Simplified data processing on large clusters[END_REF], but it has been later extended to many other types of operations, such as linear algebra operations [START_REF] Seo | Hama: An efficient matrix computation with the mapreduce framework[END_REF]. In this case, the Map function, that is applied in parallel to every pair in the input, operates on a prepared dataset. For instance, in the case of a matrix product, one could imagine to have as input dataset all compatible pairs (a i,k ,b k,j ) for all n 3 possible values of i, j and k. In this case, the output of the Map operation would be a pair consisting of the value a i,k ⇥ b k,j and the key (i, j), so that all a i,k ⇥ b k,j , for 0  k  n -1 would be associated to the same key (i, j) and therefore to the same reducer, that would in turn be responsible for computing their sum. While allowing complex computations to run over MapReduce, such solutions lead to a large redundancy in data communication: if a given processor is responsible for computing many a i,k ⇥ b k,j products, it will receives as many values of a and b, even if some of them overlap.

B. Model

In this paper, the target computing platform is a heterogeneous master/worker platform, with p computing resources labeled P 1 ,P 2 ,...,P p . The master P 0 sends out chunks to workers over a network: we can think of a star-shaped network, with the master in the center. In order to concentrate on the difficulty introduced by the non-linearity of the cost, we consider the simplest communication model, where all communications between the master and the computing resources can take place in parallel, the speed of the communication between P 0 and P i being only limited by the incoming bandwidth of node P i . Unless stated otherwise, we assume that the computing platform is fully heterogeneous. The incoming bandwidth of processor is denoted 1/c i (so that c i is the time needed to send one unit of data to P i ) and its processing speed s i =1 /w i (so that w i is the time spent by P i to process a unit of computation).

In the literature, several models have been considered. The master processor can distribute the chunks to the computing resources in a single round, (also called single installment in [START_REF] Bharadwaj | Scheduling Divisible Loads in Parallel and Distributed Systems[END_REF]), or send the chunks to the computing resources in multiple rounds: the communications will be shorter (less latency) and pipelined, and the workers will be able to compute the current chunk while receiving data for the next one. In both cases, a slave processor can start processing tasks only once it has received the whole data.

Similarly, in order to concentrate on the influence of non-linearity, return communications [START_REF] Beaumont | Scheduling divisible loads with return messages on heterogeneous master-worker platforms[END_REF], [START_REF] Beaumont | Fifo scheduling of divisible loads with return messages under the one-port model[END_REF], [START_REF] Beaumont | Link-heterogeneity vs. node-heterogeneity in clusters[END_REF] will not be taken into account and we will consider the case of a single round of communications.

C. Outline

Thereby, there exists on the one side powerful theoretical models like DLT for scheduling jobs onto heterogeneous platforms, and on the other side powerful tools like MapReduce that allow to easily deploy applications on large scale distributed platforms. However, both are best suited to workloads with linear complexity: in this paper, we prove that divisible load theory cannot directly be applied to quadratic workloads in Section II, such as proposed in [START_REF] Hung | Scheduling nonlinear computational loads[END_REF], [START_REF]Scheduling nonlinear computational loads[END_REF], [START_REF] Suresh | An analytical solution for scheduling nonlinear divisible loads for a single level tree network with a collective communication model[END_REF], [START_REF] Suresh | Scheduling nonlinear divisible loads in a single level tree network[END_REF], [START_REF] Suresh | Scheduling second-order computational load in master-slave paradigm[END_REF]. We precisely state the limits for classical DLT studies in Section III. Then, we review and propose solutions based on a careful preparation of the dataset, and study how this approach could be applied in the context of MapReduce operations in Section IV.

Recently, several papers [START_REF] Hung | Scheduling nonlinear computational loads[END_REF], [START_REF]Scheduling nonlinear computational loads[END_REF], [START_REF] Suresh | An analytical solution for scheduling nonlinear divisible loads for a single level tree network with a collective communication model[END_REF], [START_REF] Suresh | Scheduling nonlinear divisible loads in a single level tree network[END_REF], [START_REF] Suresh | Scheduling second-order computational load in master-slave paradigm[END_REF] have studied the case on non-linear divisible scheduling. For instance, let us consider that a non-linear (say N α ,α>1) cost operation is to be performed on a list of N elements. In order to analyze the impact of non-linearity and to assess the limits of the approach followed in these papers, we will concentrate in this section on fully homogeneous platforms, where all communication and computing resources have the same capabilities.

Each computing resource P i , 1  i  P has a (relative) computing power 1 w and it is associated to a bandwidth 1 c . Thus, it will take c • X time units to transmit X data units to P i and w • X α to execute these X units of load on P i . In this model, the direct translation of classical linear DLT problem to the nonlinear case DLT that is proposed in [START_REF] Hung | Scheduling nonlinear computational loads[END_REF], [START_REF]Scheduling nonlinear computational loads[END_REF] consists in finding the amount of load that should be allocated to each processor.

In general, complex platforms are considered, with fully heterogeneous computing and communication platforms, and some flavor of one-port model, what makes the resolution of above problem difficult and requires sophisticated techniques (see for instance the solutions proposed in [START_REF] Suresh | An analytical solution for scheduling nonlinear divisible loads for a single level tree network with a collective communication model[END_REF], [START_REF] Suresh | Scheduling nonlinear divisible loads in a single level tree network[END_REF], [START_REF] Suresh | Scheduling second-order computational load in master-slave paradigm[END_REF]) in order to determine in which order processing resources should be fed by the master node and how much data to process they should be given. In our very simplistic homogeneous model, with parallel communications and no return communication of the results, the above problem becomes trivial: the ordering is not important since all the processors are identical and therefore, in the optimal strategy, each P i receives N P data elements in time N P c and starts processing them immediately until time

N P c + ✓ N P ◆ α w.
In general [START_REF] Suresh | An analytical solution for scheduling nonlinear divisible loads for a single level tree network with a collective communication model[END_REF], [START_REF] Suresh | Scheduling nonlinear divisible loads in a single level tree network[END_REF], [START_REF] Suresh | Scheduling second-order computational load in master-slave paradigm[END_REF], for more heterogeneous platforms and more sophisticated communication models, obtaining the optimal allocation and a closed formula is impossible. Nevertheless, the importance of this issue is not crucial. In fact, the main problem with this approach is that when P is large, the part of the computations that is processed during this phase is negligible with respect to the overall work to be processed. Indeed, the overall work W is given by W = N α and the overall work W Partial performed by all p processors during this first phase is given by

W Partial = P ✓ N P ◆ α = N α P α-1 , so that W-W Partial W =1- 1 P α-1
that tends toward 0 when P becomes large ! Therefore, the difficult optimization problem solved in [START_REF] Suresh | An analytical solution for scheduling nonlinear divisible loads for a single level tree network with a collective communication model[END_REF], [START_REF] Suresh | Scheduling nonlinear divisible loads in a single level tree network[END_REF], [START_REF] Suresh | Scheduling second-order computational load in master-slave paradigm[END_REF] has in practice no influence on the overall computation time, since asymptotically (when P becomes large) all the work remains to be done after this first phase. Above results shows the intrinsic linear complexity of the problems that are divisible. Indeed, a divisible load can be arbitrarily split between any number of processors (say N ) in small pieces (of size 1 in this case) that can be performed independently (in time c + w), so that no dependencies should exist between data and the overall complexity is necessarily linear in the size of the data ((c + w)N in this setting).

Nevertheless, this does not mean that only linear cost complexity tasks only can be processed as divisible load tasks or using MapReduce. Indeed, the possibility remains (i) either to modify the initial data, such as proposed in [START_REF] He | Mars: a mapreduce framework on graphics processors[END_REF], [START_REF] Fatahalian | Sequoia: programming the memory hierarchy[END_REF], [START_REF] Seo | Hama: An efficient matrix computation with the mapreduce framework[END_REF] for matrix multiplication (in this case the initial N 2 size data is transformed into a N 3 size data by replicating matrix blocks before applying a MapReduce operation) (ii) or to decompose the overall operation using a long sequence of MapReduce operations, such as proposed in [START_REF] Berlińska | Scheduling divisible mapreduce computations[END_REF]. In the remaining of this paper, we will concentrate on the first approach, that consists in expressing a problem that contains data dependencies into a (larger) problem, where data dependencies have been removed.

III. DIVISIBLE LOAD THEORY FOR ALMOST LINEAR WORKLOADS

A. Sorting with Homogeneous Computing Resources

As we have seen in the previous section, DLT cannot in general be applied to workloads with super-linear complexity N α , α>1 unless the size of the initial data is increased or if several DLT operations are applied. However, there are some intermediate cases where the complexity is close to N and thus where the workload can be seen as almost divisible.

In this section, we consider the problem of sorting N numbers and we propose a parallel algorithm based on DLT to distribute this computation onto an homogeneous platform. Indeed, since the work required by sorting is W = N log N and if the initial dataset is split into p lists of size N p , the work produced when sorting the p lists in parallel in the DLT phase is given by

W Partial = p ✓ N p ◆ log ✓ N p ◆ = N log N -N log p, so that W-W Partial W = log p log N ,
that is arbitrarily close to 0 for large values of N . Therefore, contrarily to what happens is the case of tasks whose complexity is N α , in the case of sorting, almost all the work can be expressed as a divisible load task and sorting is likely to be amenable to DLT. Nevertheless, applying directly partial sorts to the p lists would not lead to a fully sorted result, and a preprocessing of the initial list is needed. More precisely, we rely on the sample sort introduced and analyzed in [START_REF] Frazer | Samplesort: A sampling approach to minimal storage tree sorting[END_REF], [START_REF] Li | Parallel sorting by over partitioning[END_REF]. The sample sort is a randomized sort, that relies on the use of a random number generator. The parallel running time is almost independent of the input distribution of keys and all complexity results will be given with high probability. The algorithm proceeds in three phases, depicted on Figure 1 (where red lists are unsorted and blue lists are sorted).

• Step 1: A set of p -1 splitter keys are picked and then sorted to partition the linear order of key values into p buckets.

• Step 2: Based on their values, the keys are sent to the appropriate bucket, where the ith bucket is stored in the ith processor.

• Step 3: The keys are sorted within each bucket (using 1 processor per bucket). Clearly, due to randomization, the buckets do not typically have exactly equal size and oversampling is used to reduce the ratio between the size of the largest bucket and its expected size N p . Using an oversampling ratio of s, a sample of sp keys are chosen at random, this sample is sorted, and then the p-1 splitters are selected by taking those keys in the sample that have ranks s, 2s, 3s, ..., (p -1)s.

Therefore, the time taken by Step 1 (using oversampling) is sp log(sp), i.e. the necessary time to sort (on the master processor) the sample of size sp. The cost of Step 2 (on the master processor) is given by N log p, since it involves, for each element to perform a binary search in order to determine its bucket.

At last, let us determine the time taken by Step 3. First, let us apply Theorem B.4 proved in [START_REF] Blelloch | A comparison of sorting algorithms for the connection machine cm-2[END_REF], with α =1+( 1 / log N ) N . The time taken by Step 3 is given by (MaxSize log MaxSize) so that, with high probability, it is bounded by

N p log N + o (N log N ) .
Therefore, if we choose s = (log N ) 2 as oversampling ratio, the time required by the preprocessing of Step 1 and Step 2 is dominated by N log p and the overall execution time is given by the time taken by Step 3, i.e. N p log N + o(N log N ) and is therefore optimal on p processors with high probability.

Therefore, in the case of sorting, it is possible, by introducing a preprocessing phase on the initial data (but keeping the same data size), to reduce the high cost operation to a fully Divisible Load Task and therefore, in the case of sorting, optimizing the data distribution phase to slave processors under more complicated communication models that the one considered in this paper, is meaningful.

B. Generalization to Heterogeneous Processing Resources

Extending above result to heterogeneous computing resources is not difficult. Indeed, let us denote as previously by 1 wi the processing power of computing resource P i , 1  i  P , so that sorting N i elements sequentially on P i takes w i N i log N i . Then, after oversampling with ratio s = (log N ) 2 , we assign to the ith bucket elements whose value is in the interval "$ P i-1

0 1/w k P P 0 1/w k (P -1) % s, $ P i 0 1/w k P P 0 1/w k (P -1) % s -1 # .
Therefore, processor P i receives (see Theorem B.4 in [START_REF] Blelloch | A comparison of sorting algorithms for the connection machine cm-2[END_REF], taking again α =1+( 1 / log N )

1
3 and s = log 2 N )) a list of at most 1/wi P P 0 1/w k elements (asymptotically when N becomes large) with high probability, and the load is well balanced with high probability between the processors (asymptotically when N becomes large).

Therefore, sorting is amenable to a divisible task at the price a preprocessing phase, even in presence of heterogeneous processors.

IV. SCHEDULING NON-LINEAR WORKLOADS ON HETEROGENEOUS PLATFORMS

As shown in Section II, Divisible Load Theory cannot be used to schedule non-linear workloads, i.e. workloads with a complexity of order O(N α ), with α>1 for a size N of the data. For these workloads, dividing the computation into small chunks (whose size sums up to N ) does not enable to perform enough work, and the same chunk of data is needed for several chunks of computations. For example, when multiplying two N ⇥ N matrices A and B, the same element a i,j needs to be multiplied to N elements of B, thus contributing to several chunks of computation.

One solution to overcome this problem is to introduce data redundancy, i.e. to replicate the data for each chunk of computation. For matrix multiplication, the a i,j element will be replicated for each computational chunk where a i,j is involved. Several distributions of the data are possible and have been implemented using MapReduce, such as row/column distribution or block distribution [START_REF] Seo | Hama: An efficient matrix computation with the mapreduce framework[END_REF], [START_REF] He | Mars: a mapreduce framework on graphics processors[END_REF]. Whatever the distribution, in MapReduce, the load-balancing is achieved by splitting the workloads in many tasks, which are then scattered across the platform. The fastest processor (or the one with smallest external load) gets more chunks than the others, so that all processors finish their share approximately at the same time (in Hadoop [START_REF] White | Hadoop: The definitive guide[END_REF], some tasks are themselves replicated at the end of the computations to minimize execution discrepancy). Assuming large and fast communication capacities, this method produces an almost perfect load-balancing. However, because of data replication, it may lead to a large communication overhead. Thus, communication links may become bottleneck resources if the replication ratio is large.

In this section, we study how to reduce the communication overhead for non-linear workloads on heterogeneous platforms. Our method relies on understanding the underlying data structure and the structure of dependencies between tasks, and then proposing a well-suited distribution of the data. Our objective is to achieve a perfect load-balancing of the workload and simultaneously to minimize the amount of communications.

We first present in Section IV-A a data distribution scheme for a basic operation, the outer product of two vectors (whose complexity is of order N 2 and data size is of order N ), before extending it in Section IV-B to the matrix multiplication (whose complexity is of order N 3 and data size is of order N 2 ). For each problem, we compare our method to classical block distributions used in MapReduce implementations, both theoretically and through simulations.

A. Outer-product: 2D data distribution

We consider the problem of computing the outerproduct a T ⇥ b of two (large) vectors a and b (of size N ), what requires the computation of all a i ⇥ b j , 8 1  i, j  N (see Figure 2(a)). As stated above, we target an heterogeneous computing platforms and we denote by s i =1 /w i the processing speed of P i . Let us also introduce the normalized processing speed of P i x i = s i / P k s k so that P i x i =1 . At last, let us assume without loss of generality that processors are sorted by non-decreasing processing speed:

s 1  s 2 •••s p .
Our objective is to minimize the overall amount of communications, i.e. the total amount of data send by the master initially holding the data (or equivalently by the set of devices holding it since we are interested in the overall volume only), under the constraint that a perfect load-balancing should be achieved among resources allocated to the outer product computation. Indeed, due to data dependencies, if we were to minimize communications without this load-balancing constraint, the optimal (and inefficient) solution would consist in making use of a single computing resource.

In what follows, we compare two approaches, Homogeneous Blocks and Heterogeneous Blocks.

1) Homogeneous Blocks: the first approach is based on a classical block distribution of the data, where the computational domain is first split into a (large) number of homogeneous chunks. Then, a demand driven model is used, where processors ask for new tasks as soon as they end processing one. With this approach, load is well balanced and faster processors get more blocks than slower one. This typically corresponds to the data distribution schemes used by MapReduce implementations. 2) Heterogeneous Blocks: The second approach consists in taking into account the heterogeneity of the processors when splitting the computational domain. Among the partitions where processors receive a fraction of the domain proportional to their speed, we search for the one that minimizes the amount of exchanged data. 1) Homogeneous Blocks approach: Let us first consider the Homogeneous Blocks approach, where the computation domain is split into squares of size D ⇥ D. Each chunk of computation consists in the outer-product of two chunks of data of size D from vectors a and b: (a i ,...,a i+D-1 ) T ⇥ (b j ,...,b j+D ). We choose to partition the computation domain into square blocks in order to minimize the communication costs: for a given computation size (D 2 ), the square is the shape that minimize the data size (2D). Heterogeneity is handled using dynamic load balancing: faster processors will get more blocks than slower ones. In order to minimize the communication cost, we choose the size of blocks so as to send a single chunk to the slowest processor. Since its relative processing speed is x 1 , the size of the atomic block will be D 2 = x 1 N 2 , so that D = p x 1 N (let us assume that N is large so that we can assume that this value is an integer). The total number of blocks is thus N/ p x 1 N 2 =1 /x 1 . Using this atomic block size, if the demand driven scheme achieves perfect load balancing, processor P i receives a number of blocks n i proportional to its computing speed

n i = x i x 1 = s i s 1 .
Again, let us assume for now that all these quantities are integer. The number of chunks distributed among processors is therefore

P i x i /x 1 =1/x 1 ,
what ensures that all blocks are processed.

The total amount of communications Comm hom generated by the Homogeneous Blocks approach is the number of blocks times the input data for a block (2D).

Comm hom =1/x 1 ⇥ 2N p x 1 =2N s P i s i s 1 .
2) Heterogeneous Blocks approach: The main drawback of the Homogeneous Blocks approach is clear in the context of strongly heterogeneous resources. When the ratio between the smallest and the largest computing speeds is large, the fastest processor will get a large number of (small) square chunks. For such a processor, the ratio between communications and computations is far from being optimal. Indeed, if these small square chunks could be grouped into a larger square chunk, data reuse would be better and the required volume of communications would be smaller.

With the Heterogeneous Blocks approach, a unique chunk is sent to each processor in order to avoid unnecessary data redundancy. We build upon a partitioning algorithm presented in [START_REF] Beaumont | Partitioning a square into rectangles: Np-completeness and approximation algorithms[END_REF]. In this paper, the problem of partitioning a 1 ⇥ 1 square into p rectangles of prescribed area a 1 ,a 2 ,...,a p is addressed. Two different objectives related to the perimeter of the rectangles in the partition are considered: minimizing the maximum half-perimeter (PERI-MAX) or the sum of the halfperimeters (PERI-SUM).

Our objective is to partition the square computational domain of the a i ⇥ b j into p rectangles, whose area are proportional to their relative computation speed (so as to enforce an optimal load balancing), and to minimize the overall volume of communications. A processor will be given a chunk of k consecutive values of a (a i ,...,a i+k ) and l values of b (b j ,...,a j+l ). The amount of communication for this processor is given by k + l, i.e. the half-perimeter of its rectangular chunk. Moreover, chunks must partition the whole domain without overlap. Using scaled computational speeds a i = x i , our problem is equivalent to minimizing the sum of the half-perimeters when partitioning the unit square into rectangles of prescribed area, i.e. the PERI-SUM problem.

In [START_REF] Beaumont | Partitioning a square into rectangles: Np-completeness and approximation algorithms[END_REF], several column-based algorithms are proposed: the square domain is first split into columns that are then divided into rectangles. In particular, a column-based partitioning algorithm for the PERI-SUM problem is proven to have a guaranteed performance. The sum of the half-perimeters Ĉ given by the algorithm is such that

Ĉ  1+ 5 4 LBComm where LBComm =2 p X i=1 p a i ,
where LBComm is a straightforward lower bound on the sum of the half-perimeters, which is larger than 2. Thus,

Ĉ  7 4
LBComm.

Note that as soon as there are enough processors, LBComm 2 so that the approximation ratio gets asymptotically close to 5/4.

Using this partitioning, scaled to the N ⇥ N computational domain, the total amount of communications Comm het can be bounded as

Comm het  7N 2 p X i=1 p x i = 7N 2 P p i=1 p s i p P p i=1 s i .
3) Comparison of Block Homogeneous and Block Heterogeneous Approaches: Using previous analysis, we can bound the ratio ρ between the amounts of communication generated by both Block Homogeneous and Block Heterogeneous approaches. When the platform is fully homogeneous, i.e. when all processors have the same computation speed s i = s 1 , both approaches provide the same solution, and our analysis simply states that ρ ≥ 4/7. However, the Block Heterogeneous approach is interesting as the platform goes heterogeneous.

ρ = Comm hom Comm het ≥ 2N q P i si s1 7N 2 P p i=1 √ si p P p i=1 si = 4 7 ⇥ P i s i p s 1 P i p s i .
A typical situation is depicted on Figure 2(b) and Figure 2(c), where a basic task corresponds to a small red matrix block. In the case of the Homogeneous Block approach, tasks are requested on demand and therefore, the memory footprint (data needed by a processor of relative speed 12) on a and b vectors will be high, whereas with the Heterogeneous Block approach, the volume of necessary data is highly reduced.

Consider for example the case where the first half of the platform is built from slow nodes (speed s 1 ) while the second one is built from nodes that are k times faster (speed ks 1 ). Then, our analysis shows that ρ ≥ 1+k

1+ √ k ≥ p k -1.

B. 3D data distribution: matrix multiplication

Let us start by briefly recalling the classical parallel matrix multiplication implementations. The whole set of computations can be seen as a 3D cube where element (i, k, j) corresponds to the basic operation a i,k b k,j . At the notable exception of recently introduced 2.5D schemes [START_REF] Solomonik | Communication-optimal parallel 2.5d matrix multiplication and lu factorization algorithms[END_REF], all implementations (see [START_REF] Goto | Anatomy of highperformance matrix multiplication[END_REF] for a recent survey), including those implemented with MapReduce [START_REF] He | Mars: a mapreduce framework on graphics processors[END_REF], [START_REF] Seo | Hama: An efficient matrix computation with the mapreduce framework[END_REF] or designed for GPUs [START_REF] Matsumoto | Multi-level optimization of matrix multiplication for GPU-equipped systems[END_REF] are based on the ScaLAPACK algorithm [START_REF] Blackford | ScaLAPACK user's guide[END_REF], that uses the outer product described in Section IV-A as building block. For the sake of simplicity, we will concentrate on the case of square matrices only. In that case, all 3 matrices (A, B and C = A ⇥ B) share the same layout, i.e. for all i, j, the same processor is responsible for storing A i,j ,B i,j and C i,j . Then, at each step k, any processor that holds some part of the kth row of A broadcasts it to all the processors of its column and any processor that holds some part of the kth column of B broadcasts it to all the processors of its row (see Figure 3). Of course, actual implementations of this algorithm use a blocked version of this scheme, and a level of virtualization is added. Indeed, since the number of blocks is much larger than the number of processors, blocks are scattered in a cyclic fashion along both grid dimensions, so that each processor is responsible for updating several blocks at each step of the algorithm. Thus, the volume of communication induced by the matrix multiplication algorithm is exactly proportional to the sum of the perimeters of the rectangles assigned to the processors, and therefore, the ratio proved in Section IV-A is valid between heterogeneity aware implementations based on Heterogeneous Block distributions and MapReduce implementations [START_REF] He | Mars: a mapreduce framework on graphics processors[END_REF], [START_REF] Seo | Hama: An efficient matrix computation with the mapreduce framework[END_REF]. In Section IV-C, we experimentally analyze this ratio through extensive simulations corresponding to several heterogeneity profiles.

C. Experimental evaluation

To assess the quality of the two data distribution policies presented above, we present here an experimental evaluation based on simulations. Both policies are compared to the lower bound on the communication LBComm introduced above: this bound corresponds to allocating to each processor a square whose area corresponds exactly to its relative processing speed, i.e. each processor P i is given a N p x i ⇥ N p x i square.

Considering that all obtained quantities are integer, the following lower bound holds for the total number of necessary communications LBComm =2N P

i p x i = 2N P i √ si p P i si .
In the Block Homogeneous strategy, we make a strong assumption: the number of chunks to give to processor P i is given by s i /s 1 , that is supposed to be an integer. In fact, these numbers have to be rounded to integers, thus leading to a possibly prohibitive load imbalance. Therefore, we propose a more realistic strategy, that splits the chunks into smaller blocks in order to avoid a large load imbalance. Let us first define the load imbalance of a given load distribution as e = t max -t min t min , where t max (respectively t min ) is the largest (resp. smallest) computation time in the platform. In the Comm hom strategy, the chunk size is chosen as large as possible to avoid unnecessary data redundancy, what may lead to a large load imbalance. To avoid this, we introduce the Comm hom/k strategy, that divides the block-size by k for increasing values of k until an acceptable load-balance is reached. In our simulations, the stopping criterion for this process is when e  1%. This strategy is expected to lead to a larger amount communications, at the benefit of a better load balancing. We perform simulations for a number of processors varying from 10 to 100. For each simulation, we generate the processing speeds using three different policies: the processing speeds either (i) are homogeneous, (ii) follow a uniform distribution in the range [START_REF] Shirazi | Scheduling and load balancing in parallel and distributed systems[END_REF]100] or (iii) follow a log-normal distribution with parameters µ =0and σ =1 . Figures 4(a), 4(b) and 4(c) present the results of these simulations. We compute the amount of communication induced by each strategy for a large matrix and we plot the ratio with the lower bound on communication. Each point represent the average ratio of a given strategy for 100 simulations with random parameters and error bars illustrate the standard deviation of this ratio.

As expected, the Block Homogeneous Comm hom strategy performs very well in an homogeneous setting: each processor gets a square corresponding to its share, so that Comm hom/k does not increase the number of chunks. Comm het requires more communications, but the increase is usually as small as 1% of the lower bound. However, as soon as computing speeds get heterogeneous, Comm hom and Comm hom/k experience a large increase in the volume of communications: with 100 processors, the more realistic Comm hom/k strategy leads to a communication amount which is 15 to 30 times the lower bound (depending on the random distribution of computing speeds). On the contrary, the Comm het strategy never requires more than 2% more than the lower bound.

Interestingly, we can notice that Comm het is much better than its theoretical guarantee: in the previous analysis, we proved that it is a 7/4-approximation compared to the lower bound, but in all our experiments, it is always within 2% of this bound.

V. C ONCLUSION

Recently, several papers have considered the case on non-linear divisible scheduling. We prove in this paper that tasks whose complexity is of order N α , for α>1 cannot be considered as divisible tasks, except if data is first replicated. The same applies to MapReduce jobs for non linear complexity tasks, that require a specific preparation of the data set. We have considered and analyzed the case of matrix multiplication (or equivalently outer product). On the one hand, we prove that if all computing resources have the same capacities, then classical implementations achieve very good results, even for non linear complexity tasks.

On the other hand, these experiments prove that in the context of heterogeneous computing resources, and when the complexity of the underlying task is not linear, taking explicitly heterogeneity into account when partitioning data and building tasks is crucial in order to minimize the overall communication volume, that can be easily reduced by a factor of 15 to 30. Nevertheless, at present time, by construction, MapReduce implementations are not aware of the execution platform and cannot therefore use an adapted data layout. Without changing the programming framework, whose simplicity is essential, adding directives in order to declare affinities between tasks and data could be partially solve this problem. For instance, in the context of the outer product or the matrix multiplication, favoring among all available tasks on the master those that share blocks with data already stored on a slave processor in the demand-driven process would improve the results. Therefore, we believe that this paper provides a sound theoretical background for studying data distribution and partitioning algorithms in the context of the execution of non-linear complexity tasks on heterogeneous platforms and opens many practical perspectives, by proving the interest of proposing new mechanisms in MapReduce to take into account affinities between tasks and data.

1 3Figure 1 .

 11 Figure 1. Sample sort with p =4processors and s =4oversampling ratio

  Using Heterogeneous Block approach.

Figure 2 .

 2 Figure 2. Data sent to a given processor for both implementation of outer-product.

Figure 3 .

 3 Figure 3. Matrix Multiplication algorithm based on the Outer Product

  amount to the lower bound Comm het Comm hom Comm hom/k (c) Computation speeds following log-normal random distribution.


Anthony Meilland

Philippe Stee

Alain Spang

Fabien Malbet

Fabrizio Massi

Dieter Schertl

The binary Be star δ Scorpii at high spectral and spatial resolution : II The circumstellar disk evolution after the periastron

Keywords: high angular resolution -Techniques: interferometric -Stars: emission-line, Be -Stars: winds, outflows -Stars: individual (δ Sco) -Stars: circumstellar matter Stars: binary

Introduction

We initiated an interferometric follow-up of the binary Be star δ Scorpii which has been regularly observed with the VLTI/AMBER (Petrov et al. 2007) and CHARA/VEGA (Mourard et al. 2009) instruments since 2007. In Meilland et al. (2011), hereafter Paper I, we managed to refine δ Scorpii's binary orbit and the orbital period was the only orbital element which was significantly different from the one determined by Tango et al. (2009). Consequently, we determined that the next periastron should take place around 2011 July 5 (± 4 days). This result was also confirmed by Tycner et al. (2011). Thanks to the combination of high spectral and spatial resolution we constrained the geometry and the kinematics of δ Scorpii's circumstellar environment. We confirmed that the line emission originates from a dense equatorial disk fully dominated by rotation. We found that the rotation appears to be Keplerian, with an inner Send offprint requests to: ame@oca.eu ⋆ Based on observations made with ESO telecopes at La Silla Paranal Observatory under GTO programme IDs 089.D-0105(A) and 089.D-0105(B) ⋆⋆ Fig. 3 is only available in electronic form via www.edpsciences.org boundary (photosphere/disk interface) rotating close to the critical velocity (V c ). The expansion was found to be approximately 0.2 km s -1 .

With the measured vsini of 175 km s -1 (Berbacca & Perinotto 1970), and the determined inclination angle of 28 ± 8 o we inferred that the star rotates between 60% and 105% of its critical velocity. Moreover, taking into account a possible underestimation of the vsini because of gravity darkening (Townsend et al. 2004), the lower limit on the stellar rotation would rise to about 70 % V c . Finally, we measured the extension in the Hα line, i.e., 4.8 ± 1.5 mas, which is about the same as the binary separation at the periastron, i.e., 5.9 mas. Thus, we suggested that a strong interaction between the already ejected matter and the companion should occur during the periastron passage. This letter is organized as follows: in Sect. 2 we briefly summarize other recent observations of δ Scorpii; in Sect. 3 we present our VLTI/AMBER observations and the data reduction. In Sect. 4, the data are analyzed and compared to those acquired in April 2010, one year before the periastron, and finally, a short conclusion is drawn in Sect. 5.

Recent observations of δ Scorpii

the star δ Scorpii (7 Scorpii, HR 5953, HD 143275) was observed by Millan-Gabet et al. (2010) using the CHARA/MIRC instrument in the H band. They resolved an elongated disk with a Gaussian FWHM of 1.18 x 0.91 mas. In the K band the disk was only marginally spatially resolved and they estimated an uncertain continuum disk FWHM of 0.7 ± 0.3 mas. When taking into account the continuum/line flux ratio they obtained much larger sizes for the line emission regions, i.e., 2.2 ± 0.4 mas for He i and 1.9 ± 0.3 mas for Brγ. In Paper I, we resolved the circumstellar disk in the Hα (FWHM= 4.8 ± 1.5 mas), Brγ (FWHM=2.9±0.5 mas), and the 2.06 µm He I (FWHM= 2.4 ± 0.3 mas) lines, as well as in the K-band continuum (FWHM ∼ 2.4 mas).

In a recent paper Che et al. (2012) imaged the disk distortion near the periastron. They revised the binary orbit and studied the primary disk properties during the periastron with NPOI and CHARA/MIRC interferometers. They were able to fit a global model where the secondary follows a Keplerian orbit, as already mentioned in Paper I, and the disk properties were found to be stable through all their observations. The estimated the new periastron passage to be UT 2011 July 3 07:00 ± 4:30 which is within the error bar of the prediction made in Paper I. They also estimated the mass of the secondary, based on the revised binary orbit and radial velocity (RV) measurements from Miroshnichenko et al. (2001), to be ∼ 6 M ⊙ with the primary mass 13.9 M ⊙ estimated from v sin i, apparent Teff, and V-band photometry measurements. The mutual angles between the disk mid-plane and orbital plane were estimated to be either ∼ 119 • or ∼ 171 • . From the fitted global model, they also found that ∼ 3% of the H-band flux comes from a fully resolved envelope. The primary disk was found to be mainly symmetric and stable, contributing to 71.4 ± 2.7% of the total H-band flux. This implies a quiescent inner disk and no ongoing material outflow after the periastron which is very different from what we obtained in the K band.

On the other hand, Rivinius et al. (2012) conducted a spectroscopic study of δ Scorpii emission lines during the periastron. They found that Hα and Hβ profiles were strongly modified at the periastron. The major modification was that the Violet/Red (V/R) ratio went from V/R>1 to V/R<1 in only a few days. They argued that this was due to the direct tidal effect of the companion on the primary disk.

We note that we looked at many Hα and Hβ observations made by amateurs available in the BeSS database1 and that the line profiles did not change significantly from the state detected by Rivinius et al. (2012) after the periastron.

Observations and data reduction

The star δ Scorpii was observed in June 2012 during two half-nights with the VLTI 1.6m auxiliary telescopes (AT) using the AMBER beam recombiner (Petrov et al. 2007). The observations were carried out using high spectral resolution mode (R=12000) centered on either Brγ or He i (2.06 µm) for each of the two half-nights. In order to enhance the signalover-noise ratio (S/N) by enabling longer exposure time, the fringe tracker FINITO was also used. The star HD 139663 (D = 2.133 ± 0.147 mas, [START_REF] Lafrasse | JMMC Stellar Diameters Catalogue[END_REF]) was used as an interferometric calibrator.

The data were reduced using the standard AMBER data reduction software amdlib v3.0.3b1 (see Tatulli et al. 2007 andChelli et al. 2009 for detailed information on the AMBER data reduction). The average-raw complex visibility and closure phase was determined using the standard method, i.e., keeping the 20 % of the frames with the higher S/N ratio. The interferometric calibration was then done using custom-made scripts2 described in Millour et al. (2008). The (u,v) plan coverage for all the observations is plotted in Fig. 1.

The disk geometry and kinematics

Qualitative comparison of the 2010 and 2012 data

The Brγ and He i absolute visibility show a strong variability with an amplitude of about 0.2. Assuming that the binary separation at the epoch of the observations should be 90.8 mas with a position angle (PA) of 156.08 o , according to Paper I orbital parameters, these variations are compatible both in terms of amplitude and frequency with the binary modulation. However, the accuracy of the absolute visibility of VLTI/AMBER data taken in HR mode is not sufficient to significantly refine the orbital parameters.

The differential visibilities and phases, closure phase, and line profile from the 2010 Brγ observation are plotted in Fig. 2. In the same figure, one single Brγ 2012 observation with similar projected baselines is overplotted. The full 2012 dataset for Brγ and He i emission lines is presented in Fig. 3 (online).

The visibility and phase variations shown in Fig. 2 are typical for a rotating disk. As discussed in detail in Meilland et al. (2012), the morphology of visibility variations for a purely rotating disk depends on the baseline orientation. Considering a non-fully resolving baseline, the variation is V-shaped for a baseline aligned with the disk's minor axis (baseline 1 of Fig 2 ) and W-shaped for a baseline aligned with the major axis (baseline Fig. 2. 2010 (in red) and 2012 (in blue) Brγ VLTI/AMBER data for a similar projected triplet of baselines.

of Fig 2)

. Baselines with intermediate orientation exhibit intermediate morphologies (see, for example, the weak W-shaped variation for the third baseline shown in Fig. 2).

The S-shaped phase variations also strongly depend on the baseline orientation. In the case of a purely rotating disk, the amplitude is maximum for a baseline aligned with the major axis and zero for a perpendicular baseline. Thus, from the simple analysis of data shown in Fig. 2 it seems that the disk orientation did not significantly change between 2010 and 2012 with the major axis being oriented roughly in the direction North-South.

However, Fig. 2 also clearly shows important changes in the visibilities, phases, and line profile between 2010 and 2012. Some of them favor the hypothesis of a significant growth of the equatorial disk between the two epochs:

• The visibility drop is clearly deeper in 2012.

• The 2012 S-shaped phase variations are larger.

• The 2012 line profile equivalent width (EW) is larger.

• the double-peak separation is narrower in 2012.

Nevertheless, other major changes seen in Fig. 2 are not related to a growth of the disk:

• The V/R ratio went from 1.2 in 2010 to 0.9 in 2012.

• The W-shaped-visibilities peak ratio is reversed.

• The 2012 S-shaped phase variations are asymmetric.

Such changes of asymmetries in the spectrum, visibilities, and phases across the Brγ line, are clear evidence of changes in the circumstellar matter distribution around the star. In the case of a rotating disk, these asymmetries usually stems from the inhomogeneity of the disk. The origin and nature of the disk's inhomogeneities and their evolution around the periastron are discussed in Sect. 4.3.

We note that all this evidence of the evolution of the disk around δ Scorpii between 2010 and 2012 are seen not only in the data plotted in Fig. 2, but also on all the Brγ and He i dataset (see Fig. 3 online).

Modeling of the 2012 data using the kinematic model

To model the wavelength dependence of the visibility, differential phase, and closure phase in the observed emission lines we used the simple kinematic model developed for fast fitting of an expanding and/or rotating thin equatorial disk already described in paper I. The model is presented in more detail in Delaa et al. (2011) and Meilland et al. (2012). We used the best-fit model parameters of the 2010 data from Paper I as a starting model. As expected, we obtained a very unsatisfactory fit with a χ 2 r of 24.7. We then assumed that the stellar and kinematics parameters should not change between 2010 and 2012, and that only the Brγ and He i disk extension and line equivalent width (EW) should have evolved. In this case we found a best-fit model for FWHM Brγ = 8.0 ± 0.5 R ⋆ (4.2 ± 0.3 mas), EW Brγ = 7.2 ± 0.3 Å, FWHM HeI = 5.2 ± 0.3 R ⋆ (2.8 ± 0.2 mas), and EW HeI = 9.1 ± 0.3 Å. However, for this set of values the χ 2 r is still high, i.e., 17.4. We then performed a fit after releasing constraints on the stellar inclination angle (i), disk major axis orientation (PA), and disk kinematics, i.e., rotational velocity at the photosphere (V rot ) and exponent of the rotational law in the disk (β). Thus, the only fixed parameters for this set of models are the distance taken from van Leeuwen (2007) and the stellar radius taken from Paper I. The best-fit model parameters are presented in Table 2 and the model visibility, phases, and line profiles are overplotted on Fig. 3 (online only).

The χ 2 r for this best-fit model, however, is still high, i.e., 14.6. By separating the fit on each line we found that the He i data have a significantly better fit (i.e., χ 2 r = 8.6) than the Brγ data (χ 2 r = 22.7). For both lines, the largest deviations between the observations and the data comes from the asymmetrical shape of the line profile as well as some visibility and phase variations. This effect will be discussed in the next section. We also note a second source of deviation. Our best-fit model clearly underestimates the visibility for baselines close to the disk's minor axis. However, using a larger disk extension, we would overestimate the visibility for the other baselines and the over-resolution of the disk would strongly modify the shape of the phase variations. Such underestimation of the minor-axis visibility was already detected for α Col, one of the Be stars observed with VLTI/AMBER (Meilland et al. 2012). However, its origins remain unknown.

The disk asymmetry

Okazaki (1997) proposed that the cyclic variation of the V/R peak ratio with a period of several years originates in the precession of one-armed density waves in the disk. This assumption was compatible with spectro-interferometric observations of γ Cas (Berio et al. 1999), κ CMa (Meilland et al. 2007), andζ Tau (Caricofi et al. 2009). This could also explain the change of the peaks ratio of the W-shaped visibility seen in Fig. 2. However, according to observations made by Rivinus et al. (2012), the V/R changes for δ Scorpii occurred during a very short timescale of a few weeks around the periastron. Thus, it seems that the hypothesis of the one-arm oscillation cannot explain the current disk properties.

We also note that, using our simple axisymmetric model, we obtained a better fit on the He i data than on the Brγ data. Considering that the Brγ emission line is formed further from the primary, it is supposed to be more perturbed by the companion tidal forces than the He i line. Consequently, the worse fit of the Brγ data might be a clue of a more important tidal perturbation. One can also argue that the broad-band measurements taken by Che et al. (2012) could not detect any tidal effect as the H-band continuum stems from an even smaller part of the disk than the He i line. Nevertheless, in the frame of an asymmetry induced by tidal effects in June-July 2011, it is also difficult to understand why the V/R remained almost constant between July 2011 and our observations in June 2012. One should expect the induced inhomogeneity to be either in co-rotation or precession, causing cyclical V/R variations with timescales ranging from weeks to years that are not seen in the available spectroscopic data.

Evolution of the disk extension

Our VLTI/AMBER data clearly gives evidence that the Brγ emission FWHM increased by about 40% (i.e., 2.3 ± 1.3 R ⋆ ) between 2010 and 2012. Assuming that this was directly cause by a physical growth of the circumstellar disk we derived a mean expansion velocity of 0.20±0.11 km s -1 . This value is similar to the one previously determined for δ Scorpii by Miroshnichenko et al. (2003) for the 2001-2003 period, i.e., 0.4 km s -1 , and the one derived in Paper I for the 2000-2005 and 2005-2007 periods, 0.24 km s -1 and 0.19 km s -1 , respectively. It is also compatible with the one determined by Kanaan et al. (2008) for the transient disk surrounding the binary Be star Achernar, i.e., 0.27 ± 0.08 km s -1 .

The He i line, however, did not exhibit the same variation as the Brγ line. Its growth, of about 8%, remains within the uncertainties of the derived extension. One possible explanation for this difference is that He i has a higher excitation energy. Consequently, it is produced by denser and hotter material than Brγ. Thus, even if the disk is increasing in size as a function of time the He i emitting region should remain compact, whereas the Brγ may be produced in the newly formed outer part of the disk.

In a recent paper, Stee et al. ( 2012) measured the disk extension in these two lines for the relatively steady disk of the Be star HD110432. The values derived for their measurement using the same methods showed that the He i emission is about 25% smaller than the Brγ emission. This value is between the one we found for δ Scorpii for our 2010 (18%) and 2012 (37%) models.

Conclusion

Using the unique high spectral and spatial resolutions of the VLTI/AMBER instrument, we studied the evolution of δ Scorpii's circumstellar environment one year before and one year after the 2011 periastron passage. At both epochs, the data clearly shows that the ejected matter is located in a thin equatorial disk in Keplerian rotation as it is for other observed Be stars. Since 2000, the disk seems to have grown with a quasi-constant velocity of about 0.2 km s -1 and was still growing during the 2010-2012 period. Such speed is compatible with the viscous excretion disk model proposed by [START_REF] Lafrasse | JMMC Stellar Diameters Catalogue[END_REF].

We found evidence that the disk was asymmetric at both epochs of our observations, but with a different morphology in 2010 and 2012. According to spectroscopic observations made by Rivinius et al. (2012), the major modifications occurred near the periastron. Thus, the observed asymmetry cannot originate from a one-arm oscillation as proposed by Okazaki (1997), but must stem from tidal effects induced by δ Scorpii's companion. However, according to available spectroscopic data and our 2012 observations, the disk asymmetry seems to be "frozen" since the periastron. Finally, long-term spectroscopic and interferometric follow-up are needed to understand both the nature of the inhomogeneity after the periastron and its link with the tidal effects of the companion.

Brγ differential visibilities differential phases, and closure phases He I differential visibilities differential phases, and closure phases Fig. 3. The full dataset for δ Sorpii 2012 VLTI/AMBER observations. The red lines represents the VLTI/AMBER data, and the green ones, the best-fit symetric kinematic model.

Fig. 1 .

 1 Fig. 1. δ Scorpii (u,v) 49lan coverage obtained in 2012 with VLTI/AMBER observations centered on Brγ (left) and 2.06 µm He i (right) emission lines.

Table 1 .

 1 VLTI/AMBER log for the 2012 δ Scorpii observations.

	Observation Date Time	Triplet	Projected Baselines L. (m) P.A.( o )
	Observations centered at 2.17 µm (Brγ)
	2012-06-22 23:46 D0-H0-G1 53/ 71/ 63 2012-06-23 00:32 D0-H0-G1 57/ 71/ 67 2012-06-23 01:21 D0-I1-H0 79/ 40/ 62 2012-06-23 02:05 D0-I1-G1 82/ 46/ 71	48/ 170/ 124 57/ 175/ 126 94/ -36/ 63 98/-143/ 132
	Observations centered at 2.06 µm (He I ) 2012-06-23 23:16 D0-H0-G1 50/ 71/ 60 43/ 167/ 123 2012-06-24 00:07 D0-I1-H0 70/ 39/ 56 87/ -40/ 53 2012-06-24 00:52 H0-I1-G1 40/ 45/ 71 142/-152/ 177 2012-06-24 01:35 D0-I1-G1 80/ 46/ 71 95/-146/ 130 2012-06-24 02:18 D0-I1-G1 82/ 47/ 71 99/-141/ 134

Table 2 .

 2 The 2010 and 2012 best-fit kinematic models.

	Parameter	(unit)	2010	2012
	R ⋆ (fixed) d (fixed) i PA	Global geometric parameters (R ⊙ ) 8.5 (pc) 150 (deg) 30 ± 8 (deg) -12 ± 7	8.5 150 26 ± 5 -14 ± 4
		Global kinematic parameters
	V rot β	(km.s -1 )	500 ± 50 0.5 ± 0.1	500 ± 50 0.5 ± 0.05
	FWHM Brγ	Brγ disk geometry (R ⋆ ) 5.5 ± 1	7.8 ± 0.3
	EW Brγ	(Å)	6.5 ± 0.5	6.2 ± 0.2
	FWHM HeI	He I disk geometry (R ⋆ ) 4.5 ± 0.5	4.9 ± 0.2
	EW HeI χ 2 r	(Å)	8.5 ± 0.5 6.2	8.8 ± 0.3 14.6

http://basebe.obspm.fr

available at https://www.oca.eu/spip.php?article562

Acknowledgements. The Programme National de Physique Stellaire (PNPS) and the Institut National en Sciences de l'Univers (INSU) are acknowledged for their financial supports. A. Meilland acknowledge financial support from the University of Western Ontario. The authors would like to thank ESO staff for their help to make the observations a success and more particularly J-B Lebouquin and F. Rantakyro.