Milena Nowakowska

Olivier Herbinet

Anthony Dufour

Pierre-Alexandre Glaude
email: pierre-alexandre.glaude@univ-lorraine.fr

Detailed kinetic study of anisole pyrolysis and oxidation to understand tar formation during biomass combustion and gasification

Keywords: Anisole, Pyrolysis, Oxidation, Tars, Biomass, Kinetic modeling

Anisole was chosen as the simplest surrogate for primary tar from lignin pyrolysis to study the gas-phase chemistry of methoxyphenol conversion. Methoxyphenols are one of the main precursors of PAH and soot in biomass combustion and gasification. These reactions are of paramount importance for the atmospheric environment, to mitigate emissions from wood combustion, and for reducing tar formation during gasification. Anisole pyrolysis and stoichiometric oxidation were studied in a jet-stirred reactor (673-1173 K, residence time 2 s, 800 Torr (106.7 kPa), under dilute conditions) coupled with gas chromatography-flame ionization detector and mass spectrometry. Decomposition of anisole starts at 750 K and a conversion degree of 50% is obtained at about 850 K under both studied conditions. The main products of reaction vary with temperature and are phenol, methane, carbon monoxide, benzene, and hydrogen. A detailed kinetic model (303 species, 1922 reactions) based on a combustion model for light aromatic compounds has been extended to anisole. The model predicts the conversion of anisole and the formation of the main products well. The reaction flux analyses show that anisole decomposes mainly to phenoxy and methyl radicals in both pyrolysis and oxidation conditions. The decomposition of phenoxy radicals is the main source of cyclopentadienyl radicals, which are the main precursor of naphthalene and heavier PAH in these conditions.

Introduction

Environmental concerns such as the control of greenhouse gas emissions have led to an increased interest in the use of renewable energy. Biomass is widely used in combustion but can also be utilized in more advanced applications such as the production of a synthesis gas (syngas, a mixture of CO and H 2 ), which can be used for the production of liquid fuels such as Fisher-Tropsch or methanol. Lignocellulosic biomass may be a promising feedstock through the gasification processes [1] and [2], but tar is also formed during gasification [3]. The tar content in the product gas is the major cumbersome and problematic parameter in the gasification processes [START_REF] Milne | Biomass Gasifier "tars": Their Nature, Formation, and Conversion[END_REF]. Tar represents a complex mixture of over 100 compounds [START_REF] Evans | [END_REF] and [6]. It leads to fouling, coke deposition, and catalyst deactivation. Hence, tar conversion or removal is one of the main challenges for the successful development of commercial gasification technologies and has been extensively studied [7] and [8].

Evans and Milne [START_REF] Evans | [END_REF] defined three main classes of tars: primary tars (low temperature, oxygenated) and secondary and tertiary tars (benzene, polycyclic aromatic hydrocarbons-PAHs, etc.). In combustion and gasification reactors, the heaviest (tertiary) tars are formed by the conversion of primary tars, for instance, in the freeboard or emulsion phase of a fluidized bed or in the flame during combustion. A better understanding of the detailed gas-phase chemistry leading from primary tars to secondary and tertiary ones is needed. The aim of this work is to deeply understand these reactions under conditions representative of biomass gasification in air-blown fluidized beds, but also PAH formation in wood combustion furnaces [9]. A significant fraction of aromatic tars and soot originates from lignin pyrolysis products mainly composed of guaiacol-type and syringol-type units [START_REF] Evans | [END_REF], [10], [11] and [12]. Carbohydrates from cellulose and hemicellulose can also contribute to aromatic formation through the production of light unsaturated species such as acetylene and propargyl radicals [START_REF] Evans | [END_REF], [13] and [14]. The simplest aromatic model compound with a methyl ether unit, which is anisole, has been chosen as a surrogate for lignin primary tar, as a first step toward developing an advanced elementary kinetic model including major primary tar surrogates.

To date, there have been very few works on tar gas-phase conversion coupling well-controlled experiments, tar quantitative analysis, and elementary modeling, especially concerning anisole. Earlier experiments for anisole pyrolysis were performed in a stirred-flow reactor by Mulcahy et al. [15], who reported that reactions between methyl radicals and anisole proceed mainly by H-abstraction on the methoxy group, which is about 12 times faster than the attack on the aromatic ring. A mechanism for the formation of the main identified products (benzaldehyde, toluene, and methane) was also proposed. During the 1980s, Schlosberg et al. [16] quantified 10 products of anisole decomposition after a 2-min pyrolysis at 723 K. A shock tube study of the decomposition of anisole between 1000 and 1580 K and at pressures from 0.4 to 0.9 atm permitted an evaluation of the rate constant of the unimolecular decomposition of phenoxy radicals (C 6 H 5 O) to CO and C 5 H 5 and of the reaction of methyl radicals with phenoxy to form cresols [17]. Mackie et al. [18] studied the pyrolysis of anisole in a perfectly stirred reactor over the temperature range 850-1000 K and at low pressure (between 16 and 120 × 10 -3 atm). The yields of the main products (CH 4 , C 2 H 6 , C 6 H 5 OH, methylcyclopentadiene, CO, cresols, benzene, and benzaldehyde) as functions of temperature at two residence times (0.14 and 0.95 s) were reported and a rate constant was proposed for the unimolecular decomposition of anisole: Rate parameters were A = (2.9 ± 1.0) × 10 15 s -1 and Ea = 64.0 ± 0.6 kcal/mol. Arends et al. [19] studied the thermal decomposition of anisole in an excess of hydrogen in a flow reactor over the temperature range 793-1020 K and at atmospheric pressure. The main identified products were similar to those found in previous studies. A first kinetic model for anisole decomposition was proposed, which was composed of 37 reactions involving 23 species. However, some assumptions on the rate constant were made to obtain good agreement with experimental results. Pecullan et al. [20] presented the first study on anisole oxidation. Experiments near 1000 K were performed in the Princeton flow reactor and made it possible to observe that anisole undergoes first-order decay under both pyrolysis and oxidation conditions. It was proved also that under the studied conditions, the distribution of reaction products such as phenol, cresols, methane, and ethane was independent of the equivalence ratio. A new kinetic model was proposed, consisting of 66 reversible reactions involving 31 species. A good prediction of experimental data was reached for anisole, carbon monoxide, methylcyclopentadiene, and total phenolic compounds. More recently, Platonov et al. [21] showed the formation of various decomposition products, including PAHs, during the pyrolysis of anisole at high temperature (1023-1173 K). Friderichsen et al. [22] used a hyperthermal nozzle and a flow tube reactor followed by analysis by mass spectrometry and FTIR spectroscopy. Phenoxy and cyclopentadienyl radicals were found to be important intermediates in the decomposition of anisole. In addition, the identification of 9,10-dihydrofulvalene suggested that this molecule can be considered as a naphthalene precursor (2 C 5 H 5 → 9,10-dihydrofulvalene → naphthalene). To observe radical intermediates, Scheer et al. [23] and [24] used a hyperthermal tubular flow reactor with about 65 μs residence time. Their results confirm that the first step in the decomposition of anisole is the loss of the methyl group to form a phenoxy radical, followed by the ejection of CO to yield a cyclopentadienyl radical. They showed also that C5H5 decomposes to the propargyl radical (CH 2 CCH) and acetylene only at high temperature (>1373 K). Afterward, propargyl radical recombination is a significant benzene formation channel. At lower temperatures, most of benzene comes from a ring expansion reaction of methylcyclopentadiene (MCPD), resulting from methyl radicals combining with cyclopentadienyl radicals. Naphthalene formation was attributed also to cyclopentadienyl radical combination. The most recent paper dealing with anisole [25] proposed a new kinetic mechanism validated against experimental data from Pecullan et al. [20] and Platonov et al. [21]. The literature review shows that anisole is a fragile compound, which reacts easily by unimolecular O-CH 3 bond breaking, which leads to methyl and phenoxy radicals. This is consistent with the bond dissociation energies [START_REF] Luo | Handbook of Bond Dissociation Energies in Organic Compounds[END_REF] and [START_REF] Muller | [END_REF] involved in anisole and presented in Figure 1. To the best of our knowledge and despite all these extensive works, there is not yet a detailed kinetic mechanism for anisole pyrolysis and oxidation validated against well-controlled thermal conditions with quantitative analysis of products. Moreover, while the very low value of the C-O bond dissociation energy allows easy decomposition, most studies have been performed at temperatures higher than 1000 K, i.e., higher than that involved in biomass pyrolysis/gasification where primary tars are formed. The present work aims at investigating experimentally and theoretically the kinetics of anisole pyrolysis and oxidation at intermediate temperatures (673-1173 K) and to predict the formation of intermediates leading to secondary and tertiary tars, such as PAHs.

Experimental methods

The experimental study was performed in a jet-stirred reactor operated at constant temperature and pressure. The experimental apparatus has already been used for numerous pyrolysis and oxidation gas phase studies. Its description can be found in previous papers [28] and [29]. The main features of the experimental and analytical methods are described below. It has been shown that under the present conditions, this reactor can be considered as a zero-dimensional ideal perfectly stirred reactor with homogenous temperature and composition. The jet-stirred reactor (Figure 2) is made of fused silica and consists of an annular preheater and of a sphere in which the reaction takes place, with an injection cross located at the center. The reactor (volume 88 cm 3 ) and the injection cross were designed to obtain four turbulent jets, providing good mixing of the gas phase [30]. The heating of the preheater and of the reactor is achieved by the mean of Thermocoax resistances rolled around the different parts and controlled by Eurotherm 3216 controllers and type K thermocouples (temperature precision better than 1%). The heating of the preheater is divided into two zones, the first one heated at a temperature 150 K lower than the reaction temperature and the second one being heated at the same temperature as the reactor. The reaction temperature was measured by another type K thermocouple located in the intra-annular part of the preheater, with its extremity at the level of the injection cross in the center of the reactor. The feed of liquid anisole (Sigma Aldrich, purity 99.7%) was controlled by a Coriolis mass flow controller (Bronkhorst). This device was associated with an evaporator in which anisole was mixed with the carrier gas (helium IC or argon 5.0, Messer) and heated to 473 K. The carrier gas flow rate was controlled by a gas mass flow controller (Bronkhorst). The precision of the mass flow controllers is 0.5%. The carrier gas for experiments was helium. A set of experiments was also performed with argon to enable the quantification of hydrogen. The pressure in the reactor was kept constant at 800 Torr (106.7 kPa) using a valve located downstream of the reactor. The residence time of anisole was 2 s for both pyrolysis and oxidation studies.

