Sébastien Peirani

Intae Jung

Joseph Silk

Christophe Pichon

Evolution of the baryon fraction in the Local Group: accretion versus feedback at low and high z

Keywords: Galaxies, Local Group -Galaxies, haloes -Dark matter -Methods, N-body simulations

Using hydrodynamical zoom simulations in the standard ΛCDM cosmology, we investigate the evolution of the distribution of baryons (gas and stars) in a local group-type universe. First, with standard star formation and supernova feedback prescriptions, we find that the mean baryonic fraction value estimated at the virial radius of the two main central objects (i.e. the Milky Way and Andromeda) is decreasing over time, and is 10-15% lower than the universal value, 0.166, at z = 0. This decrease is mainly due to the fact that the amount of accretion of dissipative gas onto the halo, especially at low redshift, is in general much lower than that of the dissipationless dark matter. Indeed, a significant part of the baryons does not collapse onto the haloes and remains in their outskirts, mainly in the form of warm-hot intergalactic medium (WHIM). Moreover, during the formation of each object, some dark matter and baryons are also expelled through merger events via tidal disruption. In contrast to baryons, expelled dark matter can be more efficiently re-accreted onto the halo, enhancing both the reduction of f b inside R v , and the increase of the mass of WHIM outside R v . Varying the efficiency of supernovae feedback at low redshift does not seem to significantly affect these trends. Alternatively, when a significant fraction of the initial gas in the main objects is released at high redshifts by more powerful sources of feedback, such as AGN from intermediate mass black holes in lower mass galaxies, the baryonic fraction at the virial radius can have a lower value (f b ∼ 0.12) at low redshift. Hence physical mechanisms able to drive the gas out of the virial radius at high redshifts will have a stronger impact on the deficit of baryons in the mass budget of Milky Way type-galaxies at present times than those that expel the gas in the longer, late phases of galaxy formation.

INTRODUCTION

The distribution of galaxies in the Universe on large scales, as revealed by sky surveys such as the CFA Redshift Survey (Geller & Huchra 1989), the 2DF Galaxy Redshift Survey (Colless et al. 2001) or the Sloan Digital Sky Survey (Gott et al. 2005) is structured into a ⋆ E-mail: peirani@iap.fr filamentary web in which the intersections correspond to the most massive objects (i.e. galaxy clusters, Bond et al. 1996). This specific distribution is now commonly and easily reproduced by cosmological hydrodynamical N-body simulations of the ΛCDM cosmology which also indicate that the distribution of dark matter, the dominant mass component of the Universe (Komatsu et al. 2011), follows the same trend. Indeed, since in the traditional picture of galaxy formation, galaxies are supposed to form when baryonic gas falls into the gravitational potential of their host dark matter halo, the distribution of dark matter is therefore expected to faithfully trace that of the baryons. However on galactic scales, observations tend to suggest that the spatial distributions of dark matter and baryons (especially in the form of gas) may display some substantial differences. In particular, it has been shown that galaxies are missing most of their baryons, -most galaxies are severely baryons-depleted relative to the cosmological fraction (see for instance Bell et al. 2003;Hoekstra et al. 2005;McGaugh 2010).

This so-called "missing baryons problem" (see Bregman 2007 for a complete review), if real, calls for two alternative scenarios. Either a significant part of the gas never collapsed into the gravitational potential wells of protogalaxies in the first place, or some of the gas has been expelled by galaxy formation feedback processes such as supernova winds. Hence solving the missing baryon problem may prove to be central in order to constrain galaxy formation models.

Past studies based on hydrodynamical simulations indicate that most of the "missing" baryons might lie in a gaseous phase (the so-called warm-hot intergalactic medium, hereafter WHIM) in the temperature range 10 5 -10 7 K and at moderate over-density (Cen & Ostriker 1999;Davé et al. 2001;Cen & Ostriker 2006). According to these studies, the WHIM is primarily shockheated during the formation of large-scale structures, while feedback mechanisms associated with star formation should also have an additional impact on this phase during galaxy formation. For instance, for low mass systems, the cosmological UV-background is supposed to reduce both the star formation and baryon content (Thoul & Weinberg 1996;Bullock, Kravtsov & Weinberg 2000;Gnedin 2000;Somerville 2002;Ricotti, Gnedin & Shull 2002;Benson et al. 2002;Read, Pontzen & Viel 2006;Hoeft et al. 2006;Okamoto, Gao & Theuns 2008;Peirani 2010;Nickerson et al. 2011, Hambrick et al. 2011a). Supernova feedback is also expected to expel some gas, especially for low mass systems, as suggested by various numerical investigations (Scannapieco et al. 2008(Scannapieco et al. , 2009;;Faucher-Giguère, Kereš & Ma 2011 and references therein). Other mechanisms have been proposed, such as turbulence in the baryonic intergalactic medium (Zhu, Feng & Fang 2011) or pre-heating by pre-virialisation (Mo et al. 2005), although Crain et al. (2007) claimed that this process is unable to prevent the collapse of gas (instead, according to these authors, non-gravitational feedback is required in order to reduce the efficiency of gas cooling and star formation in dwarf galaxies). Another less investigated potential mechanism, which may also affect the evolution of the distribution of baryons on small scales, is the effect of accretion and merger events (and in particular their associated tidal disruption) that are expected to be frequent in the framework of the hierarchical model. For instance, in an early paper, Navarro & White (1993) have shown that during a merger involving dark matter and adiabatic gas, there is a transfer of energy between these two components that leads to a situation where the gas is less tightly bound. More recently, Sinha & Holley-Bockelmann (2010, 2011) found that a few per-cent of gas can be driven into the intergalactic medium (IGM) by galaxy mergers, using either merger trees calculated from the Press-Schechter formalism or idealized simulations of mergers of galaxy clusters. Using hydrodynamical simulations, Nickerson et al. (2011) have also concluded that tidal forces may cause significant mass loss from satellites of all masses.

In the present work, we make use of cosmological "zoom" simulations with an extended treatment of the physics of baryons to study the formation of Milky Waylike galaxies. Our aim is to characterise the relative role of supernova feedback to accretion and mergers in the evolution of the distribution of baryons for objects of such masses. We will also test scenarios in which a significant fraction of gas in progenitors is expelled at high redshift by more powerful sources of feedback, such as AGN associated with massive black holes. This will allow us to quantify two distinct processes which may allow us to address the so-called missing baryon problem, should it persist.

This paper is organised as follows. In section 2, we summarize the numerical modelling, and section 3 presents our main results on the evolution of the distribution of baryons in our simulated local group universes. We give our main conclusions in the last section.

NUMERICAL MODELLING

Initial conditions and simulation

The numerical methodology used in the present paper is described in detail in Peirani (2010) to which we refer the reader for more information. For the sake of clarity, we summarize the main steps below.

We analyse three cosmological zoom simulations for a ΛCDM universe using WMAP5 parameters (Komatsu et al. 2009), namely ΩM = 0.274, ΩΛ = 0.726, Ω b = 0.0456, H0 = 70.5 km/s/Mpc, n = 0.96 and σ8 = 0.812. Each simulation was performed in a periodic box of side 100 h -1 Mpc with 2 × 2048 3 effective dark matter and gas particles in the highest resolution region (a sphere of 7 h -1 Mpc of radius). In this region of interest, the mass resolution are mDM ≈ 7.4 × 10 6 h -1 M⊙ and mgas = mstars ≈ 1.5 × 10 6 h -1 M⊙. The Plummerequivalent force softening adopted for the high mass resolution particles were 1 and 0.5 h -1 kpc for dark matter and gas particles respectively and were kept constant in comoving units.

Initial conditions have been generated from the MPgrafic code (Prunet et al. 2008) and the simulations were performed using GADGET2 (Springel 2005), with added prescriptions for metal-dependent cooling, star formation (in this work we remind that the star formation efficiency is c * =0.02), feedback from Type Ia and II supernovae (SN), UV background (starting at z = 8.5) and metal enrichment. The three simulations have common initial conditions whose phases are consistent with the local group but essentially differ in the quantity of energy released by SN derived from star particles i (Ei). As mentioned in Peirani (2010), we consider that a fraction γ of this energy is deposited in the j th neighbour gas particle by applying a radial kick to its velocity with a magnitude ∆vj = (2wjγEi/mj ), where wj is the weighting based on the smoothing kernel and mj is the mass of gas particle j. The first simulation Sim1a uses the standard value of γ = 0.1 while in the second one, Sim1b, we have considered a higher efficiency γ = 1.0 in order to investigate how our results would be affected. In the third simulation (Sim2), our aim is to study the effects of earlier high energy ejection to the ISM induced either by intermediate mass black holes or other high energy processes such as hypernovae events. Intermediate mass black holes are the likely missing link between Population III and the supermassive black holes in quasars and plausibly as part of the hierarchy of structure formation. Hypernovae are likely to be more prevalent than normal supernovae in the earliest phases of star formation: for example, the Population III IMF is most likely top-heavy, and this may also be the norm for the precursors of spheroids, as evidenced for example by the radial distribution of SNe in disturbed (including recently merged) galaxies (Habergham, Anderson & James 2010) and the frequency of low mass x-ray binaries in ultracompact dwarf galaxies (Dabringhausen et al. 2012). A higher frequency of hypernovae could plausibly provide an order-of-magnitude higher feedback efficiency than supernovae. The case of AGN is intriguing, as even higher feedback is required in deep potential wells, where supernovae are relatively inefficient, in order to account for the correlation between black hole mass and spheroid velocity dispersion (Silk & Nusser 2010;Debuhr et al. 2012). Moreover, the recent detection of a massive gas outflow in a distant quasar (z = 6.4) strongly suggests that a strong quasar activity is already in place at very early times (Maiolino et al. 2012). For this purpose, a simple modelling was used in which a much higher efficiency (γ = 50) was considered during a very short (∆t ∼ 45 Myr) at earlier times (z ∼ 8.0) and γ = 0.1 otherwise. We also justify the choice of this high efficiency by simple energetics and momentum comparisons between supernovae and AGN. For instance, for a 10 7 M⊙ black hole along with its 10 10 M⊙ in stars, the number of type II SNe produced in star formation phase (∼ 150 M⊙ per SN for a Chabrier IMF) is ∼ 7.10 7 . Thus, the total SNe energy injected is ∼ 10 59 ergs over a dynamical time, say t dyn ∼ 10 8 yr, or 10 43.5 ergs/s. Since the Eddington luminosity is L Eddington ∼ 10 45 ergs/s energetics favour AGN. Furthermore, for supernovae momentum conservation starts at shell velocity ∼ 400 km/s. Since the ejection velocity is ∼ 2.10 4 km/s, the momentum injected is ∼ 2% of initial momentum namely ∼ 10 48 gm cm/s or ∼ 10 32.5 dynes. For AGN, momentum injected is L Eddington /c or ∼ 10 34.5 dynes. Thus AGN injects 100 times more momentum flux than supernovae. Finally, there is another factor of 10-100 that favours AGN, the so-called mechanical advantage factor due to ram pressure on the expanding bubble as suggested by recent hydrodynamical simulations (Wagner, Bicknell & Umemura). Therefore the value of γ = 50 seems to be reasonable in order to account for the accumulation of the different potential sources of feedback at high redshifts. The feedback parameters are summarized in Table 1.

Physical properties of main galaxies

In each simulation, we analyse a pair of galaxies with physical characteristics similar to the Milky Way-Andromeda pair (MW-M31) and with similar galaxy environments up to 7 h -1 Mpc. Some relevant physical properties at the virial radius Rv at z = 0 of these two objects from each simulation are summarized in Table 2. Note that we keep the same definition of the virial radius RV and virial mass MV used in Peirani (2010) (see paragraph 2.3), namely RV is the radius where the enclosed mean density MV /(4πR 3 V /3) is ∆c times the critical density, and ρc(z) = 3H(z) 2 /8πG, where H(z) = H0 ΩM(1 + z) 3 + ΩΛ. In the cosmology adopted in the present study, ∆c = 97.6 at z = 0. Furthermore, the baryonic fraction is defined as follow:

f b ≡ m baryons m baryons + mDM , (1)
where m baryons and mDM refer to the masses of baryons (gas+stars) and dark matter respectively in the same specific region of the universe. When not specified in the text, we do not distinguish the cold and hot phases of the gas in the estimate of f b . According to the adopted cosmology, the universal f b value is:

f b ≡ Ω b ΩM ≈ 0.166. ( 2
)
The time evolution of the virial, dark matter (DM), gas and stellar masses of each object at the virial radius and for the three simulations is shown in Figure 1. First, we notice that there is no particular difference between the evolution of dark matter halo masses (and hence viral masses) between Sim1a and Sim1b. However, the baryonic compositions are quite different. Indeed, stellar masses are significantly reduced in Sim1b due to higher amounts of energy released by the SN. The mass of the gas component is therefore higher in Sim1b, as expected. For Sim2, the stellar masses are even more reduced. This is mainly due to the fact that a significant fraction of the gas has been expelled at high redshifts.

Indeed, while the DM masses are quite similar between the 3 runs, the virial masses for objects in Sim2 tend to be lower at high redshift, which indicate that a fraction of gas expelled at high z does not recollapse at some later stage. Interestingly, Halo2 undergoes a minor merger event at z ≃ 0.1, which corresponds for instance to a sudden variation in the evolution of its stellar mass. This same accretion event takes place at z ≃ 0.2 in Sim2. Moreover, we can also identify clearly in Sim2 another minor merger event at z ≃ 0.05 that has not occurred (yet?) in Sim1a and Sim1b. This explains why the dark matter mass (and therefore the virial mass) of Halo2 is higher in Sim2 than Sim1a or Sim1b. This also suggests that small perturbations in the properties of haloes at earlier times can have some significant impact on the whole accretion history (see for instance, Thiebaut et al. 2008). It is also worth mentioning that although the two main central objects have the same virial mass at z = 0, they have undergone different mass accretion histories. Indeed, while Halo2 grows regularly through smooth mass accretion, the evolution of the mass of Halo1 is more sporadic thanks to several major merger events. It is therefore instructive to see to which extent these two different mass evolutions could affect the evolution of the distribution of baryons in the vicinity of these two objects.

It is also worth studying how the stellar masses found in the different simulations are related to the mass of their host dark matter haloes. From Table 2, the mass of dark matter haloes varies from 1.02×10 12 M⊙ to 1.09 × 10 12 M⊙ while according to observational and numerical analysis, the expected stellar mass range inside those haloes is M * ≃ [2.5 -3.5] × 10 10 M⊙ (Guo et al. 2010;Moster et al. 2010;Brook et al. 2012b). In order to make a suitable comparison, instead of considering the stellar mass inside the virial radius Rv, we define the galaxy radius as R gal = 0.1 Rv and all stars inside this radius contribute to the galaxy stellar mass M * . By doing this, we only select stars which belong to the central main galaxy while excluding stars of accreted satellites inside Rv. In Figure 1, the evolution of each stellar mass inside R gal is shown in dashed line. One can see that at z = 0, galaxy stellar mass ranges are

M * ,Sim1a ≃ [1.1 -1.2] × 10 11 M⊙, M * ,Sim1b ≃ [9.9 -1.1] × 10 11 M⊙ and M * ,Sim2 ≃ [7.0 -7.4] × 10 11 M⊙ for Sim1a, Sim1b
and Sim2 respectively. Thus, each simulation tends to produce too many stars in the central region of haloes and the trend is more severe in Sim1a and Sim1b relative to Sim2. This problem was already pointed out from other hydrodynamical simulations which indicate that low feedbacks may result in more than an order-ofmagnitude too many stars (see for instance Piontek & Steinmetz 2011;Brook et al. 2012b). As a first conclusion, stronger and earlier feedbacks seem to be required in order to improve the stellar mass-halo mass relation (also recently suggested by Stinson et al. 2012) and, in the present study, one probably also needs the action of additional sources of feedback in the late phase of galaxy formation such as AGN to get even closer to the expected galaxy stellar masses. Moreover the fact that our model leads to too high central stellar masses suggests that the effect of all different sources of feedback has been underestimated. Thus in the following, baryonic fraction values found in our simulations may be slightly overestimated.

Finally, Figure 2 shows the evolution of the star formation rate (SFR) for each object and simulation. At z = 0, SFR values are around [1-5] M⊙/yr which is consistent with observational values derived by Robitaille & Whitney (2010) using Spitzer data. Note that these cosmic evolutions are similar for the two haloes: at high redshifts (z ≥ 2) the lower the SN efficiency the higher the stellar mass. And the higher the mass of formed stars at high redshift, the higher the subsequent SN feedback. This explains why, at low redshifts, a higher fraction of stars is produced in Sim2 relative to Sim1a and Sim1b.

THE BARYON FRACTION IN THE LG

3.1 Spatial distribution in the simulated LG Fig. 3 shows the projected distribution of baryons in our LG type universes derived from Sim1b and Sim2 and at four specific epochs (i.e. z = 5, 2, 0.5 and 0). We clearly see that the evolution of the distributions of baryons and dark matter do not follow the same trend. Indeed, the regions of the universe in red and dark blue correspond to regions where the baryonic fraction is lower than the universal value. And those specific regions are mainly located around the forming protogalaxies or galaxies at each redshift, while in the very inner part of dark matter haloes, f b is higher. This is an expected result: due to cooling, the gas collapses to the center of haloes where stars can be formed. However, note that the size of these "red" cavities increases over time which suggests that the fraction of gas that has collapsed to the center of the halo is not immediately replaced by some fresh gas from its vicinity.

Note also, at high and low redshift, the existence of relative high baryon content regions which are located either in the filaments or beyond the virial radii. Such anisotropic distribution seems to be more pronounced in Sim2 relative to Sim1b. If the high baryon fractions in filaments can be understood by the dissipative nature of gas, allowing it to cool to the dense filaments, it is crucial to characterise the mechanisms that drive high baryon content regions beyond the virial radii.

Cosmic evolution of f b

The effects of feedback at low redshift

Let us first only consider the role of supernovae feedback in the evolution of the mass budget of Milky Way-type galaxies. We therefore focus on both Sim1a and Sim1b in this sub-section.

The cosmic evolution of the baryonic fraction f b estimated at the virial radius for the two main haloes is shown in Figure 4. In each of these two simulations, they follow the same trend. At high redshift, f b is close to the universal value and sometimes is slightly higher. Cold flows provides gas to form stars. Note that the virial radius is estimated according to the local density, which is dominated by the newly formed stars. Since the resolution in these simulations is limited, especially for low mass systems, the virial radius may be underestimated, and subsequently this baryonic fraction may be in fact overestimated at high redshifts.

From z ∼ 3, f b is decreasing with cosmic time until it reaches the values of ∼ 0.15 and ∼ 0.14 at z = 0 for Halo1 and Halo2 respectively. These latter values are quite close to the universal value which is not surprising since similar trends have been already found from other hydrodynamical simulations using weak or no supernovae feedback, for objects with a mass of 10 12 h -1 M⊙ at z = 0 (Crain et al. 2007;Faucher-Giguère, Kereš & Ma 2011). More interestingly, no particular differences are seen in the evolution of f b between Sim1a and Sim1b. This suggests that higher SN feedback can reduce the star formation rate (and therefore the final stellar mass) but is rather inefficient in expelling the gas outside the virial radius at high and low redshifts for massive galaxies. This results are in agreement with Stinson et al. (2011) who show that increased feedback in haloes of such mass scale affects the star formation more than baryon content in the circumgalactic medium (CGM). Moreover, they show that higher feedback models can account for sufficient OVI and HI gas in the CGM compatible with the observed distributions. Therefore, in the present scenario, the relatively small differences in the baryon fraction between Sim1a and Sim1b could make supernovae feedback prescriptions difficult to distinguish and one probably needs to focus on stellar mass and/or properties of the intergalactic medium in order to have better constraints (see also Davé, Oppenheimer & Sivanandam 2008;Shen, Wadsley & Stinson 2010;Stinson et al. 2011Stinson et al. ,2012;;Scannapieco et al. 2012;Brook et al. 2012b). In particular, it seems that higher feedback of Sim1b compared to Sim1a may be necessary in reproducing values closer to the correct stellar masses and properties of intergalactic medium as we will see in section 3.3.

We have also studied the distribution of the integral f b (< R) and differential f b (R) as respect to the radius R at three different redshifts z = 4.7, 0.8 and 0 in Figure 5. Only results from Halo2 are shown since similar trends are obtained for Halo1. At the center of each halo, f b is high due to the presence of the galaxy. At larger radii, R 3Rv, f b (< R) tends towards the universal value as expected. Interestingly, for high redshifts (z > 4.7), f b (< R) converges to the universal value from above, while this is not the case at lower redshifts, which suggest, as mentioned above, that the gas infall is not recurrently replaced by some fresh gas from the outskirt of the haloes. Also, an excess of baryons at (1.5 -2.5)Rv can be clearly seen in the variations of the differential f b (R) for different cosmic times, which correspond to regions in yellow in Fig. 3 around haloes. Similar plots were recently derived from observations using Chandra and Suzaku facilities for an isolated elliptical galaxy with a ∼ Milky way mass (Humphrey et al. 2011) and a fossil group (Humphrey et al. 2012).

The effects of feedback at high redshift

In Sim2 however, due to earlier and important eviction of gas from other potential sources of feedback (AGN, hypernovae,...) f b reaches its lowest value f b ∼ 0.05 at z ∼ 6, then increases until z ∼ 1.5 and becomes nearly constant (f b ∼ 0.12) until the present time. This strongly suggests that sources of feedback acting at high redshift, even for a very short period, can have a stronger impact on the final mass budget of massive galaxies than those acting in the later and longer phases of galaxy formation. One can also notice from Figure 2 that the gas outflow taking place at high redshift lead to the quenching of the star formation during a few million years (5 ≤ z ≤ 8). This is in nice agreement with conclusions of Valiante, Schneider & Maiolino (2012) using semi-analytical models of the formation and evolution of high-redshift quasars. Indeed, according to their model, the observed outflow of a distant quasar at z = 6.4 is dominated by quasar feedback which "has considerably depleted the gas content of the host galaxy, leading to a down-turn in the star formation rate at z < 7 -8". Nevertheless, future observational proofs of star formation quenching at these redshifts are required in order to confirm such statements.

The role of accretion of matter

The time evolution of the baryonic fraction inside the virial radius depends on the competition between the evolution of the mass that arrives within Rv, and the mass that leaves Rv. It raises the following question: are the low values of f b found at z = 0 due to the fact that more baryons have been expelled out of the virial radius or because more dark matter has been accreted inside Rv? To answer this question, instead of studying the evolution of the mass that enter or leave the virial radius at each subsequent step, we focused on the location at a given redshift z of all particles that have been accreted up to this specific redshift. To do this, we have first identified all the particles that have been accreted at each snapshot of the simulation from z ∼ 8 to "z". Then, we have computed the fraction of those particles that are inside or outside Rv at the considered redshift. By doing this, the interpretation of our findings will be easier since each accreted particle contributes only once in the final result. Indeed, certain particles can be accreted, ejected and accreted again during the whole evolution of these haloes.

Figure 6 shows the evolution of the total accreted DM and baryon mass from z ∼ 8 to a given redshift, z, for both haloes, and for both Sim1b and Sim2 (results derived from Sim1a are not shown since they are similar to those from Sim1b). We also plot the fraction of these accreted particles that are located inside or outside Rv at this given redshift z. First, we look at the evolution relative to Halo1. One can see that a higher fraction of dark matter mass has been accreted from z = 8 to z = 0 relative to the baryon component (solid lines). Moreover, a similar fraction of dark matter and baryons are located outside the virial radius at z = 0. That means that the fraction of baryons released during the formation of Halo1 is not higher than that of dark matter, which demonstrate that the low value of f b at z = 0 follows from the fact that more dark matter particles has been accreted over cosmic time. This can be understood given the difference in nature between the dissipative gas and the non-dissipative dark matter. In particular, the gas component can be shockheated at the virial radius which considerably affects and slows down its further accretion on to the halo as already shown in hydrodynamical simulations using no feedback (see for instance Oñorbe et al. 2007). In contrast, especially after an important episode of accretion or merger event, the dark matter component can "damp out" its orbital motion until it completely relaxes. For Halo2 the same trends are observed but the mass of baryons outside the virial radius turns out to be higher that that of dark matter. However this does not necessarily mean that a higher fraction of baryons has been expelled through physical mechanisms such as feedback or merger events. Indeed, if we represent the fraction of baryons and DM released from z ∼ 0.6 to z = 0 (see Figure 6. The evolutions of the total dark matter mass (solid black line) and baryon mass (solid grey line) that was accreted from z ∼ 8 to a given redshift z (see the main text for more details) for Halo1 (first column) and Halo2 (column 2), derived from Sim1b (first line) and Sim2 (second line). In each panel, the corresponding fractions of mass of DM (black color) and baryons (grey color) that are inside Rv (dashed lines) or outside Rv (dotted line) is represented. The dotteddashed black and grey lines represent respectively the fraction of DM and baryons outside Rv , corresponding to particles accreted from z ∼ 0.6 to a given redshift z . Note that the baryon mass evolutions have been multiplied by (Ω M -Ω b )/Ω b to facilitate a direct comparison with DM mass evolutions.

Figure 6), one can see that same amounts of mass have been expelled at z = 0 by the two components. It implies in fact that a significant fraction of dark matter which has been expelled between z = 8 and z ∼ 0.6 has been re-accreted between z ∼ 0.6 and z = 0. This is not the case for baryons which are generally heated through the expelling process and therefore remain outside the virial radius afterwards. In order to illustrate this, Figure 7 presents the evolution of the distribution of baryons and dark matter for Halo1. From z = 2.6, the main progenitor undergoes two mergers with mass ratios of 1:1.8 and 1:5. Here, we have selected all particles inside the virial radius of each of the 3 objects. Note that the orbits of the two mergers are quite radial. After the two merger events, the spatial distributions of baryons and dark matter are quite different. Most of the dark matter particles are located inside the virial radius, though a small part has been expelled in the same orbital direction. The baryons tends to either collapse toward the center or are expelled all around. This phenomenon is observed at z = 0 where a non-negligible part of baryons have been expelled and remain outside the virial radius. Figure 7 also shows the evolution of the mass of DM and baryons located outside the virial radius. In the case of the dark matter component, there are 3 distinct phases. First, from the beginning of the experiment until t = 3 Gyr, both satellites are accreted inside Rv. Then in the second phase (from t = 3 to t ∼ 4.5 Gyr) some of the dark matter is expelled via tidal disruption. Finally from t ∼ 4.5 Gyr some of the dark matter particles are re-accreted. Note that evolution from t ∼ 10 Gyr is probably due to the effect of another merger event. However, the baryon component follows the two first phases but not the third one: once baryons are expelled, most of them remain outside the virial radius. In this specific example, we found that 2.1% and 5.8% of dark matter and baryon masses respectively are released at z = 0 from the three main objects. This is in good agreement with results obtained from idealized simulations (Sinha & Holley-Bockelmann 2009). It is also interesting to see that due to mergers, unbound gas diffuses far away from the galaxy; thus mergers represent a potential mechanism for enriching the IGM with metals. Note that this enrichment must occur early, as in our model, to avoid disruption of the cold and weakly enriched Lyman alpha forest. It is also worth mentioning that during a merger event, gas can be expelled via both tidal disruption and outflows from feedback due to the merger driven starburst. However, since there is no significant difference in the evolution of f b within the virial radius between Sim1a and Sim1b, this suggests that supernovae feedbacks have a lower contribution relative to tidal disruption for Milky Way mass galaxies. As observed in the evolution of the dark matter component, tidal disruption should guarantee escape whereas a starburst driven outflow is much more problematic (see for instance Powell, Slyz and Devriendt 2011).

Finally, it is interesting to point out that Halo1 undergoes major merger events between z ∼ 0.8 and z ∼ 0.5, leading to important mass outflows from z ∼ 0.5 and z = 0. In particular, Figure 6 indicates that more dark matter particles are released than baryons after the merger events. Most of this dark matter is not yet re-accreted, which explains why Halo1 tends to have values of f b higher than Halo2 at z = 0. This also explains why in Figure 3 the high f b value regions around Halo1 seem to be less pronounced in comparison to those around Halo2.

Thus, the evolution of the mean baryonic fraction value f b at the virial radius is essentially governed by the relative efficiency at which the dark matter and baryons are accreted. Indeed, the evolution of the accreted mass of both dark matter and baryons follows two different regimes: a rapid growth at high redshifts (z ≥ 1.5) and a slower one at lower redshifts. These evolutions are directly related to the total mass evolutions in Figure 6 which also presents these two regimes. This behaviour compares well to the evolution of individual dark matter haloes (see for instance Wechsler et al. 2002;Vitvitska et al. 2002, Peirani, Mohayaee & Pacheco 2004, McBride et al. 2009and Tillson, Miller & Devriendt 2011), indicating a fast growth of the mass for z > 1.5 -2.0, followed by a phase where the accretion rate is slower. Hydrodynamical simulations have also suggested that galaxy formation presents two phases: a rapid early phase (z ≥ 1.5) during which stars mainly form within In the lower panels, we show the evolution of the mass of baryons and DM outside the virial radius of the host halo (left panel) and the evolution of the baryonic fraction estimated from these 3 haloes inside and outside the virial radius (right panel). One can see that a non negligible number of DM particles are re-accreted inside Rv (from t ∼ 4.5 to t ∼ 10 Gyr) which is not the case for baryons. As a result, the final values of f b (< Rv ) and f b (> Rv) are respectively lower and higher than initial ones.

the galaxy from infalling cold gas [START_REF] Katz | The IGM/Galaxy Connection[END_REF]Kereš et al. 2005;Ocvirk, Pichon & Teyssier 2008;Dekel et al. 2009) and a second one (z ≤ 1.5) during which stars are accreted (Oser et al. 2010;Hirschmann et al. 2012). This latter phase is dominated by the accretion of diffuse matter and small satellites (Fakhouri & Ma 2010;Genel et al. 2010;Wang et al. 2011;van de Voort et al. 2011). Hence the high f b values at high redshift derived from Sim1a and Sim1b can be explained by the accretion of high gas-to-dark matter ratio from dense region such as filament via cold flows. But at lower redshift, f b is decreasing due to the accretion of low gas-to-dark matter ratio material, in particular from the diffuse region. Indeed, we found that the baryonic fraction of the diffuse accreted matter (namely f b ∼ 0.11 -0.14) is on the average much lower than the universal value: as the haloes become more massive, the temperature at the virial radius increases and the gas is shock-heated and this process slows down its accretion onto the halo. For Sim2, f b reaches its lowest value (f b ∼ 0.05) at z ∼ 5 right after some significant expulsion of gas. After that, there is a short phase where f b increases because the accreted matter has a higher baryonic content (but not necessarily higher than the universal value). Then, in the late phase of evolution of each halo, f b tends to be constant (f b ∼ 0.12) because the average accreted matter tends to have the same baryonic fraction value as the mean value inside virial radii.