Unconverted reactant and products were analyzed by means of several gas chromatographs (GC). For pyrolysis studies, the analysis was performed in two steps, as light and heavy species could not be analyzed at the same time. Heavy species (including more than five heavy atoms such as carbon and oxygen) were condensed in a trap connected to the reactor outlet and maintained at liquid nitrogen temperature. The content of the trap after addition of an internal standard (n-octane) and solvent (acetone) was analyzed in an Agilent 7890. A gas chromatograph fitted with an HP-5 capillary column, a flame ionization detector (FID), and a mass spectrometer for species identification. Light species (including less than six heavy atoms) were analyzed online using two gas chromatographs. The first one was fitted with a Plot-Q capillary column (providing good separation for hydrocarbons from methane to 1,3-cyclopentadiene) and a FID; the second one with a Carbosphere packed column (providing good separation for methane, C 2 hydrocarbons, and carbon oxides), a thermal conductivity detector (TCD), and a FID. Hydrogen was quantified using a Perichrome 1250 gas chromatograph fitted with a 5-Å molecular sieve and a TCD. The carrier gas in the JSR was switched from helium to argon for these analyses to enhance the sensitivity and to avoid the superposition of a peak of helium with the peak of H 2 . In the case of oxidation studies, the whole analysis was made online with a transfer line to GC heated to 413 K to avoid the condensation of heavy products. It has been checked that both methods (on-line and off-line) led to similar results. Note also that in the temperature range used in this study, the amount of radical species remains low in comparison with that of stable products; combination reactions cannot really affect the quantification. All the products were analyzed by gas chromatography using the instruments mentioned before. For heavy products, the calibration was made using a gas sample of known composition for anisole, benzene, toluene, phenol, and cresols. In the case of other heavy products, the carbon effective number method [START_REF] Tranchant | Manuel pratique de chromatographie en phase gazeuse[END_REF] was applied to quantify other species with anisole as a reference. The limit of detection was about 1 ppm for the heaviest species (FID), about 10 ppm for carbon oxides and oxygen (TCD), and 100 ppm for hydrogen. Uncertainty estimates were about ±5% for species that are analyzed online and ±10% for species condensed in the trap in pyrolysis experiments. In the case of phenolic compounds, i.e., phenol and cresols, their high affinity for surfaces led to higher uncertainties, around 30%.

Two series of experiments, in pyrolysis and oxidation respectively, were made in the range of temperature 673-1173 K. Measures were made every 25 K to follow the decay of the reactant mole fraction and also the evolution of the product mole fraction profiles. The inlet mole fraction of anisole was 0.01 in pyrolysis experiments, and was reduced to 0.005 in the oxidation study to avoid the fouling of the heated transfer line to the GCs. Oxidation experiments were performed under stoichiometric conditions ( = 1). In pyrolysis, the carbon mass balance was better than 95% up to 1100 K. Above this temperature, a progressive decrease in the carbon balance was observed, due to the formation of heavy PAHs and soots, which were not quantified. In oxidation, the mass balance remained correct over the whole temperature range, since carbon oxides are obtained at high temperature.

Kinetic modeling

A new detailed kinetic model for the pyrolysis and combustion of anisole has been developed based on a recent model of ethylbenzene oxidation [START_REF] Husson | [END_REF], which was used as a comprehensive reaction base for the reaction of aromatics and smaller species. This mechanism contains the different reaction pathways producing aromatic rings from lighter species, such as propargyl radical combinations and C4 and cyclopentadienyl reaction routes [33]. A submechanism for the reactions of anisole has been developed in the present work. Reactions of xylenes [34] and the reactions of formation of PAHs proposed by Slavinskaya et al. [35] have been included. Reactions with acetylene, involved in the socalled HACA mechanism, are included. Reactions of benzaldehyde, cresols, and naphthalene have been revised. The final model contains 303 species and 1922 reactions. Simulations were performed using the PSR (perfectly stirrer reactor model) module in the Chemkin software suite [36]. Thermochemical data were mainly evaluated by the software THERGAS [START_REF] Muller | [END_REF].

The primary mechanism for anisole decomposition and oxidation is composed of 24 elementary reactions (Table 1) including initiations, ipso-additions, H-abstractions, and radical decompositions. Reactions of initiation include the breaking of the two weakest bonds in the anisole molecule, i.e., O-CH 3 and OCH 2 -H bonds. Kinetic parameters of these two channels were estimated from the reverse radical combinations. The breaking of the C-O bond leading to phenoxy and CH 3 radicals (reaction 1 in Table 1) uses the same kinetic parameters for the reverse reaction as that used in the termination of benzyl radicals with CH3, forming ethylbenzene [37]. This corresponds to a value of the first-order decomposition of anisole at the upper limit in the uncertainty range of the rate constant proposed by Mackie et al. [18], with a factor of 2.2 at 1000 K. The breaking of the C-H bond (reaction 2) is modeled by kinetic parameters estimated by analogy with alkyl radicals [38]. Bimolecular initiations with oxygen molecules produce View the C 6 H 5 OCH 2  and HO 2 radicals (reaction 3). Kinetic parameters were estimated using rules proposed by Ingham and Walker for the alkanes [39], assuming an activation energy equal to the enthalpy of reaction. Ipso-additions involving H-atom, OH, and CH 3 radicals have been considered. For the reaction of anisole with H (reaction 4), the kinetic parameters of the equivalent reaction phenol + H were used [37]. Parameters for the reaction with OH (reaction 5) were estimated by analogy with the similar reaction of toluene + OH [40]. In the case of the reaction with CH3 (reaction 6), kinetic parameters were assumed to be that of the addition of a CH3 radical on a C=C bond [41]. Nine hydrogen atom abstractions on the methoxy group of anisole were also considered. Abstractions by H, OH, and HO 2 (reaction 8-10) were estimated using an Evans-Polanyi relationship proposed by Dean and Bozzelli for the reactions of alkanes [START_REF] Dean | Combustion chemistry of nitrogen[END_REF]. In the case of abstraction by H-atoms (reaction 8), this correlation leads to a high rate constant for the C-H bonds of weak bond dissociation energy [START_REF] Moss | [END_REF]. 1.0E12 0.0 0.0 Note: Rate constants are given in the form k = A T n exp(-E a /RT) where A has units of cm, mol, and s, T has units of K, and E a has units of cal/mol. The origin of the rate parameters is discussed in the text.