Properties of the IGM

The WHIM has proven to be difficult to detect and represents a serious candidate for the missing baryons.

It is therefore instructive to estimate the fraction of the WHIM in our simulated local universes. For this purpose, we have divided the Inter-Galactic Medium (IGM) into four components: (1) stars, (2) cold gas (T < 10 5 K), (3) WHIM (10 5 K < T < 10 7 K) and hot gas (T > 10 7 K). Figure 8 shows the time evolution of three baryon components (stars, cold gas and WHIM) within a sphere centred on the local Universetype haloes and of comoving radius 1.5 Mpc/h. The evolution of the hot gas (T ≥ 10 7 K) is not shown in Figure 8 since this component represents less than 1% of the total baryon mass budget in our selected region, as expected since we are not in a cluster environment. It shows that the WHIM represents ∼ 30% of the baryon budget at z = 0, in good agreement with past work based on hydrodynamical simulations (Cen & Ostriker 1999;Davé et al. 2001;Cen & Ostriker 2006;Rasera & Teyssier 2006;Tornatore et al. 2010) when AGN feedback is not taken into account. We also found that the warm and hot gas becomes the dominant baryonic component at large radii (R > 150 kpc) from each main halo, in good agreement with previous observational analysis of massive galaxies (Humphrey et al. 2011(Humphrey et al. , 2012) ) or with previous numerical simulations (van de Voort et al. 2011; van de Voort & Schaye 2012). However, when the contribution of the WHIM is not taken into account in the estimation of the baryonic fraction at the virial radius, namely when only the cold gas (T ≤ 10 5 K) and stars are considered, f b values obtained at z = 0 only decrease by ∼ 15%. This follows from the fact that stellar masses obtained in our experiments have a large contribution in the final f b values and thus we probably need the additional effect of AGN feedback to regulate the SFR and also increase the mass of WHIM inside the virial radius. Nevertheless, the distribution of baryonic fraction seen in Figure 3 is significantly affected when the WHIM is not taken into account. Indeed, Figure 9 shows that in this latter case, the baryonic fraction is lower than 0.5 F b in the direct vicinity of the virial radius.

Outside the virial radii, our analysis in the previous section has suggested that the IGM is composed of gas that either never collapses onto haloes or is released by mechanisms such as feedback or merger events. But what are the relative proportions? From Figure 6 we can estimate that 20% of the total accreted baryons is expelled from z = 8 to z = 0 for Halo1 and Halo2, which is quite significant. These results are also in good agreement with Sinha & Holley-Bockelmann (2010) who found from semi-analytic techniques that up to ∼ 25% of gas can be unbind over the course of galaxy assembly. From our analysis, we found that the gas that has been released after being accreted on to haloes is mainly located between Rv and 2Rv at z = 0. If we focus on Halo1 which is quite isolated, we found that the gas component between Rv and 2Rv consists of 55.9%, 58.6% and 60.5% of gas that has been expelled from Halo1 from Sim1a, Sim1b and Sim2 respectively. In other words ∼ 40% of this gas has never collapsed onto Halo1.

It is also worth mentioning that 95% of the gas that has been released through merger events or feedback is in the form of WHIM at z = 0 but only ∼ 65% in the case the gas that has never collapsed into the halo. Hence, the gas that has been ejected through different mechanisms during the galaxy formation and evolution significantly increases the mass of the WHIM.

DISCUSSION AND CONCLUSIONS

In this paper, we have analysed hydrodynamical simulations with an extensive treatment of the physics of baryons in order to study the evolution of the distribution of baryons in a local group-type object. In particular, we have investigated the impact of i) supernova feedback, ii) more powerful sources of feedback at high redshift such as AGN from intermediate mass black holes or hypernovae and iii) accretion/ejection and merger events in the evolution of the mass content of Milky-Way type galaxies. First, when no high early feedback is considered as is the case in Sim1a and Sim1b, we found that the evolution of the mean baryonic fraction value f b estimated at the virial radius of the two main central objects is essentially governed by the relative efficiency at which dark matter and baryons are accreted. Indeed, in the first and early phase of galaxy formation (the cold mode of accretion z ≥ 1.5) during which stars mainly form within the galaxy from infalling cold gas [START_REF] Katz | The IGM/Galaxy Connection[END_REF]Kereš et al. 2005;Ocvirk, Pichon & Teyssier 2008;Dekel et al. 2009) f b is in general higher than the universal value. But in the second phase (the hot mode of accretion), dominated by the accretion of diffuse matter or small satellites, f b tends to decrease and is 10-15% lower than the universal value, 0.166, at z = 0. This decrease is not due to the fact that more baryons (relative to DM) are expelled beyond the virial radius by physical mechanisms such as feedback or merger events. Rather this essentially follows from the relative nature of dissipative gas and non-dissipative dark matter. Specifically, the gas component can be shock-heated at the virial radius, which considerably slows down its further accretion into the halo (see for instance Oñorbe et al. 2007). However, shock-heating will presumably not be a large factor in low mass galaxies, but as mentioned in the introduction, feedbacks such as UV-background and supernovae winds are supposed to play a major role in the effective reduction of the baryon content of these galaxies. Additionally, we found that accretion and particularly merger events can release a fraction of both dark matter and baryons. But again, dark matter that has been expelled can be re-accreted more efficiently than the corresponding baryons, hereby enhancing the decrease of f b inside the virial radius. Consequently, our simulations suggest that the gas located beyond the virial radius is mainly in the form of WHIM that comes from either gas that never collapsed prior to the formation of the proto-galaxy (∼ 40%) or gas expelled by feedback and merger processes (∼ 60%).

Our comparison of Sim1a and Sim1b also suggest that globally increasing somewhat the efficiency of supernovae does not significantly affect the evolution of f b inside the virial radius. In other words, standard SN feedback is not powerful enough to expel the baryons out of the gravitational potential of dark matter haloes of such mass scales. This also means that, if future observations indicate that f b inside the virial radius of Milky Way mass galaxy at z = 0 is relatively close to the universal value, then in order to make clear a distinction between the different feedback models, the stellar mass and the properties of the intergalactic medium may provide better constraints than baryon fractions, as suggested by recent numerical works (Davé, Oppenheimer & Sivanandam 2008;Shen, Wadsley & Stinson 2010;Stinson et al. 2011Stinson et al. ,2012;;Scannapieco et al. 2012;Brook et al. 2012b). But at high redshift the situation is quite different since haloes are less massive and thus feedback mechanisms have a more effective contribution. Indeed for our other scenario (for which a higher fraction of gas is expelled at high redshift (z ∼ 8) by more powerful source of feedback such as AGN from intermediate mass black holes or hypernovae implemented in Sim2) the accretion of gas is strongly affected at high redshift, and drives the mean baryonic value inside the virial radius to a lower values at z = 0. This is one important conclusion of this paper: to get a low value of f b at present time, physical mechanisms able to expel the gas at high redshift will have a stronger impact on the deficit of baryons in the mass budget of Milky Way type-galaxy today than those that expel the gas in the longer, late phases of galaxy formation. Thus, if future observations favour a relatively low value of f b at z = 0, then this "high redshift feedback" scenario is the most probable.

However, it this worth mentioning that the effective contribution of non-gravitational heating such as galactic winds from supernovae in the evolution of f b within the virial radius of Milky Way type galaxy is still a controversial and a debated topic in the recent literature. From our modelling and experiments, we found that increasing the efficiency of supernova feedback leads to the reduction of the stellar mass but has no particular impact on the evolution of f b (< Rv) for system of such mass scales (cf Powell, Slyz & Devriendt 2011). Such statements were also advocated by Anderson & Bregman (2010) from observational constraints on the density of hot gas around the Milky Way. Moreover, using cosmological simulations, Faucher-Giguère, Kereš & Ma (2011) have investigated the net baryonic accretion rates through RV . Although the net baryonic accretion rate is sensitive to galactic outflows, especially for low mass systems, the baryon mass fractions at z = 0 for objects with a mass higher than 10 12 M⊙ and for different SN feedback prescriptions have similar values to those derived in the present study. However, other hydrodynamical simulations from Scannapieco et al. (2008Scannapieco et al. ( , 2009) ) using higher SFR efficiency (c * = 0.1) found that at z = 0, the baryon fractions inside the virial radius are ≤ 0.10 within the virial radius "indicating that a significant amount of baryons has been lost through winds".

AGN feedback is supposed to play an important role in the formation of massive objects such as giant ellipticals in galaxy groups (Tabor & Binney 1993;Ciotti & Ostriker 1997;Silk & Rees 1998). From previous numerical studies, the inclusion of AGN is expected to reduce the stellar mass (see for instance Sijacki et al. 2007;Dubois et al. 2010Dubois et al. , 2012;;Hambrick et al. 2011b) and also increase the mass of WHIM (Cen & Ostriker 2005;Tornatore et al. 2010;Durier & Pacheco 2011;Roncarelli et al. 2012). Our study does not include the effect of AGN feedback at low redshift, but we think that such sporadic events will not have a major contribution in the baryonic fraction value inside the virial radius at present time. Indeed, from observationally motivated constraints, AGN and SN at low to moderate redshift seem not to produce the expected correlations with the baryonic Tully-Fisher relationship (Anderson & Bregman 2010). Moreover, as mentioned above, observations of massive objects reported by Humphrey et al. (2011Humphrey et al. ( , 2012) ) seem to suggest that AGN feedback has no particular impact on the baryonic fraction values within the virial radius. Also, recent hydrodynamical simulations show that AGN and SN feedback may mutually weaken one another's effect by up to an order of magnitude in haloes in the mass range of 10 11.25 -10 12.5 M⊙ (Booth & Schaye 2012). At high redshift, the effects of galactic winds on the IGM are also controversial. In the present paper, our simple minded modelling indicates that galactic winds can expel some gas outside the virial radius and affect significantly the evolution of the baryon content of massive galaxies. Similar conclusions have been recently drawn by Roncarelli et al. (2012) from hydrodynamical simulations suggesting that galactic winds acting at earlier epochs can prevent the IGM to collapse in dense structures. Also, as mentioned above, semi-analytical model from Valiante, Schneider & Maiolino (2012) indicates that the observed outflow of a distant quasar at z = 6.4 (Maiolino et al. 2012) is expected to have considerably depleted the gas content of the host galaxy. However, simulations conducted by Di Matteo et al. (2011) or Shen et al. (2012) found that galactic outflows at high redshift may not be very effective at stopping the cold gas from penetrating the central regions.