The rate constant was consequently decreased to slow down this reaction, which is very sensitive at the lowest studied temperatures: the A-factor was divided by 2 and the activation energy was increased by 1 kcal/mol, which is in the range of uncertainty of the correlation used. Kinetic parameters for H-abstractions by O atoms and C 2 H 5 radicals (reactions 7, 12) were estimated by analogy with the abstraction of a tertiary hydrogen in alkanes [44] because of the similarity between the bond dissociation energy of both C-H types. Rate parameters for the H-atom abstraction of anisole by the methyl radical CH 3 (reaction 11) were proposed in the literature by Mulcahy et al. [15]. Parameters for the reaction with the cyclopentadienyl radical (reaction 13) were estimated by analogy with the reaction of allyl with iso-butane [45], and that with the phenyl radical (reaction 14) by analogy with the reaction of phenyl with propane [46]. The reaction with phenoxy radical (reaction 15) was written in the reverse direction and its rate parameter was estimated from that of the reaction C 6 H 5 OH + CH 3 = C 6 H 5 O + CH 4 [47]. H-atom abstractions produce C 6 H 5 OCH 2  radicals, which can react by β-scission (reaction 16) or isomerize to C 6 H 5 OCH 2  by an internal ipso addition (reaction 17). This latter decomposes afterward, mainly to benzaldehyde (reactions [18][19][20]. Kinetic parameters are those theoretically calculated by da Silva and Bozzelli [48]. The combination of C 6 H 5 CH 2 O with the H-atom (reaction 21) and the disproportionations between C 6 H 5 CH 2 O and the Hatom, CH 3 , and C 5 H 5 were also added (reactions [22][23][24], with kinetic parameters estimated by analogy with the same types of reactions in alkane pyrolysis [38]. Some sensitive reactions have been also updated in the model of ethylbenzene oxidation [START_REF] Husson | [END_REF] used as a reaction base for the combustion of aromatic compounds. They are summarized in Table 2. Note that this reaction base was validated and used previously for combustion conditions, in which sensitive reactions are different. The modifications were made principally in phenol, 1,3-cyclopentadiene, and benzene submechanisms.

In the phenol submechanism, the high-pressure kinetic parameters for the reaction of combination of phenoxy radicals and hydrogen atoms (reaction 1 in Table 2) have been updated from He et al. [49]. The decomposition of the phenoxy radical to cyclopentadiene and carbon monoxide (reactions 2-3) has been described in more detail. Up to now, the phenoxy radical decomposed directly into C 5 H 5 and CO in the mechanisms proposed in the literature. According to theoretical calculations made in our laboratory, and including the pressure dependence [50], a first step has been added in which the phenoxy radical isomerizes to a bicyclic radical; this bicycle further forms a cyclopentadiene radical and carbon monoxide (reaction 3). The rate constant derived from theoretical calculations proposed by Xu and Lin [51] was used for the decomposition of phenol to C 5 H 6 and CO (reaction 4) [51]. Reactions of phenoxy radical with CH 3 , which especially produce cresols, were added on the basis of the theoretical calculations by Pecullan et al. (reactions 5-8) [20].

The 1,3-cyclopentadiene submechanism was found to be very sensitive in the fate of the phenolic intermediates. For reaction C 5 H 6 = C 5 H 5 + H (reaction 9), kinetics parameters from Robinson and Lindstedt [52] recommended at 1 atm were used. H-atom abstractions by H-atoms and CH 3 radicals were also reevaluated. The abstraction by H-atoms (reaction 10) was reevaluated theoretically at the CBS-QB3 level of theory [50] and the one involving CH 3 radicals was taken from Pecullan et al. [20] (reaction 11). The addition of H-atom on cyclopentadiene leading to allyl radical and C 2 H 2 (reaction 12) was added with the rate constant from Zhong and Bozzelli [53]. H-atom abstractions by C 3 H 3 (reaction 14, producing allene or propyne), C 4 H 3 (reaction 15, producing vinylacetylene), and C 5 H 7 (reaction 16, producing cyclopentene) were added and their rate constants were estimated by analogy with the equivalent H-atom abstractions on toluene [29]. The reaction rate of C 6 H 5 CO decomposition (reaction 17) was updated from Nam et al. [54] and the kinetic parameters for reaction C 6 H 6 + H = C 5 H 4 CH 3 (reaction 18) were taken from Marinov et al. [55]. The H-atom abstraction by C 6 H 5 on C 2 H 4 to form benzene and C 2 H 3 radical (reaction 19) was included based on the data from Tokmakov and Lin [56]. Thermochemical data for phenoxy radicals, cyclopentadienyl radicals, and benzyne were also updated, using data proposed by Sirjean et al. [50], Burcat's database [START_REF] Burcat | Ideal Gas Thermochemical Database with updates from Active Thermochemical Tables[END_REF], and JetSurf [START_REF] Sirjean | A High-Temperature Chemical Kinetic Model of n-Alkane Oxidation[END_REF], respectively. are given in the form k = A T n exp(-E a /RT) where A has units of cm, mol, and s, T has units of K, and E a has units of cal/mol. The origin of the rate parameters is discussed in the text.

Results

Experimental results

The pyrolysis of anisole has been investigated as a function of the temperature under the conditions described above. The evolution of the mole fraction of anisole is displayed in Figure 3. Anisole starts to decompose at 750 K and reaches 50% conversion at 873 K. Above 950 K, the conversion is complete and the fraction of anisole becomes close to zero. The reaction leads to the formation of dozens of compounds.

Light products include hydrogen, carbon monoxide, methane, ethane, ethylene, and acetylene, profiles of which are presented in Figure 4. These compounds reach high mole fractions and are mostly produced at high temperature. In the case of methane, a shoulder in the mole fraction profile can be seen between 900 and 950 K, which indicates a change in the mechanism. Methane is usually produced by H-atom abstraction by methyl radicals. This transition corresponds to the temperature of complete conversion of anisole. At higher temperature, the reactions of the products of anisole decomposition become dominant in the mixture; H-atom abstractions on these products may be more difficult than on anisole, and other reactions, such as addition on unsaturated bonds, may compete with the methane formation. Among the C 2 species, ethane and ethylene reach a maximum mole fraction about 1000 ppm, around 975 K and 1100 K, respectively. Acetylene is produced in large amounts above 1000 K. These observations are consistent with the decomposition channel of anisole to methyl and phenoxy radicals proposed in the literature [18], [19] and [20]. Methyl radicals can react indeed by H-atom abstraction on the reactant to produce methane and by combination to yield ethane and ethylene, leading to large amounts of these molecules. Phenoxy radicals can decompose to CO and cyclopentadienyl radicals. These latter radicals lead to the formation of 1,3-cyclopentadiene and of 5-methyl,1,3-cyclopentadiene by reactions with H-atoms and methyl radicals, respectively. Aromatic products are displayed in Figure 6. Benzaldehyde, phenol, cresols, and benzofuran, which contain an O-atom like anisole, have been quantified. Benzaldehyde is produced at the lowest temperatures, corresponding to the very beginning of the decomposition reaction of anisole, and its mole fraction reaches a maximum of about 250 ppm at 900 K. Phenol, cresols, and benzofuran are produced in larger amounts, around 5000, 1300, and 750 ppm, respectively. The main aromatic hydrocarbons are benzene, toluene, o-and p-xylene, ethylbenzene, styrene, naphthalene, and acenaphthylene. Benzene is the major one, with a peak around 2700 ppm. Naphthalene is produced at the highest studied temperatures (from 1000 K) and reaches 1000 ppm. The maximum mole fraction of toluene is around 300 ppm, while other products are minor. Note that heavier compounds such as dibenzofuran or phenanthrene were also detected in small amounts but could not be quantified with proper reproducibility. Some soot deposits were obtained at the exit of the reactor at the highest temperatures. Among the products, two groups can be distinguished: first, compounds with a bell-shaped profile such as phenol or cresols; second, products exhibiting a monotone profile such as CO, hydrogen, or acetylene. The selectivity of the reaction products varies as a function of the temperature. Below 950 K, the main products are phenol, methane, carbon monoxide, cresols, benzaldehyde, and benzofuran. When the temperature increases above 950 K, anisole is completely converted and primary products also reach a high degree of reaction. The major products are then carbon monoxide, hydrogen, methane, benzene, and polyaromatic species. Figure 7 presents the selectivity of anisole and oxygenated products among the oxygenated species as a function of temperature. Phenolic compounds are phenol, o-cresol, and p-cresol; furans are benzofuran and dibenzofuran. Phenolic compounds constitute the major oxygenated compounds below 1050 K. At higher temperatures, CO becomes the major oxygenate formed by the decomposition of the phenol and its derivatives, which leads also to aromatics, mainly benzene and naphthalene (Figure 6). Note that conversely to usual combustion conditions, light unsaturated soot precursors such as acetylene, allene, and propyne seem to be present in small amounts, which cannot yield noticeable amounts of aromatics. Furthermore, aromatic and polyaromatic hydrocarbon formation starts around 900 K, at temperatures much lower than those typical for such mechanisms [33]. Benzaldehyde exhibits noticeable selectivity only at the lowest temperature, while furans remain minor products over the whole temperature range. This evolution corresponds to the maturation of tars described by Evans and Milne [START_REF] Evans | [END_REF], going from primary tars, for which anisole is the simplest surrogate from lignin, to secondary tars, i.e., the phenolic compounds, and to tertiary tars, which are aromatics and polyaromatics produced by the decomposition of the phenolic molecules and the combination of light unsaturated hydrocarbons. The color of the solutions sampled in the condensers (Figure S1 in the Supplementary material) also illustrates this observation by varying significantly with the temperature from transparent (773-823 K) to yellow (873-1023 K) and brown (>1023 K), highlighting the formation of heavy species and even soots. The stoichiometric oxidation of anisole has been also investigated as a function of the temperature under the conditions described above. The evolution of the mole fraction of the reactants is displayed in Figure 8. The temperatures corresponding to the start of the reaction, the consumption of 50% of anisole, and the total conversion of the reactant are exactly the same as obtained in thermal decomposition, which is consistent with observations in the literature [20]. Anisole involves a weak bond and decomposes easily by a unimolecular reaction even at moderate temperatures. The addition of oxygen has very little influence under these conditions. Note that the oxygen consumption starts above 900 K, at a temperature for which anisole conversion is close to complete. In contrast to pyrolysis results, it appears that the addition of oxygen does not significantly affect the nature of the main products at low temperatures, but leads to the combustion products (CO and CO 2 ) at high temperatures, instead of producing aromatics and soot. The mole fraction profiles of the light products are displayed in Figure 9. The only new C 1 -C 2 products are carbon dioxide and acetaldehyde. CO and H 2 mole fractions reach a maximum around 1000 K and are completely converted to CO 2 and water at 1100 K. The amount of C 2 compounds is lower than in pyrolysis.

Ethylene remains the major one (800 ppm at the peak, close to the maximum amount of methane), followed by ethane (300 ppm) and acetylene (250 ppm). Acetaldehyde is abundant around 950 K (150 ppm). The nonaromatic heavier compounds are presented in Figure 10. Acrolein (75 ppm around 950 K) and 1-butene (12 ppm at 925 K) have been quantified in oxidation, while propane was not detected anymore. We can observe that the higher reactivity induced by the presence of oxygen in the mixture above 900 K leads to larger amounts of C 3 and C 4 compounds in comparison with pyrolysis. The peak of propene is now close to 45 ppm, while 1,3-butadiene and vinylacetylene maximum fractions are respectively three times and twice those obtained in pyrolysis. This can be explained by the greater decomposition of the aromatic intermediates allowed by their reactions with oxygenated radicals. Conversely, 1,3-cyclopentadiene and 5-methyl,1,3-cyclopentadiene amounts have dropped off. Figure 11 displays the profiles of aromatic products. The same products as in pyrolysis have been quantified, except m-guaiacol and hydroxybenzaldehyde, which were detected only in oxidation, and o-and p-xylenes, ethylbenzene, and acenaphthylene, which were not detected in oxidation. Phenol is the major aromatic product (around 1300 ppm at 900 K), while the maximum amounts of cresols and benzofuran decrease to 310 ppm and 150 ppm around 900-950 K, respectively. The amount of benzaldehyde, which is produced at low temperatures, is similar to that obtained in pyrolysis. Among the nonoxygenated aromatics, the most abundant product is benzene (370 ppm). In comparison with pyrolysis, the amount of naphthalene drops off dramatically to 20 ppm, which is consistent with the decrease in the amount of precursors such as C 5 H 6 . Traces of 2,3-dihydrobenzofuran and xanthone were identified but not quantified. Among the reaction products, the larger selectivity was obtained for carbon monoxide, phenol, carbon dioxide, methane, benzaldehyde, cresol, and hydrogen at 873 K. At 1023 K, anisole is mainly converted to carbon monoxide, phenol, and carbon dioxide. The maximum mole fraction of naphthalene is reached at about 1000 K and this compound is then completely converted under oxidative conditions, whereas its maximum was obtained at 1150 K under pyrolysis conditions without any further decrease in its profile.

Comparison between experimental data and simulations

The detailed model developed in this study aims at predicting the conversion of anisole and the formation of the main products on a large range of residence times and thermal conditions. The predictions using the model were compared to our experimental data (Figures 3,4, 5 and 6 in pyrolysis, Figures 8,9, 10 and 11 in oxidation). To extend the validity of the model, results obtained in the Princeton flow reactor by Pecullan et al. [20] were also simulated. In both thermal decomposition and oxidation, the model predicts the anisole reactivity perfectly (Figures 3 and8). The model simulates well the complete conversion of the reactant achieved at 950 K under both conditions and the profiles of the primary products. At higher temperatures, some discrepancies between simulations and experiments increase since the reaction is not that of anisole any more, but that of a very complex mixture of phenolic, aromatic, or furanic products. The amounts of most minor products are, however, estimated within a factor of 2. Note that the model uses as a reaction base a mechanism for light aromatics (benzene, toluene, ethylbenzene) validated under various condition of combustion [29], [START_REF] Husson | [END_REF] and [59]. Consequently, arbitrary adjustments of reactions in this submechanism have been avoided since they would yield little further information, but might change results for aromatic combustion drastically.

In pyrolysis, the behavior of the light products is well reproduced (Figure 4). Hydrogen and methane, which are produced by H-atom abstractions by H-atoms and methyl radicals, respectively, are correctly simulated at the lowest temperature, showing that the rate constants of the H-atom abstractions on anisole are correctly estimated. When the temperature increases, hydrogen is slightly overpredicted and methane underpredicted, reflecting the fact that H-atom abstractions on anisole become a minor channel in the formation of these light molecules above 830 K, conversely to H-atom abstraction on cresols, benzaldehyde, phenol, and 5-methyl,1,3-cyclopentadiene. Many of these latter rate constants are just estimates by means of structure-reactivity relationships. Ethylene is a little overestimated, while ethane is slightly underestimated, but the order of magnitude is correct. A major source of ethylene is the reaction between methyl radicals and HOC 6 H 4 CH 2 radicals that derive from cresol molecules, yielding ethylene and phenol. The product distribution and the rate constant of this reaction contained in the xylene submechanism [34] are estimated and uncertain. The amounts of C 3 and C 4 light products are simulated at least within a factor of 2 but the profile shapes fit the experimental data well (Figure 5). Vinylacetylene is the only one that is strongly underestimated. However, these intermediates are formed in small quantities, mostly above 900 K, and do not correspond to primary products of anisole decomposition but rather to products from other reactions, such as reactions of C 2 and the decomposition of aromatic rings. This induces larger uncertainties in their modeling. 1,3-Cyclopentadiene, which is produced by the decomposition of phenoxy radicals, and 5-methyl,1,3-cyclopentadiene, obtained from the combination of a methyl and a cyclopentadienyl radical, are also a little underestimated.

The simulation of aromatic products presented in Figure 6 is almost satisfactory. Considering the well-known experimental difficulties encountered in the quantification of phenols [18] and [START_REF] Pecullan | Pyrolysis and Oxidation Kinetics of Anisole and Phenol[END_REF], the oxygenated aromatics, which are formed during the early stage of anisole decomposition, are rather well reproduced, at least below 950 K. The mole fraction of phenol is, however, underpredicted, while that of cresols is overpredicted around 900 K. The formation of these species is competing through the fate of the phenoxy radical. We have used the values proposed by Pecullan et al. for the phenoxy + methyl radical multichannel reaction [20]. These rate constants theoretically evaluated are the only ones available in the literature, but lead to these slight discrepancies. The underprediction of benzofuran at low temperatures suggests, however, the lack of a route of formation for this compound. Among the cyclic hydrocarbons, benzene is correctly predicted at the lowest temperatures, but the model does not reproduce the experimental results above 900 K very well: its formation is first overpredicted and its amount decreases dramatically above 1000 K, unlike the experiments. Conversely, naphthalene, which is formed mainly from the combination of two cyclopentadienyl radicals, is overvalued by a factor of 3 above 1000 K. Benzene is mostly formed in our model by reactions of 5-methyl,1,3-cyclopentadiene. 5-Methyl,1,3-cyclopentadiene comes at the lowest temperatures from the reaction of phenoxy radicals with CH3 proposed by Pecullan et al. [20] and, when the temperature increases, from the combination of methyl and cyclopentadienyl radicals. The discrepancy in the simulation of benzene and naphthalene at high temperature is linked to the modeling of reactions of 1,3-cyclopentadiene and cyclopentadienyl radicals, which are especially tricky in pyrolysis in a complex mixture of unsaturated hydrocarbons and radicals. All the possible reactions are still not completely understood and are under discussion in the literature [50], [START_REF] Kislov | [END_REF], [62], [START_REF] Cavallotti | Proc. Combust. Inst[END_REF] and [START_REF] Scheer | [END_REF]. Especially, we have tried to describe more accurately the formation of naphthalene, including recent theoretical studies [62] and [START_REF] Cavallotti | Proc. Combust. Inst[END_REF]. However, the simulated results changed only marginally and the lack of information about the bimolecular reactions of the new intermediates added more uncertainties. Eventually, the reaction of two C 5 H 5 radicals yielding naphthalene and two H-atoms was preferred, with the rate constant proposed by Melius et al. [65]. Note, however, that the sum of the simulated mole fractions of benzene and naphthalene fits the experimental value well. Even if the formation of some PAHs has been taken into account in the model [35], the lack of soot formation may also contribute to an undervalued consumption of naphthalene.

Despite the small amount of xylenes (below 15 ppm), the model overpredicts the quantity of these species by a factor of 1.5 (o-xylene) and 3 (p-xylene). For ethylbenzene and styrene, the model overvalues the quantities a little compared to the experimental results. Beyond the compounds which have been quantified experimentally, the model predicts also some polyaromatic compounds such as pyrene, indene, and phenanthrene, thanks to the addition of the submechanism from Slavinskaya et al. [35]. Their simulated profiles are presented in Figure S2 in the Supplemental material. Pyrene, often used in models as the precursor for soot inception, reaches 400 ppm at 1100 K.

In the presence of oxygen, the conversion of the reactants is presented in Figure 8. The model prediction follows the experimental profiles perfectly. The major light products (Figures 9 and10), such as CO, CO 2 , and methane, are very well represented by the model. The other C 1 -C 2 species (Figure 9), which are produced in small quantities (20-800 ppm), are well reproduced, except in the case of acetaldehyde, an oxidation product of the small unsaturated species such as propene. This underprediction can be related to the underprediction of propene or 1,3-butadiene (Figure 10). As in pyrolysis, these latter molecules are produced by the decomposition of the aromatic rings, in which many uncertainties remain in the C 2 -C 4 product distribution, especially at the low temperatures investigated here. 1,3-Cyclopentadiene has a simulated profile slightly shifted toward the higher temperatures. Note that 1,3-cyclopentadiene is a very reactive species, which is in a quasi-equilibrium with resonance-stabilized cyclopentadienyl radicals. A small uncertainty in the production or consumption fluxes leads to a large variation of its mole fraction. At low anisole conversion, i.e., without the strong influence of secondary reactions, 1,3-cyclopentadiene is rather well reproduced (see Figure 10). Heavy products displayed in Figure 11 are very well simulated. The simulated profiles of benzaldehyde, benzofuran, benzene, toluene, and naphthalene are well estimated. This means that the most important channels of formation and consumption of these compounds are well handled by the mechanism. Only cresols are overestimated by a factor of 2. As written above, the rate constants theoretically determined by Pecullan et al. [20] for the reactions of phenoxy and methyl radicals lead to an overestimation of this channel. Hydroxybenzaldehyde is produced by the oxidation of the methyl group of cresols and is then slightly overpredicted at temperatures above 900 K. Among the phenolic species, m-guaiacol is the only one to be undervalued by the model, indicating that a formation pathway may be missing. m-Guaiacol comes in modeling from the combination of methyl and HOC 6 H 4 O radicals, this latter radical being produced by a H-atom abstraction and an oxidation reaction on the aromatic cycle of phenol. The discrepancies between experiments and simulations found in pyrolysis for hydrocarbons above 1000 K do not appear in oxidation. Since the complex molecular growth yielding soot precursors competes now with some fast oxidation reactions with O-atoms, OH radicals, or O 2 molecules, which oxidize aromatic molecules and radicals into smaller species. As stated above, the distribution of the products of these latter reactions is more uncertain, leading to a less accurate modeling of C 3 -C 4 molecules. The simulations also show a strong decrease in the xylene formation (Figure S3 in the Supplemental material), below 0.1 ppm. This explains why these species were not found experimentally anymore. The heavy polycyclic molecules are also present in oxidation; however, the predicted quantities of phenanthrene and pyrene are below 1 ppm.

In order to validate the primary model of decomposition of anisole, the experimental results obtained by Pecullan et al. [20] at atmospheric pressure near 1000 K have been also simulated. A short residence time of about 100 ms made it possible to form principally the primary products from anisole decomposition. Figure 12 displays experimental data and simulations in pyrolysis as a function of residence time for an inlet mole fraction of 1077 ppm of anisole and a temperature of 1003 K. A shift of 7 ms was done in the simulation in order to take into account the uncertainty of the zero point in the flow reactor. The model reproduces perfectly the anisole consumption and the formation of CO, 1,3-cyclopentadiene, and benzene. More differences are obtained with cresols and phenol, but these products are more difficult to quantify experimentally [START_REF] Pecullan | Pyrolysis and Oxidation Kinetics of Anisole and Phenol[END_REF], and the uncertainty in the rate constant of the reaction producing cresols was discussed above. Other validations against these data, especially in oxidation, are presented in the Supplementary material (Figure S5).

Analysis and discussion

The reaction pathways involved in the conversion of anisole are illustrated in Figure 13 at 800 K and Figure 14 at 1000 K, under the conditions of the jet-stirred experiments. At 800 K, the simulated conversion is 12% and anisole reacts mainly by the unimolecular bond breaking, yielding methyl and phenoxy radicals (42% of the total anisole consumption). At this low temperature, the decomposition of the resonance-stabilized phenoxy radicals to CO and cyclopentadienyl radicals does not occur. The phenoxy radicals react rather by H-atom abstraction on anisole or primary products to produce phenol or with CH 3 to produce cresols and to a lesser extent CO and 5-methyl,1,3-cyclopentadiene. Phenol and cresols come almost completely from phenoxy radicals. Anisole reacts also by H-atom abstraction by H-atom, methyl, phenyl, and phenoxy radicals (39% of the total anisole consumption). These reactions yield C 6 H 5 OCH 2 radicals, which isomerize entirely to C 6 H 5 CH 2 O. This latter decomposes mainly to benzaldehyde. Minor routes of consumption of anisole are ipso-additions of H-atoms yielding benzene (14%) and of methyl radicals yielding toluene (5%). When the temperature increases, the bimolecular reactions of anisole does not compete any more with the bond breaking and becomes negligible. At 1000 K, the conversion is almost complete and 99.9% of anisole is converted into phenoxy and methyl radicals. Phenoxy radicals react still with methyl radicals to produce cresols and by H-atom abstraction or by combination with H-atoms to produce phenol, but another important reaction channel is now the decomposition to CO and C 5 H 5 radical through a bicyclic intermediate. Sensitivity analyses performed in pyrolysis at 800 and 1000 K confirm the relative simplicity of the initial channels of anisole reaction (Figure 15; coefficients at 1000 K have been divided by 10). At 800 K, the promoting reactions for anisole conversion are the unimolecular initiation and the metatheses with H-atom, methyl, and phenoxy radicals. The only noticeable inhibiting termination is the combination of phenoxy radicals with H-atoms, producing phenol. At high temperatures, the unimolecular initiation becomes by far the most sensitive step. Since the radical bimolecular reactions are then negligible in anisole consumption, no inhibiting process really appears. Very similar results were obtained in oxidation, which explain the similar reactivity of anisole in pyrolysis and oxidation whatever the equivalence ratio [20]. One of the main concerns in biomass combustion or gasification processes is the formation of tertiary tars, PAHs, and ultimately soot. During the decomposition of the lignin, large amounts of methoxyphenols are released. Anisole is a simplified surrogate of many methoxy-phenols formed during lignin pyrolysis [66]. More complex methoxyphenols such as guaiacol involve the same methoxy group and should react in the same way. Anisole decomposes to phenoxy radicals, which leads to large amounts of cyclopentadienyl radicals. At 1000 K, C 5 H 5 produces mainly naphthalene in pyrolysis. In oxidation, this PAH is still formed by the combination of two C 5 H 5 radicals, but its formation is in competition with the oxidation of the radicals, which reaches 75% of C 5 H 5 consumption, as shown in the flow analysis of cyclopentadienyl radicals presented in the Supplemental material (Figure S8). Figure 16 displays the sensitive reactions to the formation of naphthalene at 1000 K. In pyrolysis, the most promoting reaction is the decomposition of the bicyclic intermediate (bicycle in Figure 16) formed in the phenoxy radicals to CO + C 5 H 5 pathway. This reaction involves a high activation energy and is the primary source of C 5 H 5 radicals. The formation of the bicycle from phenoxy radicals is then also promoted. The combination of two methyl radicals to form ethane also promotes naphthalene formation by removing CH 3 from the mixture, since this latter radical can also combine with C 5 H 5 to yield 5-methyl,1,3-cyclopentadiene. The inhibiting reactions are mainly the combinations of phenoxy radicals with H-atoms and methyl radicals, which compete with the decomposition to C 5 H 5 , and that of C 5 H 5 and methyl radicals. In oxidation, the limiting step in the formation of naphthalene is no longer the formation of C 5 H 5 radical from phenoxy radicals, but the combination of C 5 H 5 , which competes with its oxidation by some reactions with O 2 and HO 2 . The reactions of phenol with O-atoms and OH radicals to produce hydroxyphenoxy and hydroxyphenyl radicals also inhibit the formation of naphthalene by competition with the reactions of formation of the phenoxy radical, a precursor of C 5 H 5 .

Figure 16: Sensitive reactions in the production of naphthalene at 1000 K during pyrolysis and oxidation of anisole.

Conclusions

In biomass combustion and gasification reactors, the first steps of reaction are a thermal decomposition, which generates hundreds of primary tars from cellulose, hemicelluloses, and lignin. These species are further converted in the gas phase in what is called tar maturation, i.e., their evolution toward species with lower O/C and H/C ratios. Lignin pyrolysis produces methoxyphenols, including heavy oligomers. Using anisole as a surrogate for lignin primary tars, its thermal decomposition and its oxidation under stoichiometric conditions at moderate temperatures have been studied. It was found that the decomposition of such a fragile molecule occurs mainly by the breaking of the methoxy O-C bond, producing large amounts of phenoxy radicals. This leads to some observations.

• The decomposition of anisole, and thereafter that of methoxyphenols, is a first-order reaction that occurs easily even at moderate temperature, without any influence of the addition of oxygen.

• The main primary products are phenol and phenolic molecules (secondary tars).

• Further reactions of the primary products lead mainly to cyclopentadienyl radicals and 1,3-cyclopentadiene, which are well-known PAH and soot precursors.

• An early formation of first aromatic rings and of naphthalene occurs, which does not involve the usual mechanisms of PAH formation in combustion based on C 2 or C 3 species [67].

The HACA mechanism looks like a minor route for PAH formation below 1000 K in biomass gasification or mild combustion reactors. Indeed, PAHs are already formed from 900 K under pyrolysis [68] or oxidation [69] conditions and CH 4 and C 2 compounds are poorly converted under slightly oxidizing conditions [70]. PAHs in biomass gasification or mild combustion reactors would be mainly formed (from 900 to 1100 K) from phenoxy radicals, producing cyclopentadienyl radicals and further PAHs and soot. This finding, derived from this study on anisole, needs to be confirmed for a broader range of temperatures, with other model compounds (such as furans from carbohydrates), and on real biomass gas.

C 6 H 5

 65 OCH 3 C 6 H 5 O + CH 3

Figure 1 :

 1 Figure 1: Structure of anisole and bond dissociation energies.

Figure 2 :

 2 Figure 2: Jet-stirred reactor (JSR) used to study the gas-phase conversion of anisole.

Figure 3 :

 3 Figure 3: Thermal decomposition of anisole. Points refer to experiments and lines to modeling.

Figure 4 :

 4 Figure 4: Evolution of the mole fraction of light products of anisole pyrolysis at residence time 2 s and pressure 800 Torr. Points refer to experiments and lines to modeling.

Figure 5 :

 5 Figure 5: Evolution of the mole fraction of nonaromatic C 3 -C 6 products of anisole pyrolysis at residence time 2 s and pressure 800 Torr. Points refer to experiments and lines to modeling.

Figure 6 :

 6 Figure 6: Evolution of the mole fraction of aromatic products of anisole pyrolysis at residence time 2 s and pressure 800 Torr. Points refer to experiments and lines to modeling.

Figure 7 :

 7 Figure 7: Selectivity of selected compounds among the oxygenated species (mol%) as a function of temperature during the pyrolysis of anisole.

Figure 8 :

 8 Figure 8: Evolution of the mole fraction of the reactants during the stoichiometric oxidation of anisole at residence time 2 s and pressure 800 Torr. Points refer to experiments and lines to modeling.

Figure 9 :

 9 Figure 9: Evolution of the mole fraction of light products during the stoichiometric oxidation of anisole at residence time 2 s and pressure 800 Torr. Points refer to experiments and lines to modeling.

Figure 10 :

 10 Figure 10: Evolution of the mole fraction of nonaromatic C 3 -C 6 products during the stoichiometric oxidation of anisole at residence time 2 s and pressure 800 Torr. Points refer to experiments and lines to modeling.

Figure 11 :

 11 Figure 11: Evolution of the mole fraction of aromatic products during the stoichiometric oxidation of anisole at residence time 2 s and pressure 800 Torr. Points refer to experiments and lines to modeling.

Figure 12 :

 12 Figure 12: Mole fraction profiles of reactant and selected products during the pyrolysis of anisole in the Princeton flow reactor [18]. X anisole = 1077 ppm, T = 1003 K, P = 1 atm. Points refer to experiments and lines to the present modeling.

  C 5 H 5 recombines mostly with itself to form naphthalene, with CH 3 to produce 5-methyl,1,3-cyclopentadiene, and with H-atoms leading to 1,3-cyclopentadiene. The reaction of phenoxy radicals with o-cresol to form benzofuran is a very minor channel.

Figure 13 :

 13 Figure 13: Reaction flux analysis for anisole pyrolysis at 800 K, 800 Torr, residence time of 2 s corresponding to 12% anisole conversion. Reaction fluxes (numbers associated with arrows) are relative to the consumption of a given species.

Figure 14 :

 14 Figure 14: Reaction flux analysis for anisole pyrolysis. Conditions: 1000 K, 800 Torr, residence time of 2 s corresponding to 99.9% anisole conversion. Reaction fluxes (numbers associated with arrows) are relative to the consumption of a given species.

Figure 15 :

 15 Figure 15: Sensitive reactions in the consumption of anisole in pyrolysis at 800 and 1000 K.Coefficients at 1000 K have been divided by 10 for presentation.

Table 1 :

 1 Primary mechanism for anisole pyrolysis and oxidation.

	Reaction	A	n	E a
	Initiations
	(1) C 6 H 5 O + CH 3 = anisole	5.0E12	0	0
	(2) C 6 H 5 OCH 2 + H = anisole	1.0E14	0	0
	Bimolecular initiation
	(3) Anisole + O 2 = C 6 H 5 OCH 2 + OOH	2.1E13	0	47,900
	Ipso-additions
	(4) Anisole + H = C 6 H 6 + CH 3 O	2.23E13	0	7930
	(5) Anisole + OH = C 6 H 5 OH + CH 3 O	7.8E2	2.88	3220
	(6) Anisole + CH 3 = toluene + CH 3 O	9.64E10	0.0	8000
	H-abstractions
	(7) O + anisole = OH + C 6 H 5 OCH 2	3.0E13	0.0	3280
	(8) Anisole + H = H 2 + C 6 H 5 OCH 2	3.6E8	1.5	5730
	(9) Anisole + OH = H 2 O + C 6 H 5 OCH 2	3.6E6	2.0	-850
	(10) Anisole + OOH = H 2 O 2 + C 6 H 5 OCH 2	4.2E3	2.69	16,590
	(11) Anisole + CH 3 = CH 4 + C 6 H 5 OCH 2	5.01E11	0.0	10,500
	(12) Anisole + C 2 H 5 = C 2 H 6 + C 6 H 5 OCH 2	3.0E11	0.0	9200
	(13) Anisole + C 5 H 5 = C 5 H 6 + C 6 H 5 OCH 2	5.40	3.3	17,200
	(14) Anisole + C 6 H 5 = C 6 H 6 + C 6 H 5 OCH 2	3E11	0.0	3850
	(15) C 6 H 5 OH + C 6 H 5 OCH 2 = anisole + C 6 H 5 O	1.8E11	0.0	7700
	Radical decompositions
	(16) C 6 H 5 OCH 2 = HCHO + C 6 H 5	2.61E14	0.05	39,280
	(17) C 6 H 5 OCH 2 = C 6 H 5 CH 2 O	8.78E10	0.44	19,030
	(18) C 6 H 5 CH 2 O = C 6 H 5 CHO + H	5.26E28	-5.08	22,250
	(19) C 6 H 5 CH 2 O = C 6 H 6 + CHO	2.37E32	-6.09	28,810
	(20) C 6 H 5 CH 2 O = C 6 H 5 + HCHO	7.21E33	-6.21	36,850
	(21) C 6 H 5 CH 2 O + H = C 6 H 5 CH 2 OH	1.0E14	0.0	0.0
	(22) C 6 H 5 CH 2 O + CH 3 = C 6 H 5 CHO + CH 4	1.0E12	0.0	0.0
	(23) C 6 H 5 CH 2 O + H = C 6 H 5 CHO + H 2	1.0E13	0.0	0.0
	(24) C 6 H 5 CH 2 O + C 6 H 5 = C 6 H 5 CHO + C 6 H 6

Table 2 :

 2 Modifications in the mechanism of ethylbenzene oxidation[START_REF] Husson | [END_REF].

	Reaction	A	n	E a
	Phenol submechanism
	(1) C 6 H 5 O + H(+M) = C 6 H 5 OH(+M)	2.00E14	0	0
	(2) C 6 H 5 O = bicycle	6.17E78	-18.794	90,500
	(3) bicycle = C 5 H 5 + CO	8.29E20	-2.422	11,110
	(4) C 6 H 5 OH = C 5 H 6 + CO	2.67E14	-0.61	74,100
	(5) C 6 H 5 O + CH 3 = MCPD + CO	2.05E+75	-18.29	38,900
	(6) C 6 H 5 O + CH 3 = OC 6 H 4 CH 3 + H	3.33E+39	-7.78	31,700
	(7) C 6 H 5 O + CH 3 = HOC 6 H 4 CH 3	2.31E+73	-17.37	38,800
	(8) C 6 H 5 O + CH 3 = C 6 H 4 OH + CH 3	1.29E-31	13.2	15,600
	1,3-Cyclopentadiene submechanism
	(9) C 5 H 6 = C 5 H 5 + H	3.23E49	-10.01	101,130
	(10) C 5 H 6 + H = C 5 H 5 + H 2	6.67E07	1.874	4259
	(11) C 5 H 6 + CH 3 = C 5 H 5 + CH 4	3.11E11	0	5500
	(12) C 5 H 6 + H = C 3 H 5 + C 2 H 2	7.74E36	-6.81	32,920
	(13) C 5 H 6 + C 3 H 3 = C 5 H 5 + aC 3 H 4	0.8E12	0	15,100
	(14) C 5 H 6 + C 3 H 3 = C 5 H 5 + pC 3 H 4	0.8E12	0	15,100
	(15) C 5 H 6 + C 4 H 3 = C 5 H 5 + C 4 H 4	1.6E12	0	15,100
	(16) C 5 H 6 + C 5 H 7 = C 5 H 5 + C 5 H 8	1.6E12	0	15,100
	Benzene submechanism
	(17) C 6 H 5 CO = C 6 H 5 + CO	5.27E14	0	29,010
	(18) C 6 H 6 + H = C 5 H 5 CH 3	2.39E27	-3.92	29,200
	(19) C 6 H 5 + C 2 H 4 = C 6 H 6 + C 2 H 3	9.45E-3	4.47	4470
	Note: Rate constants

Acknowledgments

This work was partially supported by the Energy program of CNRS through the project "CRAKIN" and the French ANR through the project "GAMECO" (coordinated by EDF).


Sonia Cafieri
email: sonia.cafieri@enac.fr

Alberto Costa
email: costa@sutd.edu.sg

Pierre Hansen
email: pierre.hansen@gerad.ca

Adding Cohesion Constraints to Models for Modularity Maximization in Networks

Finding communities in complex networks is a topic of much current research and has applications in many domains. On the one hand, criteria for doing so have been proposed, the most studied of which is modularity. On the other hand, properties to be satisfied by each community of a partition have been suggested. It has recently been observed that one of the best known such properties, i.e., Radicchi et al.'s weak condition [

Introduction

Networks are often used for representation and study of complex systems with applications in many domains. A network, or graph, G = (V, E) consists of a set of vertices V and a set of edges E. Vertices are represented by points and associated with entities of the system under study. Edges are represented by lines joining pairs of vertices and are associated with relations between the entities. The shape of these lines does not matter, but only their presence or absence. Newman recently gave a detailed introduction to networks in [START_REF] Newman | Networks: An introduction[END_REF].

In many complex systems there are sets of entities which share some common characteristics and/or are likely to have some common function. Henceforth, they will be called communities or clusters or modules. In the associated network, these communities correspond to sets of vertices for which the number of inner edges, that is edges joining two vertices of the community, is larger than the (possibly weighted) number of outer edges, that is the number of edges joining two vertices one of which belongs to the community and the other not. Detection of communities is a currently central and much studied problem in the theory and application of network science, with a vast literature, see [START_REF] Fortunato | [END_REF] for an in-depth survey and references therein. Usually, one seeks a partition of the given set of entities into disjoint communities, that is each entity must belong to one and only one community. Sometimes this last condition is relaxed: the aim is to find a covering of the set of entities, that is some entities may belong to several communities. We do not consider overlapping communities in this paper. The quality of the partition obtained can be judged in several ways:

• Some heuristics do not involve a criterion to be optimized. An example is Girvan and Newman's [3] edge removal heuristic in which edges with maximum betweeness are iteratively removed, yielding partitions into an increasing number of communities. Then the quality of the results obtained can only be judged a posteriori, usually based upon some substantive information. • A variety of criteria to be minimized or maximized have been suggested by several authors. The best known of them is modularity, introduced by Girvan and Newman [3]. Modularity of a community is defined as the difference between the number of inner edges and the expected number of edges in a random configuration model which keeps the distribution of vertex degrees unchanged. A large number of heuristics [4][5][6][7][8][9][10][11][12][13] provide in moderate computing time a near optimal solution, and a few exact algorithms provide an optimal solution [14][15][16]. Strengths and weaknesses of modularity are discussed in [START_REF] Fortunato | [END_REF][17][18][19]. Another, more recent, criterion is modularity density [20]. • Instead of considering an objective function, one may specify conditions to be satisfied by each community of a partition. The first two such conditions, a strong and a weak one, were proposed by Radicchi et al. in [21]. Two further conditions, a semi-strong and a extra-weak one, were suggested by Hu et al., in [22]. The almost-strong condition was introduced in [23]. Precise definitions and some properties of these five conditions will be given in the next section. We call these five conditions cohesion conditions.

Several authors have checked if some of these conditions are violated by some community of an optimal modularity maximizing partition [24,25]. Moreover, in [START_REF] Wang | Lecture Notes in Operation Research, The Third International Symposium on Optimization and Systems Biology (OSB09)[END_REF] it is proposed to add the weak constraint to models for modularity maximization and for modularity density maximization.

The aims of the present paper are: (i) to study to what extent optimal partitions for modularity maximization of real world problems satisfy the five conditions mentioned above and (ii) to examine the effect of imposing these constraints one at a time in modularity maximization models. These models are solved exactly. Observe that using exact optimization methods does allow separation of the effect of imposing a constraint from the possible error due to the use of a heuristic (other reasons which justify using exact optimization together with heuristics for their mutual improvement are given in the introduction of [16]).

The paper is organized as follows. Definitions of five cohesion conditions are given in Sect. 2 and some of their properties are discussed. It is then examined in Sect. 3 to what extent the optimal modularity maximizing partitions of 11 real-world problems from the literature do satisfy these conditions. Modifications to models for modularity maximization due to adding each of the five cohesion constraints are studied in Sect. 4. Computational results are presented in Sect. 5 which also includes a detailed discussion of the effect of imposing the weak condition on the optimal partitions of three well known problems. Conclusions are drawn in Sect. 6.

Cohesion Conditions

As mentioned above, an important approach to communities detection in network is based on the satisfaction of reasonable a priori conditions to have a community. Radicchi et al. [21] proposed two such conditions defining communities in a strong and a weak sense, respectively. Hu et al. [22] propose two further conditions defining communities in (what we call) a semi-strong and in an extra-weak sense, respectively. A further condition, definining communities in an almost-strong sense, was introduced in [23]. In this section, we give mathematical definitions of these five conditions and discuss some of the properties of corresponding communities.

Let G = (V, E) be a network with vertex set V and edge set E. Recall that the degree k i of a vertex v i belonging to V is the number of its neighbors (or adjacent vertices). Let S ⊆ V be a subset of vertices. Then the degree k i can be separated into two components k in i (S) and k out i (S), that is the number of neighbors of v i inside S and the number of neighbors of v i outside S. Let M be the number of communities of a partition S 1 , S 2 , . . . , S M of V and let A = (A ij ) be the adjacency matrix of G, where A ij = 1 if an edge joins vertices v i and v j and A ij = 0 otherwise.

• Strong Cohesion Condition (SCC): [21]

A set of vertices S forms a community in the strong sense if and only if every one of its vertices has more neighbors within the community than outside:

∀v i ∈ S k in i (S) > k out i (S). (1)
• Almost-Strong Cohesion Condition (ASCC): [23] A set of vertices S forms a community in the almost-strong sense if and only if every one of its vertices with degree different from 2 has more neighbors within the community than outside, and every vertex having degree 2 has at least one neighbor in the same community:

∀v i ∈ S | k i = 2 k in i (S) > k out i (S) (2)
∀v i ∈ S | k i = 2 k in i (S) > 0. (3)
• Semi-Strong Cohesion Condition (SSCC): [22] A set of vertices S forms a community in the semi-strong sense if and only if every one of its vertices has more neighbors within the community than the maximum number of neighbors within any other community:

∀v i ∈ S k in i (S) > max t=1,2,...,M, S =St vj ∈St A ij . (4)
Note that we use a strict inequality in the last formula instead of a non-strict one, as proposed by Hu et al. [22]. This is done for two reasons: first, to express all five conditions in an uniform way and, second, because in this form the strong cohesion condition implies the semi-strong one. It is not the case otherwise as shown on Fig. 1(a). Indeed, the partition into two communities C 1 = {1, 2, 3} and C 2 = {4, 5, 6} satisfies the semi-strong condition with equality but not the strong condition.

• Weak Cohesion Condition (WCC): [21]

A set of vertices S forms a community in the weak sense if and only if the sum of internal degrees within S is larger than the sum of external degrees, that is the number of edges joining S to the rest of the network V \ S:

vi∈S k in i (S) > vi∈S k out i (S). (5)
This is equivalent to the condition that the number of edges within S is at least half the number of edges in the cut of S. As already observed in [21] the strong cohesion condition implies the weak cohesion condition, but the converse is not true.

• Extra-Weak Cohesion Condition (EWCC): [22] A set of vertices S forms a community in the extra-weak sense if and only if the sum of internal degrees within S is larger than the maximum number of edges joining a vertex of S to a vertex in some other community in the rest of the network:

vi∈S k in i (S) > max t=1,2,...,M, S =St vi∈S vj ∈St A ij . (6)
Again, we use a strict inequality in the last formula instead of a non-strict one as suggested in [22].

Let us now consider implications between cohesion conditions. Clearly, the strong cohesion condition implies the weak cohesion condition as inequality ( 5) is equal to the sum of the inequalities (1) for all vertices v i ∈ S. The strong cohesion condition implies the semi-strong cohesion condition as they have the same left-hand side and max t=1,2,...,M vj ∈St A ij ≤ k out i (S) ∀v i ∈ S. Note that both conditions are identical when the number of communities M is equal to 2. Similarly, the weak cohesion condition does imply the extra-weak condition as left-hand sides are the same and max t=1,2,...,M, S =St vi∈S vj ∈St A ij ≤ vi∈S k out i (S). Again, the right-hand sides are equal when M = 2. The strong conditions imply the almost-strong as the latter is equivalent to the former, except for vertices of degree 2 for which the almost-strong condition is simply a weakened version of the strong (the strict inequality being replaced by a non-strict one). Finally, the semistrong cohesion condition implies the extra-weak cohesion condition as summing the left-hand side of ( 4) gives the left-hand side of ( 6) and summing the right hand-side of (4) gives vi∈S max t=1,2,...,M, S =St vj ∈St A ij = max t=1,2,...,M, S =St vi∈S vj ∈St A ij , that is the right hand side of ( 6). These implications are represented as follows:

ASCC ⇐ SCC ⇒ SSCC ⇓ ⇓ W CC ⇒ EW CC
These implications are asymmetric, i.e., the reverse implication does not hold. Consider for example the graph on Fig. 1(a). This graph with 6 vertices and 9 edges admits a partition into two communities C 1 = {1, 2, 3} and C 2 = {4, 5, 6} satisfying the conditions WCC and EWCC, but not the conditions SCC and SSCC. Moreover, no other partition in two communities satisfies one or the other of these conditions. So, the implications SCC ⇒ W CC and SSCC ⇒ EW CC are asymmetric. Let us next consider the graph of Fig. 1(b). This graph with 9 vertices and 14 edges admits a partition into three communities C 1 = {1, 2}, C 2 = {3, 4, 5} and C 3 = {6, 7, 8, 9} satisfying the condition EWCC, but not the condition WCC.

Cohesion Conditions in modularity maximization

We next check if the optimal solutions obtained by modularity maximization for a set of real-world problems do satisfy or not, and to which degree, the five cohesion conditions described above. To that effect, we use

G-2014-11

Les Cahiers du GERAD the optimal solutions of 11 test problems whose maximum modularity partitions have been obtained exactly using the methods presented in [14][15][16]. The five cohesion conditions have been coded and checked for each cluster of each such partition. Results are summarized in Table 1, which shows that:

• The strong cohesion condition (SCC) appears to be difficult to satisfy: while it holds for at least one of the communities of each optimal partition for all datasets but one, it is satisfied by all communities of one partition only. On average, the strong condition holds for 37.93% of the communities taken together. • The almost-strong cohesion condition (ASCC) is also difficult to satisfy. It is satisfied by zero to three more clusters than the strong condition. On average, this condition is satisfied by 53.45% of the communities. • The semi-strong cohesion condition (SSCC) again appears to be difficult to satisfy. Indeed, once more for one only of the problems this condition holds for all communities in the optimal partition. The average number of communities which satisfy the semi-strong condition is 51.72% of them. • The weak cohesion condition (WCC) appears to be much easier to satisfy than the former three. Indeed, it is satisfied by all communities of 8 out of 11 of the optimal partitions. When this is not the case, that is for dolphins, p53 protein, and political books, only one community does not satisfy the weak cohesion condition. The average number of communities which satisfy this condition is 94.83% of them. The difference between the maximum modularity partitions obtained without and with the weak cohesion condition will be further studied in Sect. 5. • The extra-weak cohesion condition (EWCC) appears to be easy to satisfy. Indeed, it it verified by all but one of the optimal partitions, that is by political books, and in that case it does not hold for one community only. The average number of communities which satisfy the extra-weak condition is 98.28% of them.

In the next section we proceed to add the five types of cohesion constraints, one at a time, to existing models for modularity maximization.

Table 1: Networks from real-world problems: network dimension (n = number of vertices, m = number of edges), number of clusters (M ) found by an exact algorithm for modularity maximization and number of clusters verifying the strong condition (M strong), the almost-strong condition (M almost-strong), the semistrong condition (M semi-strong), the weak condition (M weak) and the extra-weak condition (M extraweak). Percentage of communities satisfying each of these conditions are given in the last line.

Adding cohesion constraints to modularity maximization models

Exact modularity maximization without further constraints has been studied in [14][15][16]. Exact methods are applied using two different formulations. The first one is a reformulation of the problem as a clique partitioning problem [37], the second one is a mixed-integer quadratic optimization problem [14,16].

Adding cohesion constraints is more or less difficult according to the approach selected. In clique partitioning assignment of entities to communities is not explicitly considered: it only appears as a consequence of the optimal solution, so adding cohesion constraints does not seem to be easy in that framework. Contrarywise, the quadratic mixed-integer formulation already uses variables to denote assignment of vertices or of edges to communities. So, in that case cohesion constraints can be added easily.

We now describe these new constraints in more detail. Before that we recall the main elements in the model for modularity maximization of [14]. The objective function is the modularity Q of the network under study, that is the sum of modularities of its communities:

Q = M s=1 m s m - D s 2m 2 , (7)
where m s denotes the number of edges in community s, that is the subgraph induced by S s , D s denotes the sum of degrees k i of the vertices of community s, and M is the number of communities which is not a priori known. Binary variables are then used to identify the communities to which each vertex and each edge belongs:

X rs = 1 if edge r belongs to community s 0 otherwise for r ∈ E and s = 1, 2, . . . M and

Y is = 1 if vertex v i belongs to community s 0 otherwise for v i ∈ V .
Constraints involving X and Y variables express that an edge belongs to a community only if both its end vertices belong to that community:

∀r = {v i , v j } ∈ E, ∀s ∈ {1, . . . , M } X rs ≤ Y is ∀r = {v i , v j } ∈ E, ∀s ∈ {1, . . . , M } X rs ≤ Y js .
In fact, due to maximization, this edge will belong to this community if the two conditions given above are satisfied (hence it is not necessary to add the constraint X rs ≥ Y is + Y js -1 and X rs ≥ 0). Further constraints are presented in [14] and allow bounding the cardinality of each community from below or above, using binary variables F s such that F s is equal to 1 if the community s is nonempty, 0 otherwise. Finally, symmetry breaking constraints are used, as:

vi∈V |i≤M, s=1,...,i Y is = 1, (8)
and, following Plastria [38]:

∀s ∈ {3, . . . , M -1} ∀i ∈ {s, . . . , n} j=2,...,i-1 ℓ∈{1,...,s-1} Y jℓ - ℓ∈{1,...,s} Y iℓ ≤ i -3. ( 9
)
The resulting model, derived from that presented in [14], is the following:

max M s=1 m s m - D s 2m 2 (10) s.t. ∀r = {v i , v j } ∈ E, ∀s ∈ {1, . . . , M } X rs ≤ Y is (11) ∀r = {v i , v j } ∈ E, ∀s ∈ {1, . . . , M } X rs ≤ Y js (12) ∀s ∈ {1, . . . , M } m s = r∈E X rs ( 13
)
∀v i ∈ V M s=1 Y is = 1 (14) ∀s ∈ {1, . . . , M } D s = vi∈V k i Y is (15) ∀s ∈ {2, . . . , M } F s ≤ F s-1 (16
Y iℓ ≤ i -3 ( 20
)
∀v i ∈ V, ∀s ∈ {1, . . . , M } Y is ∈ {0, 1} (21)
∀s ∈ {1, . . . , M } F s ∈ {0, 1}

∀r ∈ E, ∀s ∈ {1, . . . , M } X rs ∈ R.

To this basic model we add cohesion constraints as follows. Notice that, for each cohesion condition, the constant 1 in the right-hand side transforms the strict inequality into a non-strict one.

• SCC:

∀s ∈ {1, . . . , M }, v i ∈ V vj ∈V :j =i A ij Y js ≥ Y is ⌊ k i 2 ⌋ + 1 . (24)
Indeed, it follows from the definition of SCC that:

∀s ∈ {1, . . . , M }, v i ∈ V vj ∈V :j =i A ij Y js ≥ k i - vj ∈V :j =i A ij Y js + 1, (25)
which expresses the fact that the in-degree of vertex v i is strictly greater than the out-degree, i.e., than the degree minus the in-degree. It is valid for unweighted graphs where all coefficients A ij are integer.

A few algebraic manipulations lead to:

∀s ∈ {1, . . . , M }, v i ∈ V vj ∈V :j =i A ij Y js ≥ ⌊ k i 2 ⌋ -(1 -Y is ) ⌊ k i 2 ⌋ + Y is , (26)
as it is easily checked for both cases Y is = 1 and Y is = 0. Formula (24) follows.

• ASCC:

The almost-strong cohesion condition can be expressed in terms of the variables Y as follows:

∀s ∈ {1, . . . , M }, v i ∈ V | k i = 2 vj ∈V :j =i A ij Y js ≥ Y is ⌊ k i 2 ⌋ + 1 (27) ∀s ∈ {1, . . . , M }, v i ∈ V | k i = 2 vj ∈V :j =i A ij Y js ≥ Y is . (28)
• SSCC:

The semi-strong cohesion condition can be expressed in terms of the variables Y as follows:

∀s, t ∈ {1, . . . , M } | s = t, v i ∈ V j∈V :j =i A ij Y js ≥ vj ∈V :j =i A ij Y jt + 1 -(1 -Y is )(k i + 1). ( 29
)
Indeed, consider the following two cases: (i) Y is = 1, that is vertex v i belongs to community s. Then the left-hand side term of ( 29) is equal to the in-degree of v i , the first term of the right-hand side represents the part of the out-degree of v i corresponding to edges with extremities in community s and t = s. As the last term disappears, the condition expresses that this partial out-degree must be strictly smaller than the in-degree of v i . As similar conditions hold for all other communities, it is clear that such a relation holds for the community for which the partial out-degree of v i is largest. (ii) Y is = 0. Then the right-hand side of ( 29) is non-positive and the condition is verified.

• WCC:

The weak cohesion condition can be written as follows in terms of variables X and Y :

∀s ∈ {1, . . . , M } 4 r∈E X rs ≥ vi∈V k i Y is + 1. ( 30
)
Indeed, the sum of in-degrees for community s may be written as 2 r∈E X rs (where the factor 2 is due to edges having both vertices in the community) and must be greater than the sum of out-degrees of community s, that is the sum of all the degrees minus the sum of in-degrees for vertices of that community:

vi∈V k i Y is -2 r∈E X rs . • EWCC:
The extra-weak cohesion condition was proposed in [22] without a mathematical expression. It can be written as follows:

∀s, t ∈ {1, . . . , M } | s = t, 2 r∈E X rs ≥ r={vi,vj }∈E (Y is Y jt + Y js Y it ) + 1. ( 31
)
The left-hand side of ( 31) is equal to twice the number of edges in community s. The right-hand side is equal to the number of edges with an end vertex in s and the other in community t. This expression can be linearized introducing ∀r

= {v i , v j } ∈ E non-negative variables Z r = Y is Y jt and Z ′ r = Y js Y it : ∀s, t ∈ {1, . . . , M } | s = t, 2 r∈E X rs ≥ r∈E (Z r + Z ′ r ) + 1, ( 32
)
and adding linearization constraints ∀s, t ∈ {1, . . . , M } | s = t:

Z r ≤ Y is (33) Z r ≤ Y jt (34) Z r ≥ Y is + Y jt -1 (35) Z ′ r ≤ Y js (36)
Z ′ r ≤ Y it (37) Z ′ r ≥ Y js + Y it -1. ( 38
)

Results

We first compare the results of the simulated annealing heuristic of Wang et al. [START_REF] Wang | Lecture Notes in Operation Research, The Third International Symposium on Optimization and Systems Biology (OSB09)[END_REF] with those of our exact algorithms [16]. Three problems are common to both lists of problems solved in [16] and [START_REF] Wang | Lecture Notes in Operation Research, The Third International Symposium on Optimization and Systems Biology (OSB09)[END_REF]: karate, dolphins, political books. In [START_REF] Wang | Lecture Notes in Operation Research, The Third International Symposium on Optimization and Systems Biology (OSB09)[END_REF], results are presented for direct optimization and for considering the weak cohesion condition through a constrained model. In both cases, the value of modularity Q is given together with the number of communities of the optimal partition found and the number of such communities which satisfy the weak cohesion condition. As simulated annealing is a probabilistic optimization heuristic, results depend on the seed used in its random number generator. The Q values are given with limited precision, i.e., only 2 significant digits. Comparing first the results of a direct optimization with the optimal partitions found in [16], we observe that:

• the number of communities and the proposed solutions differ from the optimal ones in several cases: dolphins (4 instead of 5), political books (4 instead of 5); • the optimal solution is found for the smallest problem, that is karate.

G-2014-11

Les Cahiers du GERAD Turning now to optimization subject to weak cohesion constraints, again the same three problems are common to the list of those solved by the algorithms of the present paper (see Table 2) and the constrained optimization one of [START_REF] Wang | Lecture Notes in Operation Research, The Third International Symposium on Optimization and Systems Biology (OSB09)[END_REF]. Results are as follows. For karate, the number of communities found in [START_REF] Wang | Lecture Notes in Operation Research, The Third International Symposium on Optimization and Systems Biology (OSB09)[END_REF] is optimal but the solution is not as its Q value is 0.40 and the optimal solution is 0.41979 with the weak condition satisfied by all communities. For dolphins and political books, the number of communities found is optimal, while the modularity value is 0.52 in [START_REF] Wang | Lecture Notes in Operation Research, The Third International Symposium on Optimization and Systems Biology (OSB09)[END_REF] versus 0.526799 with the algorithm of the present paper for dolphins, and 0.53 versus 0.526938 for political books (so, up to the precision of two digits, the solution found in [START_REF] Wang | Lecture Notes in Operation Research, The Third International Symposium on Optimization and Systems Biology (OSB09)[END_REF] is optimal). In Table 2 we report the results obtained by unconstrained modularity maximization and modularity maximization subject to the weak constraint and the extra-weak constraint respectively. It appears that the only three cases in which results obtained with the weak constraints differ from the unconstrained case are dolphins, p53 protein and political books (represented in bold in Table 2).

We now discuss these three cases in more detail.

• dolphins: there are 5 communities with a modularity of 0.528519 in the unconstrained solution while there are 4 communities with a modularity of 0.526799 in the constrained solution. Detailed memberships of both solutions are given in Tables 3 and4. It appears that the smallest community, that is C 3 , is split among three of the four other communities: dolphin 40 goes to community C w 2 , dolphins 4, 9, and 60 go to C w 3 , and dolphin 37 to C w 4 . However, these three assignments entail further changes: dolphins 21 and 45 move from C 1 to C w 4 while dolphins 54 and 62 move from C 5 to C w 1 . Thus, nine dolphins change community when the weak condition constraint is added. • p53 protein: there are 7 communities with a modularity of 0.535134 in the unconstrained solution while there are 6 communities with a modularity of 0.534488 in the constrained solution. Both partitions are depicted in Fig. 3. It appears that a small cluster with five entities in blue in the center of the figure does not satisfy the weak cohesion condition. This community is split, four entities going to the cluster below it and one entity to the cluster above. • political books: there are 5 communities with a modularity of 0.527237 in the unconstrained solution while there are 4 communities with a modularity of 0.526938 in the constrained solution. Both partitions are depicted in Fig. 4. It appears that the smallest community of the unconstrained solution, that is the community containing books 49, 50, and 58 depicted in orange on the right side of the figure, does not satisfy the weak cohesion condition. All of these three books move to the community above it.

It appears from Table 2 that the only case in which the modularity maximization constrained with the extra-weak constraint differs from the unconstrained solution is political books, for which the constrained optimal partition contains 4 communities.

The average reduction in modularity when adding extra-weak constraints is of 0.000027, i.e., it is very moderate. Reduction in modularity when imposing the weak condition is slightly larger, i.e., 0.000242.

The effect of imposing the other cohesion conditions is much larger, as shown in Table 5. For the strong and semi-strong conditions, average modularity is reduced from 0.521334 to 0.354608; for the almoststrong conditions from 0.521334 to 0.486862. We also observe drastic reductions in the average number of communities, that is reduced from 5.090909 to 2.727273 for the strong and semi-strong conditions, from 5.090909 to 3.90909 for the almost-strong condition. The effect of imposing the strong and the semi-strong conditions are the same for the considered datasets, for which the same optimal partitions are obtained. The effect of introducing almost-strong constraints instead of strong ones is less drastic, as already observed in [23].

Conclusions

The effect of adding five kinds of cohesion conditions to a modularity maximization problem has been studied. After reviewing the strong, almost-strong, semi-strong, weak, and extra-weak cohesion conditions, we check whether they were satisfied by all communities in maximum modularity partitions of some known real world problems. It appears that the strong and semi-strong cohesion conditions are very strict. Indeed, all of the communities of an optimal partition satisfy these conditions for one dataset only. The almost-strong cohesion condition is slightly less strict, but still quite stringent. The weak cohesion condition is more intuitive and satisfied by all communities of 8 out of the 11 optimal partitions. The extra-weak cohesion condition appears to be very lax as there is only a single community of only one of the optimal unconstrained partitions which does not satisfy it.

Figure 1 :

 1 Figure 1: (Color online) (a) Partition satisfying the conditions WCC and EWCC, but not the conditions SCC and SSCC. (b) Partition satisfying the condition EWCC but not the condition WCC.

Figure 2 :

 2 Figure 2: (Color online) Partitions obtained on the dolphins network using the original modularity maximization model (left) and the modularity maximization model with the weak cohesion constraint (right).

Figure 3 :

 3 Figure 3: (Color online) Partitions obtained on the p53 network using the original modularity maximization model (left) and the modularity maximization model with the weak cohesion constraint (right).

Figure 4 :

 4 Figure 4: (Color online) Partitions obtained on the political books network using the original modularity maximization model (left) and the modularity maximization model with the weak cohesion constraint (right).

Table 2 :

 2 Modularity maximization with weak and extra-weak cohesion conditions. Network dimension (n = number of vertices, m = number of edges), number of clusters found (M , M w , M ew ) and corresponding modularity value (Q, Q w , Q ew ) for the standard modularity maximization problem and the modularity maximization problem with weak and extra-weak constraints.

	network			modularity maximization	weak	extra-weak
	dataset	n	m	M	Q	Mw	Qw	Mew	Qew
	strike	24	38	4	0.561981	4	0.561981	4	0.561981
	karate	34	78	4	0.41979	4	0.41979	4	0.41979
	Korea1	35	69	5	0.477736	5	0.477736	5	0.477736
	Korea2	35	84	5	0.450822	5	0.450822	5	0.450822
	sawmill	36	62	4	0.550078	4	0.550078	4	0.550078
	dolphins small	40	70	4	0.620714	4	0.620714	4	0.620714
	graph	60	114	7	0.502655	7	0.502655	7	0.502655
	dolphins	62	159	5	0.528519	4	0.526799	5	0.528519
	Les Misérables	77	254	6	0.560008	6	0.560008	6	0.560008
	p53 protein	104	226	7	0.535134	6	0.534488	7	0.535134
	political books	105	441	5	0.527237	4	0.526938	4	0.526938
	average			5.090909	0.521334	4.818182	0.521092	5	0.521307

Table 5 :

 5 Modularity maximization with strong, almost-strong, and semi-strong condition. Network dimension (n = number of vertices, m = number of edges), number of clusters found (M s , M as , M ss ) and corresponding modularity value (Q s , Q as , Q ss ) for the standard modularity maximization problem and the the modularity maximization problem with strong, almost-strong and semi-strong constraints.

	network			modularity maximization		strong	almost-strong	semi-strong
	dataset	n	m	M	Q	Ms	Qs	Mas	Qas	Mss	Qss
	strike	24	38	4	0.561981	2	0.257271	3	0.54813	2	0.257271
	karate	34	78	4	0.41979	2	0.132807	4	0.402038	2	0.132807
	Korea1	35	69	5	0.477736	4	0.383638	4	0.383638	4	0.383638
	Korea2	35	84	5	0.450822	3	0.424036	4	0.432469	3	0.424036
	sawmill	36	62	4	0.550078	4	0.550078	4	0.550078	4	0.550078
	dolphins small	40	70	4	0.620714	3	0.573571	4	0.620714	3	0.573571
	graph	60 114	7	0.502655	1	0	4	0.438135	1	0
	dolphins	62 159	5	0.528519	2	0.359242	3	0.480598	2	0.359242
	Les Misérables	77 254	6	0.560008	4	0.437868	6	0.52921	4	0.437868
	p53 protein	104 226	7	0.535134	2	0.284204	4	0.472502	2	0.284204
	political books 105 441	5	0.527237	3	0.497969	3	0.497969	3	0.497969
	average			5.090909	0.521334	2.727273 0.354608 3.909091 0.486862 2.727273 0.354608

We then show how these five cohesion constraints can be added, one at a time, to a quadratic convex mixedinteger optimization formulation proposed for modularity maximization. For three real world datasets, we discuss in detail the modifications in the optimal solutions due to the imposition of weak cohesion constraints.

Adding cohesion constraints to a modularity maximization model appears to be feasible and in particular for the weak condition yields intuitive and appealing results. While the present work is based on the model of Xu et al. [14], we will try to explore the addition of cohesion constraints to other models for modularity maximization.

S. Cafieri has been supported by the French National Research Agency (ANR) through grant ANR 12-JS02-009-01 ATOMIC. A. Costa has been supported by IDC Grant IDG21300102. G-2014-11 9 the standard modularity maximization model on dolphins dataset.