In summary, our study indicates that in order to reach lower f b values at z = 0 for Milky Way type galaxies, the eviction of cold gas by feedbacks during the first phase of galaxy formation at high redshifts proves to be crucial. Our study also suggests that these high redshift feedback mechanisms more efficiently reduce galaxy stellar mass in the central part of haloes and thus show better agreement to galaxy observable such as the stellar mass-halo mass relation. Moreover, efficient high redshift feedback processes seem also to be required in order to slow down the growth of galaxies at high redshift and thus to reproduce the observed number density evolution as recently pointed out by Weinmann et al. (2012). The presence of strong early feedbacks is also motivated in disk galaxy formation theory. Indeed, it has been suggested that early outflows could eject low angular momentum gas allowing disk galaxies to form (Binney, Gerhard and Silk (2001). Such hypothesis was confirmed by hydrodynamical simulations which show that outflows can indeed remove preferentially low angular momentum material which can be re-accreted latter with higher angular momentum and forms a disk (Brook et al. 2011(Brook et al. , 2012a)). If such statements are correct, numerical and observational efforts have to be focused towards characterising the respective role of each feedback process on the IGM at high redshift. In particular, are the baryons expelled from protogalaxies or do they never collapse into the gravitational potential well of these objects? In the present paper, we only focus on the first mechanism, but it is not excluded that processes such as photoionization play a major role especially at high redshifts for low mass systems. In this regard, it is interesting to compare our results with those of Guedes et al. (2011) who have studied the formation of late-type spirals from the ERIS simulation. Specifically, they claimed that the collapse of baryons is heavily suppressed at high redshift by the UV background. And it is particularly instructive to see that from their highest mass resolution simulation, the evolution of the baryonic fraction (see their Figure 6 in Guedes et al. 2011) is very similar to those we obtained from Sim2. They also found a similar value at z = 0 namely f b ∼ 0.7 < f b >. Our simulations do include UV background but our mass resolution may be too small to account for this process properly.

More detailed and well-resolved hydrodynamical simulations are required to investigate the evolution of the mass content of galaxies at high and low redshifts. But in order to achieve this goal, one must first improve the constraints of the expected energy injection to the IGM from SN, AGN, UV background (etc...) via observations.

Figure 1 .

 1 Figure1. The evolution of the virial dark matter (DM), gas and stellar mass for Halo1 (first column) and Halo2 (second column) as a function of the time or redshift from Sim1a (green curves), Sim1b (red curves) and Sim2 (blue curves). In the stellar mass panels (fourth line) we also show in dashed lines the evolution of the stellar mass inside the galaxy radius defined as R gal = 0.1 R V . Note that Halo1 undergoes several major merger episodes while Halo2 grows through smooth accretion. Moreover, different feedback process prescriptions lead to different baryon mass contents inside Rv .

Figure 2 .

 2 Figure 2. The evolution of the star formation rate (SFR) for Halo1 (left panel) and Halo2 (right panel) as labelled.

Figure 3 .

 3 Figure3. The projected baryonic fraction at z = 5 (first column) z = 2 (second column), z = 0.5 (third column) and z = 0 (fourth column) from Sim1b (first line) and Sim2 (second line). Dashed circles show virial radii. High f b value regions are clearly visible around galaxies (or proto-galaxies) and in the filaments. This trend seems to be more pronounced in Sim2 although the difference between the two simulations is vanishing at low redshift.

Figure 4 .

 4 Figure 4. The evolution of the baryonic fraction f b estimated at the virial radius for Halo1 (left panel) and Halo2 (right panel) as a function of cosmic time. In each panel, red, green and blue lines correspond to values derived from Sim1a, Sim1b and Sim2 respectively. The horizontal dashed line corresponds to the universal fraction. Note that the early feedback model induces a lower (f b ∼ 0.12) baryon fraction at present time.

Figure 5 .

 5 Figure 5. Distributions of the differential baryonic fraction, f b (R) (first column), and the integral baryonic fraction f b (< R) (second column) with respect to the radius, R, for Halo2 and for three different redshifts, namely z = 4.7 (green lines), z = 0.8 (blues lines) and z = 0 (red lines).

Figure 7 .

 7 Figure7. An illustration of the evolution of the distribution of baryons and dark matter during and after given merger events. Panels at z = 0 is 1 × 1 Mpc across while the other ones are 0.3 × 0.3 Mpc across. The white dashed lines represent the virial radii. In the lower panels, we show the evolution of the mass of baryons and DM outside the virial radius of the host halo (left panel) and the evolution of the baryonic fraction estimated from these 3 haloes inside and outside the virial radius (right panel). One can see that a non negligible number of DM particles are re-accreted inside Rv (from t ∼ 4.5 to t ∼ 10 Gyr) which is not the case for baryons. As a result, the final values of f b (< Rv ) and f b (> Rv) are respectively lower and higher than initial ones.

Figure 8 .

 8 Figure8. The evolution of the IGM components within a sphere centred on the local Universe type haloes and of comoving radius 1.5 Mpc/h: cold gas (T < 10 5 K, dashed lines), WHIM (10 5 K < T < 10 7 K, solid lines) and stars (dotted lines). Green, red and blue lines correspond to Sim1a, Sim1b and Sim2 respectively. bf The evolution of the hot gas (T ≥ 10 7 K) is not shown since this component represents less than 1% of the total baryon mass budget in the selected region.

Figure 9 .

 9 Figure 9. The projected baryonic fraction at z = 0 derived from Sim1a when the WHIM is not taken into account in the estimation of f b .

Table 1 .

 1 Fractions γ of the energy released by supernovae

		Sim1a	Sim1b	Sim2
	γ	0.1	1.0	50 (8.1 ≥ z ≥ 7.7) 0.1 elsewhere

Table 2 .

 2 Physical properties of the 2 main central objects of each simulation derived at the virial radius Rv and z = 0.

			Halo1			Halo2
		Sim1a	Sim1b	Sim2	Sim1a	Sim1b	Sim2
		(γ = 0.1) (γ = 1.0)		(γ = 0.1)	(γ = 1.0)
	m DM (×10 11 h -1 M ⊙ )	7.58	7.58	7.23	7.43	7.40	8.17
	mgas (×10 11 h -1 M ⊙ )	0.30	0.44	0.27	0.38	0.44	0.45
	mstars (×10 11 h -1 M ⊙ )	1.08	0.93	0.71	0.84	0.78	0.68
	Mv (×10 11 h -1 M ⊙ )	8.96	8.95	8.21	8.65	8.62	9.30
	Rv (Mpc)	0.282	0.282	0.274	0.279	0.279	0.286
	f b	0.154	0.153	0.119	0.141	0.141	0.121

c 0000 RAS

c 0000 RAS, MNRAS 000, 000-000

Acknowledgements I. J. acknowledges support from the Global Internship program by the National Research Foundation of Korea (NRF). We warmly thank the referee for an insightful review that considerably improved the quality of the original manuscript. We warmly thank D. Le Borgne, J. Devriendt, Y. Dubois, R. Gavazzi, T. Kimm, G. Mamon, J.A. de Freitas Pacheco and T. Sousbie for interesting discussions. We also thank D. Munro for freely distributing his Yorick programming language (available at http://yorick.sourceforge.net/) which was used during the course of this work. This work was carried within the framework of the Horizon project (http://www.projet-horizon.fr).


Hubert Garavel

Radu Mateescu

Wendelin Serwe
email: wendelin.serwe@inria.fr

Radu Garavel

Mateescu

PASM / PDMC

Keywords: asynchronous systems, distributed verification, labeled transition system, model checking, on-the-fly verification, process calculus 1 {Hubert

Distributed verification uses the resources of several computers to speed up the verification and, even more importantly, to access large amounts of memory beyond the capabilities of a single computer. In this paper, we describe the distributed verification tools provided by the CADP (Construction and Analysis of Distributed Processes) toolbox, especially focusing on its most recent tools for management, inspection, and on-the-fly exploration of distributed state spaces. We also report about large-scale experiments carried out using these tools on Grid'5000 using up to 512 distributed processes.

Introduction

When analyzing concurrent systems using explicit-state verification methods in an action-based setting, the semantic models for state spaces are Labeled Transition Systems (LTSs). One approach consists in first building the LTS of the concurrent system under study; this LTS can then be minimized modulo a bisimulation relation to increase the efficiency of further analyses, such as model checking, equivalence checking, visual checking, etc. An alternative approach is to perform these analyses on-the-fly, i.e., during the construction of the LTS, so as to detect errors without constructing the entire LTS first. The latter approach is more suitable for early verification steps, where bugs are frequent and can be quickly detected, while the former approach is more efficient in the late phases of the design process, where the entire LTS must be explored to ensure its correctness.

Due to the state explosion phenomenon (prohibitive size of the LTS for systems containing many concurrent processes and/or complex data types), LTS generation can become a bottleneck in the verification process. In such case, if the LTS is too large to be constructed on a single machine, one may resort to distributed computing infrastructures, such as clusters and grids, which increase by several orders of magnitude the amount of memory available. The CADP verification toolbox [START_REF] Garavel | Cadp 2011: A toolbox for the construction and analysis of distributed processes[END_REF] exploits this possibility by providing several tools for distributed verification, in particular DISTRIBUTOR and BCG MERGE [START_REF] Garavel | Parallel state space construction for modelchecking[END_REF][START_REF] Garavel | Distributor and bcg merge: Tools for distributed explicit state space generation[END_REF]. These tools respectively enable to construct a partitioned LTS (i.e., split into several fragments, each stored in a separate file, possibly on a different machine) and to convert it into a monolithic LTS (i.e., stored in a single file). To scale up the verification capabilities, it is sometimes beneficial to avoid the construction of a monolithic LTS and instead work as long as possible with a partitioned LTS.

In this paper, we present the CADP tools (some of which have been recently added) for manipulating partitioned LTSs. All these tools are based on distributed algorithms and are implemented using standard network communication primitives available on most clusters and grids. We also report about large-scale multi-cluster experiments on distributed LTS generation and onthe-fly reduction. These experiments have been carried out on the Grid'5000 computing infrastructure [START_REF] Cappello | Grid'5000: A Large Scale, Reconfigurable, Controlable and Monitorable Grid Platform[END_REF] using up to 512 distributed processes, and provide insight about the performance gains and scalability when the number of distributed processes increases.

The paper is organized as follows. Section 2 briefly describes the software library of CADP for handling network communications. Section 3 defines the PBG format for partitioned LTSs. Sections 4 and 5 present the new tools for creating and manipulating PBG files. Section 6 describes a new tool enabling an on-the-fly exploration of PBG files. Section 7 gives experimental measures about the use of these new tools on the Grid'5000 computing infrastructure [START_REF] Cappello | Grid'5000: A Large Scale, Reconfigurable, Controlable and Monitorable Grid Platform[END_REF] for the distributed generation and on-the-fly reduction of large LTSs. Section 8 gives a brief overview of related work. Finally, Section 9 concludes the paper and suggests directions for future work.

The Network Communication Library

A typical distributed application developed using CADP consists of N + 1 POSIX processes, namely N workers, possibly running on remote machines or on different processors/cores of multi-processor/multi-core machines, and one master, which supervises the execution of workers, runs on the user frontal machine, and interacts (inputs/outputs) with the user. The communications rsh = ssh -q rcp = scp -q connect_timeout = 30 edel-42.grenoble.grid5000.fr directory=/home/serwe/ between these processes rely on a dedicated library named caesar network 1 (or network 1 for short), which currently implements two communication topologies: star (all workers connected only to the master) and fully-connected (complete graph between workers and master).

The configuration of the grid is described as a text file in the GCF (Grid Configuration File) format [START_REF] Cadp | [END_REF], which specifies the computing nodes on which workers have to execute, the way to access them (connection protocols and parameters for file transfer), and the directories in which workers have to run. A GCF format is application-independent and can be used for launching several distributed applications running on the same computing nodes. Figure 1(a) shows an example of an GCF file with four workers.

To ensure maximal portability, the network 1 library is purposely based on standard operating system primitives (namely, TCP/IP sockets) and standard remote access protocols (namely, rsh or ssh), so that the distributed tools of CADP do not require additional communication libraries (such as MPI) to be specifically installed on each remote machine. Furthermore, the network 1 library supports seamless multi-core, intra-cluster, and inter-cluster communication: apart from performance aspects, the end-user sees no functional difference between an application running on several cores of a single machine and an application running on several machines belonging to one or many clusters.

The PBG Format

To represent LTSs, the CADP toolbox provides the BCG (Binary Coded Graphs) file format and its associated software libraries. A BCG file stores the states, labels, and transitions of an LTS in a compact way using binary encoding and dedicated compression schemes that enable efficient representation and manipulation. BCG files can be handled using the existing CADP tools (e.g., inspection, visualization, label renaming, bisimulation minimization, on-the-fly exploration, etc.) or using custom tools developed using the CADP libraries for reading and writing BCG files.

When dealing with distributed verification tools and LTSs stored on several machines, a single BCG file is no longer sufficient. The PBG (Partitioned BCG Graph) format [START_REF] Cadp | [END_REF][START_REF] Cadp | [END_REF] addresses this problem. This format is an outcome of the SENVA co-operation2 between the former SEN2 team of CWI and the former VASY team of Inria Grenoble. Specifically designed for the purpose of distributed verification, the PBG format implements the theoretical concept of Partitioned LTS introduced in [START_REF] Garavel | Parallel state space construction for modelchecking[END_REF] and provides a unified access to an LTS distributed over a set of remote machines. A PBG file gathers a collection of BCG files, called fragments (one fragment per worker), which can be stored either in separate directories located on the (possibly remote) machines on which workers execute, or on a common file system shared (e.g., using NFS or Samba) by all workers. Taken altogether, these fragments form a partition of the LTS, the states and transitions of which are distributed across the various fragments as specified in [START_REF] Garavel | Parallel state space construction for modelchecking[END_REF], each fragment storing a set of states and the transitions going into these states. Note that, taken individually, each fragment is meaningless; for instance, it may be a disconnected graph, which is never the case with an LTS representing the reachable state space of a concurrent system.

Concretely, a PBG file is a text file containing references to the fragments and the GCF file used for constructing the partitioned LTS. The example of a PBG file shown in Figure 1(b) corresponds to an LTS partitioned in four fragments, the first of which contains the initial state of the LTS. For each fragment, the PBG file lists the number of states of the fragment and the files associated to this fragment (i.e., a BCG file and a log file containing error messages that may have been issued when building the fragment). There exists a simple code library for reading and writing PBG files.

It is worth noticing that the PBG format uses a number of fragments linear in the number of computing nodes; this is better than the competing SVC format [START_REF] Blom | Compressed and Distributed File Formats for Labeled Transition Systems[END_REF], whose number of fragments is quadratic in the number of computing nodes.

Tools for PBG Creation

A partitioned LTS in the PBG format can be generated using the DISTRIB-UTOR tool [START_REF] Garavel | Parallel state space construction for modelchecking[END_REF][START_REF] Garavel | Distributor and bcg merge: Tools for distributed explicit state space generation[END_REF] of CADP. The tool works by launching several workers on (local and remote) computing nodes specified by a GCF file. Each worker is in charge of generating a fragment of the LTS, which is stored as a BCG file; a static hash function determines which state is explored by which worker. Upon termination, DISTRIBUTOR produces a PBG file gathering all these Fig. 2. Overview tab monitoring LTS generation using twelve workers distributed over six nodes of six clusters (adonis, edel, genepi, granduc, sagittaire, and suno) of Grid'5000, geographically located in Grenoble, Luxembourg, Lyon, and Sophia-Antipolis. In column "Variation", a green (respectively, orange) box indicates that the number of remaining states is increasing (respectively, decreasing or stable), and a red box indicates that the worker is idle; boxes are red for both workers on adonis-6 and the first worker on suno-9, orange for the second worker on edel-71, and green for all other workers.

fragments. The progression of the distributed computation can be monitored in real-time (see Figure 2).

Besides the obvious advantages brought by distributed LTS generation (speeding up the generation and increasing the amount of memory far beyond what is made available by a single machine), DISTRIBUTOR also provides onthe-fly reductions (τ -compression and τ -confluence) that preserve branching bisimulation and may be useful when dealing with large LTSs that cannot be generated and minimized on a sequential machine.

Tools for PBG Inspection and Manipulation

Several tools are currently available in CADP for handling PBG files:

• PBG MERGE (previously called BCG MERGE [START_REF] Garavel | Distributor and bcg merge: Tools for distributed explicit state space generation[END_REF]) converts a partitioned LTS represented in the PBG format into a monolithic LTS stored in a BCG file. The LTS fragments are merged into a single file, in which states are given a contiguous numbering that improves compactness of the resulting BCG file.

• PBG CP, PBG MV, and PBG RM are new tools for copying, moving, and removing PBG files, keeping in mind that the fragments of these PBG files may be disseminated on a number of remote machines, possibly located in different countries. These tools facilitate standard operations on PBG files and maintain consistency during these operations.

• PBG INFO is a new tool for inspecting PBG files. It currently provides several functionalities, such as consistency checking (i.e., existence and readability of all fragment files), calculation the size of the corresponding LTS (number of states and transitions), display of the list of labels, and concatenation of remote log files (this is useful, e.g., to understand the reason why a PBG generation fails, and to compute global statistics about CPU and memory usage by the workers).

Tools for PBG On-the-fly Exploration

On-the-fly verification is an approach to fight state explosion by verifying an LTS during its construction rather than constructing the LTS first and verifying it afterward. So doing, on-the-fly verification may enable an early detection of errors even in presence of state explosion. The OPEN/CAESAR [START_REF] Garavel | Open/Caesar: An open software architecture for verification, simulation, and testing[END_REF] environment of CADP provides a modular architecture for on-the-fly verification tools, which separates the language-dependent aspects (translation of a concurrent system description into an LTS) from the language-independent aspects (forward exploration of an LTS, e.g., for verification). OPEN/CAESAR defines a generic API for representing an LTS by its transition relation. OPEN/CAESAR also contains a set of libraries implementing various primitives and data structures dedicated to on-the-fly graph exploration (hash tables, stacks, etc.).

CADP currently provides OPEN/CAESAR-compliant compilers for several high-level description languages (LOTOS, LNT, and FSP) and low-level state machine formats (BCG, EXP, and SEQ). These compilers implement the OPEN/CAESAR API to explore the corresponding LTSs on the fly.

CADP also provides a set of on-the-fly verification tools based on the OPEN/CAESAR API: model checking MCL formulas (EVALUATOR), bisimulation checking (BISIMULATOR), partial order reduction (REDUC-TOR), random exploration (EXECUTOR), regular sequence searching (EX-HIBITOR), steady-state Markov chain simulation (CUNCTATOR), etc. All these tools can be applied to any description in an input language for which an OPEN/CAESAR-compliant compiler exists.

The newly developed PBG OPEN tool is an OPEN/CAESAR-compliant compiler for the PBG format. The main advantage of PBG OPEN is that it can use the memory of several machines to store the transition relation of a partitioned LTS. Therefore, PBG OPEN can explore on-the-fly large partitioned LTSs that could not be explored using other tool combinations.

PBG OPEN (see Figure 3) is a distributed tool consisting of a master and several workers, each associated to a BCG fragment referenced in the PBG file. Each worker is responsible for opening its fragment, initializing The initialization phase consists in normalizing the transition labels of the fragments by assigning unique label numbers across all workers. This is necessary because the same label may be numbered differently in different fragments. During initialization, each worker sends its list of labels to the master, which assigns a unique number to each label and then sends back to each worker globally unique numbers for these labels. This preliminary step avoids the performance overhead that would occur if the master had to renumber labels in all transitions.

Each transition is represented by a triple s, a, s ′ , where s, a, and s ′ are numbers encoding the source state, label, and target state; a is now a globally unique label number common to all fragments. Notice that fragments store the incoming transitions, i.e., a given worker stores all triples s, a, s ′ for a given s ′ , whereas the triples s, a, s ′ for a given s may be distributed between workers.

When requested by an on-the-fly exploration tool to compute all transitions going out of a given state s, the master forwards the request to all workers. Each worker retrieves (using the BCG primitives) its transitions going out of state s and sends them back to the master. The master explores the transitions received from the workers in their order of arrival, which is nondeterministic.

Given that some on-the-fly analysis tools (e.g., the REDUCTOR tool) often explore several times the transitions going out of the same state, PBG OPEN implements memoization using a cache. After receiving the list L of transitions going out of a state s, the master stores the couple s, L in the cache. When the cache is full, a couple s ′ , L ′ present in the cache is selected according to the LRU (Least Recently Used ) strategy, and is replaced by the new couple s, L . In general, neither the number of transitions in L nor the maximal size of L are known in advance (they could be computed using a preliminary LTS traversal, but this would be too costly). Two variants of the cache have been implemented in PBG OPEN using the cache 1 library initially developed for state space caching [START_REF] Mateescu | Hierarchical adaptive state space caching based on level sampling[END_REF]. These variants differ in the way the couples s, L are stored: (a) the variable-size variant stores source states s in the cache and transition lists L in the heap, outside of the cache; (b) the fixed-size variant stores both source states s and transition lists L inside the cache; if L is long, it may replace one, several, or even all entries of the cache; if L is too long for the cache, s, L will not be stored in the cache. Variant (a) is simpler to implement, but variant (b) guarantees a statically bounded amount of memory.

Applications and Experiments

We experimented these tools on examples taken from a case study on the formal verification and performance analysis of various mutual exclusion protocols [START_REF] Mateescu | Model Checking and Performance Evaluation with CADP Illustrated on Shared-Memory Mutual Exclusion Protocols[END_REF] specified in the LNT language. 3 We instantiated each protocol for four or five processes competing for the critical section and generated the corresponding LTSs4 directly: so doing, we observed larger state spaces (see Table 1) than those obtained by the compositional approach described in [START_REF] Mateescu | Model Checking and Performance Evaluation with CADP Illustrated on Shared-Memory Mutual Exclusion Protocols[END_REF].

The experiments ran on the Grid'5000 clusters geographically located in Grenoble, Luxembourg, Lyon, Reims, and Sophia-Antipolis, all of which were equipped with two processors per node (i.e., machine or server) and the same operating system (Debian Linux 6.0 "Squeeze"). The clusters differ however in the number of nodes, the type of processor, the number of cores per processor, and the amount of RAM per node (see Table 2 for details). In each cluster, the nodes are connected by 1 GBit/s links; all the clusters in Grenoble are connected to the same switch. Communication between sites uses a dedicated 10 GBit/s link. Each site provides a unique resource manager for all clusters of the site. Whenever possible, we used a separate node for the master and dedicated one core to each worker, i.e., we launched at most n workers on a node with a total of n cores. Execution time and memory consumption were measured using memtime 5 , and, if possible, averaging several executions. In execution time, we include setting up the workers (creating working directories, copying files, etc.). Concerning memory consumption, we measure for each process (master and workers) the maximal amount of memory required by this process during its execution; we call total memory the sum of these maximal amounts for all processes, and we call peak memory the greatest value among all these maximal amounts.

Because access to computing nodes of Grid'5000 is granted by resource managers, it took extra efforts to run a lot of experiments under the same conditions because we have no control on distribution of allocated nodes over the different clusters and/or switches, overall load on file servers and communication network, etc. For instance, we observed significant variation in execution time; to ensure consistent figures, we excluded the 20% extreme values before computing the average. Also, due to the rather small set of examples sharing similar characteristics (such as the average number of outgoing transitions or the number of labels), the experiments reported here illustrate tendencies; results may vary for other types of LTSs and/or grid configurations.

Performance Study of Distributed State Space Generation

For each mutual exclusion protocol example, we used DISTRIBUTOR to generate the corresponding PBG using a single cluster, multiple clusters at the same site, and multiple clusters at different sites. Figures 4 and5 give memory consumption and execution time speedup for up to 512 workers; Figure 4 is for distributed execution and Figure 5 is for multi-core execution on a single server. Note that many axes of the figures use a logarithmic scale. workers reduces the peak memory consumption6 , i.e., the maximum amount of memory required by the master or any of the workers: when using two workers, we observe reductions up to 40%. On the three largest examples ("Peterson", "MCS", and "Dijkstra"), this memory reduction enables LTS generation on the genepi cluster that provides only 8 GB RAM per node (using two or four nodes).

Peak memory consumption.

For each example there is an "optimal" number of workers which minimizes peak memory consumption: using more workers increases peak memory consumption. This can be explained by the constant memory requirements of the CADP communication library (table of nodes, communication buffers, sockets, etc). Figure 4(a) shows that the smaller the LTS, the smaller the optimal number of workers (for very small examples, distributed verification makes no sense and sequential LTS generation performs betters).

Figure 4(a) was obtained by experimenting with workers distributed over the three clusters in Grenoble; the figures observed using either a single cluster or clusters on different sites are the same.

Speedup. Figures 4(b) to 4(d)

show the speedup for different combinations of clusters. Workers use their local disks rather than a shared file system. We take as the reference for the speedup the sequential GENERATOR tool, running on a machine of the genepi cluster for Figure 4(b) and running on a machine of the edel cluster for Figures 4(c the large examples "Dijkstra", "MCS", and "Peterson" require more memory than available on single a machine: thus we take as reference the execution of DISTRIBUTOR with the smallest number of workers sufficient to generate the LTS (namely, two for "Peterson" and four for "Dijkstra" and "MCS").

The observed speedups depend on the communication cost, but the tendencies are similar. All three graphs show an almost linear speedup as long as the number of workers is lower than the "optimum" value mentioned above. If the number of workers becomes too large, speedup drops: for small examples and many workers, distributed execution even becomes longer than sequential execution (speedup lower than one) due to the overhead of setting up the distributed application (currently implemented by sequentially copying all necessary files to the remote machines).

Concerning execution time (rather than speedup): the difference in execution time between the slowest and fastest examples for a given number of workers is reduced when the number of workers increases. For instance, this difference is more than one hour for sequential execution and only two minutes with 256 workers.

Multi-core execution. Figure 5 shows the speedup and total memory observed when running DISTRIBUTOR on one 24 core server, with at most one worker per core.

Concerning speedup, the almost linear speedups of Figure 5(a) show that DISTRIBUTOR also performs well on multi-cores even if the socket-based network 1 library used by DISTRIBUTOR has not been specifically optimized for shared-memory architectures. When increasing the number of workers beyond the number of cores, the speedup drops drastically; because this is expected, it is not shown in Figure 5(a).

Concerning total memory consumption, Figure 5(b) shows that the memory consumption increases linearly with the number of workers, which can be explained by the constant memory cost per worker. To measure the overall performance of PBG OPEN on a particularly demanding on-the-fly application, we selected the REDUCTOR tool of CADP, which performs τ -confluence reduction using a sequential on-the-fly algorithm. We first generated PBG files using DIS-TRIBUTOR and then applied PBG OPEN and REDUCTOR to reduce these distributed state spaces with respect to τ -confluence, after hiding all actions except entry and leave of the critical and non-critical sections. In such "stress tests" conducted on all but the largest examples7 , the complete state space, including very large components of τ -transitions, must be explored on-the-fly (as can be deduced from the small size of the state spaces after minimization for branching bisimulation -see the last two columns of Table 1).

Figure 6 summarizes the results of the experiments using the granduc cluster and with up to 128 workers. Figure 6(a) shows the overall peak memory consumption, i.e., the maximum amount of memory used by DIS-TRIBUTOR and PBG OPEN/REDUCTOR; Figure 6(b) shows the overall execution time, i.e., the sum of the execution time of DISTRIBUTOR and PBG OPEN/REDUCTOR (we did not compute the speedup because these experiments concern a sequential algorithm operating on a distributed state space). For these experiments, the mechanism of PBG was because the effects caching will be studied in section 7.3. The numbers for one worker correspond to execution of a sequential done directly at the level of the LNT source language, using the LNT.OPEN/REDUCTOR tools DISTRIBUTOR.

As for medium-size examples, with two and workers, consumption is lower than with LNT.OPEN/REDUCTOR because PBG OPEN uses eight bytes per state, which is less than the state sizes with LNT mentioned in the second column of Table 1. Because the master executing the sequential on-the-fly algorithm of REDUCTOR is the bottleneck, the peak memory consumption is almost independent from the number of workers; however, memory consumption increases for all examples if the number of workers gets too high. The reduction of the peak memory comes at the price of a (significant) increase in execution time due to the communication latency when collecting outgoing transitions from the distributed fragments: we observed a drop in CPU usage of REDUCTOR from 100% for LNT.OPEN down to 40% for PBG OPEN (even down to 1% if the fragments are located on geographically distant sites).

Alternative τ -confluence reduction. We then compared the aforementioned τconfluence reduction using DISTRIBUTOR, OPEN, REDUCTOR against two approaches, namely:

• the sequential of REDUCTOR at the level of the LNT source language (i.e., using and on a single node),

• the distributed generation of PBG file using DISTRIBUTOR, followed by an LTS using PBG MERGE, and execution of the resulting file (i.e, using BCG OPEN and REDUCTOR on a single node).

Table 3 shows execution time and peak memory consumption of these three approaches, confirming the claims of section 6: for all but the smallest examples, PBG OPEN requires the least amount of memory. The direct connection to the source language requires the most memory, because states are larger (in particular when the LNT program contains complex data structures) than for the exploration of a PBG or BCG (where each state is represented by a number). BCG OPEN, which loads the whole LTS in memory, requires more memory than PBG OPEN, in which each worker uses its own memory to load only a fragment of the partitioned LTS.

Double τ -confluence reduction.

For some examples, we also experimented the combination of two successive on-the-fly τ -confluence reductions, applying REDUCTOR on a PBG generated using DISTRIBUTOR with activated onthe-fly τ -confluence reduction. This double reduction yields LTSs that are up to a factor 1.9 smaller than those obtained by calling REDUCTOR only (but still about 500 times larger than the minimized LTS).

Figure 7 shows the results of these experiments with up to 64 nodes of the granduc cluster. Contrary to Figure 6(a), we observe that the peak memory consumption decreases when the number of workers increases. However, the peak memory consumption is higher than for the combination of REDUCTOR and DISTRIBUTOR without τ -confluence; the only exception is "CLH" for 64 or more workers. Indeed, because activating the on-the-fly τ -confluence reduction of DISTRIBUTOR yields a significantly smaller PBG, the bottleneck concerning peak memory consumption is not the master running REDUC-TOR, but the workers running DISTRIBUTOR with activated τ -confluence reduction. We observed that for very large components of τ -transitions and few workers, DISTRIBUTOR with activated τ -confluence reduction requires almost as much memory as a single sequential execution of REDUCTOR, because the complete component is explored by (at least) one worker.

Note that the size (i.e., number of states and transitions) of the reduced partitioned LTS depends on the order in which states are explored; this order depends on both the number of workers and communication latencies.

Performance Study of Caching in PBG OPEN

Fixed-Size Caching. Figures 8 and9 show the effects of increasing the cache size in PBG OPEN up to one million entries, when using the fixed-size cache As expected, caching reduces execution time if the cache is large enough to hold all required transitions. This can be seen on small examples, in particular those obtained by DISTRIBUTOR with τ -confluence reduction: these examples have significantly smaller components of τ -transitions, so one observes reduction in execution time already for small cache sizes.

On the hand, caching peak memory consumption, except for those examples where DISTRIBUTOR requires more memory than RE-DUCTOR. For small LTSs and many cache entries, memory consumption may become even larger than for sequential execution -namely, when the cache stores the entire LTS.

Comparison of Fixed-and Variable-Size Caching. Comparing both cache implementations with rigorously the same cache size is not simple because of the varying sizes of the lists of outgoing transitions. It is difficult to predict the memory consumption of the variable-size cache, as well as the number of states that can be stored in a fixed-size cache (because each fixed-size cache entry can store either a state or a transition). To simplify our comparisons, we chose to use the same number of cache entries (which is only an approximation for comparing caches with the same amount of memory).

Figure 10 compares both implementations for four selected examples, using two workers running PBG OPEN/REDUCTOR on the granduc cluster. We observed a significant difference only for "Peterson tree", where the variablesize cache is faster, but requires more memory.

Related Work

There are other approaches to parallel and distributed verification. For instance, DIVINE [START_REF] Barnat | DiVinE: Parallel distributed model checker[END_REF] supports distributed LTL model checking on both multi- core architectures [START_REF] Barnat | Scalable Multi-core LTL Model-Checking[END_REF], clusters, and grids [START_REF] Verstoep | Efficient Large-Scale Model Checking[END_REF]. The parallel extension of Spin [START_REF] Holzmann | Parallelizing the Spin Model Checker[END_REF] also supports multi-core LTL model checking, including partial order reduction. Another example is PREACH [START_REF] Bingham | Industrial Strength Distributed Explicit State Model Checking[END_REF], which enables distributed reachability analysis of Murϕ models. PREACH is built on top of the original sequential Murϕ code, implementing distribution and communication in the distributed functional programming language Erlang. Note that flow control credit mechanism that is crucial for scalability of PREACH is built-in in our communication library, as it is provided by TCP/IP.

A key difference between these approaches and our approach is that CADP is action-based rather than state-based. Also, CADP offers tools for generating and handling distributed LTSs, which enables a wide spectrum of verification techniques to be used, including both model checking and equivalence checking (bisimulations). Let us also mention the work of Eric Madelaine and colleagues, who used DISTRIBUTOR and BCG MERGE for the generation and on-the-fly reduction of large LTSs using the PACAGrid8 infrastructure [14,2,1].

Conclusion

We presented the latest distributed verification tools recently added to the CADP toolbox in order to manipulate partitioned LTSs represented as PBG files. We experimented these new tools, together with the DISTRIBUTOR and BCG MERGE tools previously available in CADP, on a large-scale grid involving several clusters geographically located in different places and different countries. Our experiments were intended to push the PBG machinery to its limits by using hundreds of workers and to study how this influences performance and scalability.

Our experiments confirm the finding of [START_REF] Garavel | Parallel state space construction for modelchecking[END_REF] that distributed state space generation using DISTRIBUTOR scales well up to the point where "too many cooks spoil the broth".

We observed that the on-the-fly exploration of partitioned LTSs using PBG OPEN avoids the memory bottleneck faced when using a single machine, and thus enables to take advantage of clusters and grids to handle much larger problems (as the peak amount of memory needed is roughly divided by the number of computing nodes).

Conversely, our experiments showed that PBG OPEN may be slower than a sequential implementation (between two and nine times) when the on-thefly application (e.g., REDUCTOR) is an intrinsically sequential state space exploration algorithm. We observed a significant impact of communication latencies between computing nodes, a problem that we addressed by introducing caches in PBG OPEN.

We did not measure the consumed (peak) bandwidth: similar to the observations reported in [START_REF] Bingham | Industrial Strength Distributed Explicit State Model Checking[END_REF], we never found bandwidth to be a bottleneck. The heterogeneity of the clusters clearly impacts the experiments. On the one hand, the worker with the smallest amount of available memory determines if the computation succeeds (because the static hash function distributes states uniformly over the workers). On the other hand, measuring speedup becomes imprecise if one worker runs on a faster or slower processor than the sequential reference execution. However, we found it interesting to experiment how a static work distribution can cope with this heterogeneity.

This work can be pursued along several directions. Firstly, the network communication library of CADP could be optimized to support parallel (rather than sequential) file transfers during the initialization, which would remove the overhead observed when using hundreds of workers. Secondly, static loadbalancing techniques could be investigated by specializing the static hash function used for assigning states to workers. Thirdly, one could experiment PBG OPEN with a truly parallel OPEN/CAESAR application (such as DISTRIBUTOR with its τ -confluence reduction algorithm, or the distributed version of the EVALUATOR 4.0 [START_REF] Mateescu | A model checking language for concurrent value-passing systems[END_REF] on-the-fly model checker) rather than a sequential application (such as the REDUCTOR used in the present paper).

Fig. 1 .

 1 Fig. 1. Sample GCF and PBG files describing a partitioned LTS

Fig. 3 .

 3 Fig.3. Architecture of an OPEN/CAESAR on-the-fly application built using PBG OPEN and the partial order reduction tool REDUCTOR of CADP label information, and answering the master requests to compute the outgoing transitions of states belonging to that fragment.The initialization phase consists in normalizing the transition labels of the fragments by assigning unique label numbers across all workers. This is necessary because the same label may be numbered differently in different fragments. During initialization, each worker sends its list of labels to the master, which assigns a unique number to each label and then sends back to each worker globally unique numbers for these labels. This preliminary step avoids the performance overhead that would occur if the master had to renumber labels in all transitions.Each transition is represented by a triple s, a, s ′ , where s, a, and s ′ are numbers encoding the source state, label, and target state; a is now a globally unique label number common to all fragments. Notice that fragments store the incoming transitions, i.e., a given worker stores all triples s, a, s ′ for a given s ′ , whereas the triples s, a, s ′ for a given s may be distributed between workers.When requested by an on-the-fly exploration tool to compute all transitions going out of a given state s, the master forwards the request to all workers. Each worker retrieves (using the BCG primitives) its transitions going out of state s and sends them back to the master. The master explores the transitions received from the workers in their order of arrival, which is nondeterministic.Given that some on-the-fly analysis tools (e.g., the REDUCTOR tool) often explore several times the transitions going out of the same state, PBG OPEN implements memoization using a cache. After receiving the list L of transitions going out of a state s, the master stores the couple s, L in the cache. When the cache is full, a couple s ′ , L ′ present in the cache is selected according to the LRU (Least Recently Used ) strategy, and is replaced by the new couple s, L . In general, neither the number of transitions in L nor the maximal size of L are known in advance (they could be computed using a preliminary LTS traversal, but this would be too costly). Two variants of the cache have been implemented in PBG OPEN using the cache 1 library initially developed for state space caching[START_REF] Mateescu | Hierarchical adaptive state space caching based on level sampling[END_REF]. These variants differ in the way

Figure 4 (Fig. 4 .

 44 Fig. 4. Distributed state space generation; all axes in logarithmic scale

  ) and 4(d). However, on Figure4(b),

Fig. 5 .

 5 Fig. 5. State space generation on a single machine of the stremi cluster

Fig. 6 .

 6 Fig. 6. On-the-fly τ -confluence reduction using DISTRIBUTOR and PBG OPEN/REDUCTOR

Fig. 7 .Fig. 8 .

 78 Fig.7. Two successive on-the-fly τ -confluence reductions using DISTRIBUTOR with option τ -confluence reduction followed by PBG OPEN/REDUCTOR

Fig. 9 .

 9 Fig. 9. Cache with increasing number of cache entries: borderline, granduc, and suno clusters, eight workers implementation on a cluster (four workers) and several clusters on different sites (eight workers); figures for other numbers of workers show the same tendencies.As expected, caching reduces execution time if the cache is large enough to hold all required transitions. This can be seen on small examples, in particular those obtained by DISTRIBUTOR with τ -confluence reduction: these examples have significantly smaller components of τ -transitions, so one observes reduction in execution time already for small cache sizes.On the hand, caching peak memory consumption, except for those examples where DISTRIBUTOR requires more memory than RE-DUCTOR. For small LTSs and many cache entries, memory consumption may become even larger than for sequential execution -namely, when the cache stores the entire LTS.

Fig. 10 .

 10 Fig. 10. Comparison of cache implementations for increasing number of cache entries

See http://vasy.inria.fr/senva

These specifications will be included as examples in the next stable release of CADP.

In fact, these LTSs are interactive Markov chains[START_REF] Hermanns | Interactive Markov Chains and the Quest for Quantified Quality[END_REF], which can be seen as particular forms of LTSs.

For the total memory consumption, see Figure5(b).

The long execution times exceed the 62 hour limit of Grid'5000 jobs and are the reason why we did not experiment with the larger examples.

See http://proactive.inria.fr/pacagrid

See http://www.grid5000.fr

Acknowledgments.

The experiments presented in this paper were carried out using the Grid'5000 experimental testbed 9 built by Inria with support from CNRS, RENATER, several Universities, and other funding bodies. We are also grateful to Iker Bellicot, Nicolas Descoubes, Jerôme Fereyre, Yann Genevois, and Rémi Hérilier for their contribution in testing and bug-hunting of the CADP distributed verification tools.