S Di Tommaso

P Rotureau

B Sirjean

R Fournet

W Benaissa

P Gruez

C Adamo

A mechanistic and experimental study on the diethyl ether oxidation

Keywords:

This work presents the results of the theoretical investigations on autoxidation process of diethyl ether (DEE), a chemical largely used as solvent in laboratories and considered to be responsible for various accidents. Based on Density Functional Theory (DFT) calculations, the aims of this study were the identification of all the most probable reaction paths involved in DEE oxidation (at ambient temperature and under conditions that reflect normal storage conditions) and the characterization of products and all potential hazardous intermediates, such as peroxides. Results indicate that industrial hazards could be related to hydroperoxide formation and accumulation during the chain propagation step. A detailed kinetics model of DEE oxidation in the gas phase was then developed from all energetic and kinetics parameters collected during the mechanistic study. Outputs of the kinetics model, in terms of time of evolution of product concentrations, have been then compared with the experimentally measured concentration of products (notably hydroperoxides) issued from tests on DEE oxidation conducted under accelerated conditions with autoclaves.

Introduction

A great number of organic compounds spontaneously decompose by a free-radical reaction of the carbon chain with molecular oxygen, in a self-propagating process of auto-oxidation that may generate a large variety of peroxide molecules. Many laboratory accidents can be ascribed to the presence of such peroxides in stored chemicals [1][2][3][4], often under poor conditions, for over long periods. Ethers are typical solvents showing such a behavior. Dimethyl ether (DME) is the most investigated ether of this class both experimentally [5,6] and theoretically [7], since it has a relevant role, due to its good ignition properties, in the domain of biofuels. In particular, Andersen and coworkers [7] proposed a DFT study on the process of low temperature oxidation of this ether. Experimental studies also exist on the oxidation of the other short-chains ethers suggesting oxidation products or intermediates or giving some hints on the kinetics of the process. It is the case of diethyl ether (DEE) [8][9][10][11][12], diisopropyl ether (DIPE) [11,13] or methyl tert-butyl ether (MTBE) [14]. However few investigations have been devoted to the mechanistic elucidation of their oxidation under usual storage conditions [6,[START_REF] Di Tommaso | [END_REF]. Moreover, the chemistry of the oxidation of saturated branched ethers has some similarities with the well known of alkanes [START_REF] Semenov | Chemical Kinetics and Chain Reactions[END_REF], as suggested for DME by Andersen and Carter [6]. This has also been demonstrated for diethyl ether in our latest extensive theoretical study on its low temperature and low pressure oxidation [START_REF] Di Tommaso | [END_REF]. In this scientific context, where incomplete experimental data are available and only a sketchy picture of the process emerges, the present work proposes an in-depth study of DEE oxidation carried out at both theoretical and experimental levels to better understand the process and manage the accidental risks linked to diethyl ether use and storage. In the following, starting from a brief presentation of the results of the mechanistic study (based on quantum chemistry calculations) of all the possible reactions involved in the process, the development and the validation of the kinetics model obtained for the gas phase oxidation of DEE will be discussed. At the end, the issues of kinetics model, in terms of concentrations of products and time of reactions, will be compared with experimental data obtained on DEE oxidation tests carried out in accelerated conditions.

Materials and Methods

Molecular modeling

All the Density Functional Theory (DFT) calculations have been performed using the B3LYP functional [17] and Gaussian 03 program [18]. The 6-31+G(d,p) basis set was used to optimize structures and for subsequent frequency calculations to characterize stationary points as minima or transition states. Intrinsic reaction coordinate (IRC) calculations were also performed to verify that identified products and reactants were correctly connected [19]. For radicals (open shell systems), unrestricted DFT calculations were carried out and the obtained results were checked for wave-function instability. Enthalpy values relative to reactions discussed in the following are collected as activation barriers (ΔH ≠ ) and product stabilization (ΔH) in table I. The relative enthalpies are referred to a single reaction step, that is reactants-products or reactants-transition state.

Rate constants calculations

Rate constants for each elementary reaction were calculated using classical transition state theory (TST) [20]. The tunneling effect was taken into account for reactions involving inter or intramolecular hydrogen transfer using the transmission coefficient in the form proposed by Wigner [21] or Skodje and Truhlar [22] or Eckart [23]. Activation enthalpies have been computed as the difference between enthalpies of formation (with respect to isolated atoms) of all reactants and related transition states. The kinetic parameters were obtained by fitting the rate constant values obtained from TST at several temperatures between 500 and 2000 K with the equation Where A, n, and Ea are the parameters of the modified Arrhenius equation and k ∞ is the highpressure-limit rate constant.

Experimental tools

Storage of chemicals is generally performed at room temperature. Under these conditions, the kinetics are slow and real time tests can take several weeks. As the kinetics are accelerated by temperature, an experimental setting has been designed in order to store DEE at high temperature. Contrary to Naito et al. [11] who used Accelerated Rate Calorimeter (ARC) bombs of 10 mL, experiments were carried out using six 50 mL autoclaves: each stainless steel vessel had previously undergone a passivating step and equipped with a temperature sensor and a pressure transducer. The six autoclaves are inserted in a heating block controlling the temperature (cf. figure 1).

Diethyl ether (99.7%) was obtained from VWR and was used without further addition of a stabilizer. Fourteen milliliters of DEE were introduced into each autoclave and pressurized by synthesized air from a gas cylinder at an approximately pressure of 10 bar absolute. The samples were then heated to two temperatures: 80°C and 100°C. Temperatures and pressures inside the autoclaves were recorded throughout the experiment. At determined storage times, the samples were recovered in order to determine the amount of peroxide formed. Concentration of peroxide was measured using a potentiometric titrator and a solution of potassium permanganate. This method gives the global amount of peroxide function (equivalent H 2 O 2 ).

Results and discussions

Molecular modeling

In figure 2 is depicted a scheme summarizing the main reaction pathways of the diethyl ether (R) oxidation issued of DFT mechanistic study of the process [START_REF] Di Tommaso | [END_REF]. Analogously to the accepted alkane autoignition mechanism [START_REF] Semenov | Chemical Kinetics and Chain Reactions[END_REF][START_REF] Denisov | Oxidation and Antioxidants in Organic Chemistry and Biology[END_REF] the oxidation process starts with an initiation step in which a radical R is produced. This species can simply undergo decomposition. Alternatively, a chain propagating cycle is established if the R radical reacts with molecular oxygen to give a ROO species. This peroxy radical can then react with a new DEE molecule to produce a hydroperoxide ROOH and regenerate a R radical. Chain termination ensues when two radical intermediates couple to each other (e.g. dimerization) with loss of the radical character. By contrast, radical isomerization may occur without interrupting the chain reaction. This schematic picture is complicated by the possible decomposition of the hydroperoxides that represent the stoichiometric products of the ideal chain-propagating cycle.

Figure 2.Global scheme of the proposed mechanism of DEE (R) oxidation.

The experimentally observed products are also reported Among the various chain initiators that have been proposed in literature, OH radical is, from our calculations, the most efficient one from an energy point of view [START_REF] Di Tommaso | [END_REF]. Competition among different pathways starts with the formation of the radical R . In conditions of low (null) oxygen concentration, its decomposition, characterized by a high rate constant (k=1.29s -1 at 298 K, as reported in figure 3), leads to the formation of acetaldehyde, formaldehyde and ethanol. In presence of oxygen, the production of peroxy radicals (ROO ) becomes important and the corresponding reaction is strongly exothermic (ΔH=-31.8 kcal/mol) and barrierless. The resulting ROO radicals can evolve along three competitive reaction pathways: hydroperoxide production, isomerization and dimerization (see figure 2). Hydroperoxide production is endothermic by 12.8 kcal/mol and has an activation barrier of approximately 13.3 kcal/mol. Nonetheless, it becomes probable in the usual conditions of autoxidation, where DEE concentration is high. Indeed, this reaction channel leads to some of the major products identified experimentally (acetic acid, acetaldehyde, ethanol) [10]. On the other hand, the isomerization reaction features significant product stabilization and activation energies comparable to that of hydroperoxide production. Moreover, it could account for the formation of some of the observed products (acetaldehyde, acetic acid, peracetic acid) [10].

However the two pathways are characterized by different entropic effects, so that the final Gibbs free energies are 23.7 and 19.3 kcal/mol for the hydroperoxide production and isomerization, respectively. These different values of Gibbs free energy enthalpy lead to a very important difference in rate constant values, as reported in figure 3. Finally, dimerization is thermodynamically favored due to significant product stabilization. However, the activation energies along the following steps of these paths are higher than in other reactions, the lowest energy barrier being 32.3 kcal/mol. Furthermore, dimerizations are bimolecular reactions, and their kinetics are therefore highly sensitive to the concentration of peroxy radicals in the reaction medium. In fact the dimerization pathway has been suggested in literature as the major reaction channel for DEE oxidation in atmospheric conditions [12]. In solution, dimerization happens if the two peroxy radicals are close, so that the reaction is limited by the diffusion of two highly reactive species: this channel could be considered as improbable.

Therefore, among all the characterized pathways, decomposition of R and isomerization could be the most relevant in solution.

Moreover, from a hazard point of view, different hydroperoxides have been identified as products or sub products of principal reaction pathways. Considering the high activation barriers for decomposition of these compounds, their accumulation during DEE oxidation could be a (major) risk source during the storage of the chemical.

Figure 3.Sketch of the principal reaction pathways identified. For each reaction the activation free Gibbs energy and the rate constant are specified

Detailed chemical kinetics model

The approach used for the development of the chemical kinetics model of the process is summarized in figure 4.

Figure 4. Scheme of the approach used for the development of the kinetics model

The first step was the automatic generation of a detailed mechanism of the oxidation process of DEE with the software EXGAS [START_REF] Warth | [END_REF]. The mechanism produced contains thermochemistry of each species involved in the process in the form of coefficients of NASA polynomial (calculated by the software THERGAS [26] based on Benson group additivity methods) and the kinetics of each reaction considered (calculated by the software KINGAS). All thermodynamic and kinetics data issued from the DFT mechanistic study [START_REF] Di Tommaso | [END_REF] were introduced in this mechanism. The final DEE oxidation model contains 158 species and 355 reactions. Simulations were performed using the SENKIN code of the CHEMKIN II [START_REF] Kee | CHEMKIN II. A Fortran Chemical Kinetics Package for the Analysis of a Gas-phase Chemical Kinetics[END_REF] program suite assuming homogeneous gas-phase reactions in an isochoric static reactor. Temperature and pressure conditions can be selected as inputs to match target experimental data. The mechanism produced was developed and validated by the comparison with experimental literature data in terms of evolution of pressure and product concentration. In particular, because of the lack of experimental data on the DEE oxidation in the liquid phase at room temperature, the model was developed for the process in the gas phase using experimental data (concentrations of products in function of time at 153°C) from experimental work of Waddington [10]. It was then validated from comparison with the Lemay and Ouellet study [9] performed at a different temperature than Waddington. Figure 5 reports the comparison between experimental data obtained by Lemay and Ouellet [9] in their study on DEE gas phase oxidation at 168°C and the curves generated from the CHEMKIN simulation (at T=172°C) in the developed kinetics model.

Figure 5. Evolution of partial pressures of oxygen, acids and peroxides at T=172°C by CHEMKIN simulation. Experimental curves [9] have been displaced by about 10 minutes to facilitate the comparison

Some considerations can be made on the basis of this graph:

A displacement in the observation of reactivity of the mixture between experiment and simulations is observed. In particular DEE is consumed earlier (approximately 10 minutes) in the experiment at a lower temperature (168°C instead of 172°C in simulation). We have to note that experimental parameters are not accurately detailed in the studies [9,10] used as references for the development of the kinetics model. If we consider temperature, for instance, in both cases given values are referred to the temperature of the oven in which the reactor is suspended and not to the one in the reactor, and there is no reference to the adiabaticity of the reaction vessel.

After the beginning of reactivity, the trends are well reproduced by the kinetics model as well as the selectivity of the mixture of products.

In the first period of the process, the induction time, is probably characterized by heterogeneous reactions with the reactor surface, which are not taken into account in simulations. Therefore we can conclude that the model is validated for DEE oxidation in the gas phase in terms of reproduction of product concentrations. Because of ill-defined experimental conditions in the reactor and interferences of heterogeneous wall reactions, it is difficult to make firm conclusions about the difference of reactivity observed between experiments and simulations. Consequently, temperature and residence time adjustments remain acceptable.

Experimental results

The peroxide concentration vs. time profiles are given in figure 6. These results show two tendencies: initially the peroxide concentration increases to a maximum and then decreases to reach very low values. At 80°C, the maximum is around 6000 ppm after 45 hours. By increasing the temperature from 80°C to 100°C, the maximum is reached after only 5 hours and the induction time (time before peroxides are detected) is much shorter.

Figure 6.Peroxides concentration evolution at two operating temperatures

Figure 7 gives the pressure evolution for the experiment carried out at 80°C: pressure decreases in the first period, when peroxides are formed, and increases when the peroxide concentration decreases. According to Naito et al. [11], it is assumed that the pressure decrease corresponds to the oxygen consumption by the ether, the beginning of the autoxidation chain reaction, and the pressure increase corresponds to the formation of volatile products from peroxide decomposition. It has been indeed verified that the pressure drop is approximately equal to the initial concentration of oxygen in the autoclave. In this configuration and at this temperature, less than 1 mL of DEE should be vaporized in the autoclave. Nevertheless, the DEE is in excess relative to oxygen in gas phase. Thus the reaction could then occur either in the gas phase or liquid phase.

Further experiments need to be performed to determine more precisely the kinetics and the mechanism:

Other temperatures should be tested (for instance 60°C and 110°C) in order to extrapolate the kinetics to lower temperature. Other products identified in the molecular modeling study should be quantified (ethyl acetate, ethanol, ethyl formate, methanol, acetaldehyde, methanol…). The influence of the initial volume of DEE could give some indications about the phase where the reaction occurs. Various initial oxygen concentrations should be tested to confirm the correlation between the oxygen consumption and the peroxide production.

Nonetheless, the overall objective of this work is to couple the experimental results with the results produced by molecular and kinetic modeling to obtain a kinetics model of the oxidation of DEE in storage conditions, allowing the detection of its instability at ambient temperature.

Kinetics model vs. experiments

With this aim, data generated by CHEMKIN simulations in the kinetics model developed and validated for the gas phase oxidation of diethyl ether were compared with the experimental peroxide concentration curves described in the previous paragraph. In particular, all simulation parameters (initial reactor pressure, reactor volume, reactant concentrations) were chosen to correspond to actual experimental conditions, and the concentration of hydroperoxides was recorded. The simulation was made using a temperature of 107°C and compared with the peroxide concentrations measured during the experiment at 100°C. This comparison is reported in figure 8.

Figure 8. Comparison between peroxides concentration obtained from kinetics simulation at 107°C and experiments (at 100°C). Experimental curve has been displaced by about 12 hours to facilitate the comparison

The results from the simulation (at 107°C) predict a peak ROOH concentration coming 12 hours after the actual time of the peak observed in the experiment (at 100°C). In figure 7, the experimental results have been shifted in time to facilitate a comparison of the ROOH concentration. However, as seen before for the validation of the gas phase model, the ROOH concentrations predicted by the simulation and observed experimentally are quite similar. These preliminary results encourage further improvement of the kinetics model using new experimental data obtained either under accelerated conditions or at ambient temperature with the goal to better reproduce the true behavior of DEE at storage conditions.

Conclusions

Two different but complementary theoretical approaches have been used to study the oxidation process of diethyl ether: an in-depth theoretical study based on quantum chemistry calculations and a kinetics study. Concerning the reaction mechanism, both approaches show that in presence of molecular oxygen, reaction pathways involving O 2 are favored with respect to unimolecular decomposition. The kinetics model has been obtained and validated for the oxidation of DEE in the gas phase from experimental data obtained from the literature: the model reproduces well the selectivity of product blend but is less efficient for reaction time. This problem in predicting the time of reaction was also observed when comparing peroxide concentrations obtained at approximately 100°C from simulation and from preliminary oxidation tests performed with autoclaves at INERIS. However, the kinetics model is able to predict suitably the concentration of peroxides produced during the oxidation process of DEE. We can expect significant improvement in the kinetics model by taking into account the results of ongoing tests to quantify other molecular species. The objective of experimental study is also to extrapolate measurements of DEE oxidation products towards conditions of normal storage. From a risk point of view, this model is expected to determine accurately both peroxide concentrations and time of maximal concentration in order to predict potential instability and better manage industrial and laboratory risks of this chemical.

Figure 1 .

 1 Figure 1.Experimental setting for experiments at accelerated conditions

Figure 7 .

 7 Figure 7.Peroxides concentration and pressure evolution at 80°C

Table I .

 I Activation barriers in terms of enthalpies and Gibbs free energies (respectively ΔH ≠ and ΔG ≠ in kcal/mol), products stabilization (ΔH, kcal/mol) and rate constants (k, s -1 ) obtained for the three main reactions of the DEE oxidation mechanism.

	REACTION	ΔH ≠	ΔH	ΔG	k
	R decomposition	19.2	4.0	17.7	1.29
	ROO isomerization	17.2	14.8	19.3	0.88
	ROOH production	13.3	12.8	23.7	2.3 10


Guillaume Fayet
email: guillaume.fayet@ineris.fr

Patricia Rotureau

Benoît Minisini

Benoit Minisini

Decomposition mechanisms of trinitroalkyl compounds: a theoretical study from aliphatic to aromatic nitro compounds

Keywords: Nitro compounds, Decomposition mechanisms, Density Functional Theory

teaching and research institutions in France or abroad, or from public or private research centers.

Decomposition mechanisms of trinitroalkyl compounds :

a theoretical study from aliphatic to aromatic nitro compounds

INTRODUCTION

Nitro compounds are substances recognized as potentially hazardous due to the presence of nitro groups. Besides, some of them are commonly used in high explosive compositions 1 . So, particular attention is needed for the storage, transport and handling of such compounds.

For this reason, knowledge on the hazardous properties of nitro compounds is needed. In recent years, works have also been conducted for the understanding and prediction of physico-chemical hazards based on theoretical approaches [2][START_REF] Kamlet | Sixth Symposium (International) on Detonation[END_REF][START_REF] Morrill | [END_REF][5][6][7][8][9] . Indeed, these computational methods became powerful with the improvement of computer capabilities, allowing access to new information at molecular level.

For instance, quantitative structure property relationship models (QSPR) were developed to predict the hazardous properties of energetic compounds from descriptors of their molecular structures 2,[START_REF] Kamlet | Sixth Symposium (International) on Detonation[END_REF]5,[10][11][12][13][14] . In particular, models were developed in a previous work for the prediction of the impact sensitivities of nitro compounds 11 . High predictive powers were obtained for nitramine and nitroaliphatic compounds with R²=0.85 and 0.88, respectively, as evidenced on external validation sets. However, early efforts to predict impact sensitivities of aromatic compounds reveal fairly poor, achieving R² values of 0.54 only. This failure was proposed to be related to a too large diversity of compounds in the training data set regarding its small size. Indeed, some subgroups of molecules among the aromatic nitro compounds of the data set were identified to present special decomposition features. This was in particular the case for compounds presenting -(CH 2 ) n -C(NO 2 ) 3 aliphatic fragments. Similar analysis had been already proposed to distinguish between ortho and non ortho substituted nitrobenzene derivative for the prediction of their heats of decomposition 15 . Indeed, as in the case of o-nitrotoluene derivatives 16 , the interaction of nitro groups with the ortho groups can induce a specific decomposition pathway. Besides, a model with high predictive capabilities was achieved by considering only the compounds presenting no substituent in the ortho position to the nitro group that followed the same decomposition process which is the dissociation of the carbon-nitro bond 17 .

So, one purpose of the present study was to identify, in the same way, key features for compounds presenting -(CH 2 ) n -C(NO 2 ) 3 fragments that would have caused troubles in the development of an accurate QSPR model for the impact sensitivity of aromatic nitro compounds 11 . For this reason, the mechanisms involved in the decomposition of compounds presenting C(NO 2 ) 3 fragments were studied, starting with the aliphatic 1,1,1-trinitrobutane to aromatic compounds presenting -(CH 2 ) n -C(NO 2 ) 3 fragments.

The thermal decomposition of nitro compounds has been studied in numerous experimental [18][19][20][21] and computational studies 16,17,[22][23][24][25][26] . Different decomposition paths upon the type of nitro compounds under investigation 27 were thus suggested, and in particular for nitroalkanes 28 . The main reactions extracted from literature focused on the decomposition of nitroalkanes are presented in Figure 1.

Reaction 1 consists in the homolytic dissociation of the C-NO 2 bond. This reaction is commonly considered as the main reaction path for most of the nitro compounds 27 . It was in particular considered experimentally by Cottrell et al. as the primary mechanism of the thermal decomposition of trinitromethane in the range 380-430°C 29 . Using theoretical analysis 30 , Kiselev came to the same conclusion. In various works 25,31 , this reaction was even only used to characterize the decomposition of nitroalkanes.

Reaction 2 issues from a nitro-nitrite isomerization that was considered in various experimental or theoretical studies 18,24,26,30,32 . It was notably characterized theoretically by Khrapkovskii et al. for a series of nitroaliphatic compounds 32 . Nevertheless, most of these studies considered it as non competitive with reactions 1 and 3. Reaction 3 is the HONO elimination. This channel has been considered in various studies of nitroaliphatic compounds [33][34][35] notably suggested by pyrolysis experiments on nitropropanes 36 . More recently, Pauler et al. 37 proved theoretically that this mechanism was favored versus reaction 1 at low temperatures (reaction 1 being favored at higher temperatures (> 150°C)) for dinitropropyl compounds, with computational results in accordance with experiments.

Reaction 4 is the mechanism by elimination of O that was early proposed by Taylor for the decomposition of nitromethane 38 without being noticed in later studies. Reaction 5 is the breaking of the C-C bond. Even if no recent studies highlighted this channel among the main mechanisms involved in the decomposition of nitroalkanes, it represents the typical mechanism involved on the decomposition of alkanes [39][40][41] .

Focusing on trinitroalkanes, Kisilev 30 studied the decomposition of trinitomethane using quantum chemical calculations. The C-N bond breaking was proposed to be the primary mechanism. However, the reactions 1 and 2 were only considered whereas the HONO elimination (reaction 3) was proposed to be critical in the decomposition of trinitroalkanes in the experimental works of Nazin et al. 42,43 . Moreover, no study considered the five reactions 1 to 5 in a same study.

In the present study, an extended theoretical study was performed, using the density functional theory, to clarify the decomposition mechanisms of trinitroalkyl compounds. Reactions 1-5 were computed for the 1,1,1-trinitobutane, chosen as a reference molecule, and a series of derivatives to evidence the influence of structural features, e.g. presence of substituents in various positions, or the carbon chain length. These reactions were also compared to the mechanisms involved on the nitro groups on the aromatic cycle for nitroaromatic compounds. It allowed to evidence the specific mechanism of decomposition of compounds presenting trinitroalkyl fragments (-R-C(NO 2 ) 3 ) among aromatic nitro compounds since they would be more likely considered as nitroaliphatics regarding the initiation of their decomposition process.

COMPUTATIONAL DETAILS

All calculations were carried out using the density functional theory with the PBE0 hybrid functional 44 as implemented in Gaussian09 software [START_REF] Frisch | Gaussian09[END_REF] . All structures (reactants, products, local minimum intermediates and transition states) were optimized using the 6-31+G(d,p) basis set. The nature of stationary points was checked on harmonic frequencies computed at same level of calculation. The spin-unrestricted formalism was employed for the calculations on open shell molecules and spin contamination appeared to be negligible for all computed radical species (with S² values close to 0.75).

The PBE0/6-31+G(d,p) level of theory was already successfully used for nitro systems 16,17 . Good agreement was highlighted in table S1 for the calculation of the bond dissociation energies of three C(NO 2 ) 3 compounds with experimental data. Moreover, table S1 (in Supporting Information) exhibited the accuracy of this method for the trinitromethane (see in Figure 2) on a geometrical point of view when compared to experimental data [START_REF] Sadova | [END_REF]47 and results obtained with larger basis set 48 .

In this study, all energetic values were considered in terms of Gibbs Free Energies (ΔG) since entropic contributions were of importance as already demonstrated for nitroaromatic compounds on similar reactions in previous works 16,17 . This ΔG values were issued from harmonic frequencies using the standard statistical approaches, as implemented in Gaussian09 software.

To quantify the effect of the substituent groups on the different energetic values, different global reactivity descriptors were evaluated, e.g. the chemical hardness () 49 , the electronegativity () 50 and the electrophilicity index (ω) 51 .

χ = (I+A)/2 (1) η = (I-A)/2 (2) ω = χ² / 2η (3)
where I and A are the potential ionization and the electron affinity of the molecules. A second approximation was finally used since the opposite of HOMO ( HOMO ) and LUMO ( LUMO ) energies were used to evaluate the vertical potential ionization, I  - HOMO , and vertical electron affinity, A  - LUMO;

respectively.

RESULTS AND DISCUSSION

Decomposition of the 1,1,1-trinitrobutane

The five decomposition mechanisms identified in Figure 1 were computed for the 1,1,1trinitrobutane (TNBt), as a reference molecule. Table 2 and Figure 3 summarize the energetic profiles obtained for these mechanisms. All calculated molecular structures (reactive, products, transition states and local minima) are represented in Figure 4.

The first mechanism consisted in the direct breaking of the carbon-nitrogen bond related to nitro groups. For the 1,1,1-trinitrobutane (S0), such dissociation was involved for two positions, a and b, as represented in Figure 1. However, the loss of both nitro groups led to the same radical (denoted Rad1) with a computed dissociation energy ΔG diss =21.3 kcal/mol. So, the breaking of the C-NO 2 bond appeared to be strongly energetically lower than for nitro groups linked to benzene rings (e.g. ΔG diss =58.2 kcal/mol for nitrobenzene 17 and 55.1 kcal/mol for o-nitrotoluene 16 ).

The second mechanism was initiated by isomerisation of the 1,1,1-trinitrobutane (S0) to a local minimum (LM1), corresponding to a nitrite form, via the transition state TS1. Then, the breaking of the O-NO bond in the nitrite fragment led to the formation of two radicals:

CH 3 -(CH 2 ) 2 -C(NO 2 ) 2 -O  (Rad2
) and NO  . The rate-limiting step of the mechanism was the NO 2 -ONO isomerisation through TS1. The Gibbs free energies ΔG TS1 are 60.0 kcal/mol and 58.1 kcal/mol for nitro groups a and b, respectively, i.e. higher than for the direct C-NO 2 bond dissociation and so non competitive with reaction 1 for the thermal decomposition (with a difference of 36.8 kcal/mol).

Elimination of HONO (reaction 3) was computed according to three molecular pathways. First ones issued from interaction of nitro groups a and b with hydrogen on the closest carbon (carbon 2) through transition states TS2a and TS2b. Then, the concerted breaking of the C-H and C-N bonds formed, in both cases, nitrous acid (HONO) and 1,1-dinitrobut-1-ene (CH 3 -CH 2 -CH=C(NO 2 ) 2 ), denoted Prod1. Activation Gibbs free energies (ΔG TS2a and ΔG TS2b ) of 62.3 kcal/mol and 37.7 kcal/mol were computed for nitro groups a and b, respectively. This implies that the reaction on the lateral nitro group is favoured for such mechanism. The third HONO elimination path was obtained from the 1,1,1-trinitrobutane in a cis conformation (S0') which is only 2.5 kcal/mol higher in ΔG than S0. Then, the nitro group b interacts with hydrogen on the terminal CH 3 (i.e. carbon 4) through the transition state TS3, leading to the formation of HONO and of a bi-radical product

• CH 2 -(CH 2 ) 2 -C • (NO 2 ) 2 ,
denoted Prod2. This reaction path was classified as non competitive with the pathway through TS2b with an activation energy of 86.0 kcal/mol.

Concerning the elimination of an oxygen from 1,1,1-trinitrobutane (reaction 4), the dissociation of nitrogen-oxygen bonds in nitro groups led to the formation of CH 3 -(CH 2 ) 2 -C(NO 2 ) 2 -NO  (Rad3). This mechanism was computed for both nitro groups a and b with respective activation energies ΔG Rad3 of 147.1 kcal/mol and 147.3 kcal/mol. So, both of them were not competitive.

The last mechanism consisted in the breaking of the carbon-carbon bond that formed trinitromethyl radical C  (NO 2 ) 3 and propyl radical (Rad4). The dissociation Gibbs free energy of this reaction was computed to be 74.7 kcal/mol. This mechanism was not competitive with the reactions involving nitro groups (i.e. reactions 1, 2 and 3) confirming the specificity of nitro compounds compared to simple alkanes.

Finally, Figure 3 demonstrated that the two most favored reaction paths for the decomposition of 1,1,1-trinitrobutane were the direct breaking of the carbon-nitrogen bond (reaction 1) and the HONO elimination (reaction 3) for nitro group b with activation Gibbs free energies of 21.3 kcal/mol and 37.7 kcal/mol, respectively. All other mechanisms exceeded an activation energy of 58 kcal/mol. This observation is in agreement with previous experimental statements that generally considered these two mechanisms as the most competitive reaction paths for the decomposition of nitroalkanes, as previously mentioned. For this reason, the next steps of this study focused on these two main mechanisms.

It has also to be noticed that in the present case, entropic contribution was of great importance. Indeed, considering Gibbs free energies (ΔG) allowed the identification of the breaking of the carbonnitrogen bond as the favored pathway whereas HONO elimination would have been considered as competitive based on ΔE or ΔH results. Indeed, a difference lower than 2 kcal/mol in the activation energies between the two reactions was calculated (considering ΔE or ΔH), as can be seen in table 3, according to the results issued in previous studies on nitroalkanes 22,24 .

The next step of the study concerned the influence of structural molecular features on the competition between the two main mechanisms highlighted for 1,1,1-trinitrobutane, i.e. the breaking of the C-NO 2 bond and the HONO elimination.

Influence of the length of the carbon chain

To investigate the influence of the length of the carbon chain, trinitroalkanes from trinitromethane (1 carbon atom) to trinitrooctane (8 carbon atoms) were considered. The main two decomposition mechanisms highlighted in the previous section (reactions 1 and 3) were computed for these 8 compounds and summarized in table 3.

Slight increases of activation Gibbs free energies were observed for both reactions for the two smallest molecules but no change was highlighted from trinitropropane (with 3 carbon atoms) for both mechanisms with ΔG diss =21.1-21.8 kcal/mol and ΔG TS2 =37.5-38.0 kcal/mol (for compounds with 3 to 8 carbon atoms). So, no significant influence of the length of the carbon chain was observed and the carbon-nitrogen bond dissociation remained the main mechanism of decomposition. The C-NO 2 bond dissociation is the main mechanism of decomposition of all tested compounds with ΔG diss ranging from 21.1 kcal/mol to 26.6 kcal/mol.

Besides, no variation in geometric parameters appeared in the various stationary points when increasing the carbon chain from trinitrobutane to trinitrooctane, as summarized in Supporting Information (Tables S2-S5), with only small deviation for the smallest molecules (trinitromethane and trinitroethane) due to steric effects, in accordance with energetic considerations. Indeed, in the case of trinitromethane, the orientation of the nitro groups is only managed by the steric hindrance between nitro groups. All ONCH dihedral angles (corresponding to the rotation of nitro groups) were calculated to be 42.8°. In the trinitroethane structure, these angles became 43.3° due to the steric hindrance of CH 3 instead of H. Then, from trinitropropane, the configuration of the carbon chain involved a different state of symmetry between nitro groups from a C2v like configuration to a non symmetric conformation with dihedral angles between 29.0° and 52.8°. This new conformation is the result of the combined effect of the steric hindrance of the carbon chain and the weak interactions between hydrogen atoms on carbons 2 and 3 and oxygen atoms in nitro groups.

Influence of the presence of double bonds

Then, the presence of double bonds along the carbon chain was investigated. This corresponds to the case of nitroalkene compounds that were considered as a particular class among C-nitro compounds by Korolev et al. 27 . Three compounds were considered by including double bonds between carbons 2 and 3 (denoted 2=3), between carbons 3 and 4 (3=4) or both (2=3=4). Reactions 1 and 3 were computed, energetic profiles were reported in table 4 and geometric structures were summarized in Supporting Information (tables S6-S9). Two different features were highlighted as discussed in the following.

For the 3=4 compound, energetic profiles were very similar than those obtained for 1,1,1trinitrobutane with ΔG diss =22.6 kcal/mol for the breaking of the C-N bond and ΔG TS2 =35.4 kcal/mol for the HONO elimination. For the 2=3 and 2=3=4 compounds, lower dissociation energies were computed with ΔG diss =9.2 and 11.4 kcal/mol, respectively, whereas activation energies of HONO elimination increased at ΔG TS2 =48.4 and 44.0 kcal/mol, respectively.

The difference between both cases was due to the closeness of the double bond to the nitro groups.

In the structures presenting a 2=3 double bond (i.e. in 2=3 and 2=3=4 structures), carbon 2 is in sp2 hybridation state leading to a different geometric configuration compared to 1,1,1-trinitrobutane and 3=4 compound. Notably, in the Rad1 structures, a weak interaction between the hydrogen on carbon 2 and an oxygen atom of a nitro group led to a planar configuration with NCCC dihedral angle close to 0° as shown in Figure 5 for the 2=3 and 2=3=4 structures. However, the angle was about 100° for the 3=4 structure as for the Rad1 radical for 1,1,1-trinitrobutane.

Finally, the order of preference between reaction 1 and 3 was not changed between the molecules that do or do not present double bonds between carbons 2 and 3. On the contrary, the activation energy of the most favourable channel (C-H dissociation) was lowered by the presence of these bonds with ΔG about 10 kcal/mol (for TNBt and 3=4) compared to more than 20 kcal/mol for compounds without 2=3 double bonds (for 2=3 and 2=3=4). However, present study highlighted no influence of double bonds so far as these bonds are not located in the close surrounding of the nitro groups for these two decomposition channels. This observation is consistent with common statements that considered nitroalkenes as a particular family of compounds 26,27 so far as nitro groups are close to double bonds. It has to be noticed that in such configuration, further consideration would be necessary to investigate other specific decomposition process 26 .

Influence of substitution on carbon 2, 3 and 4.

In the following step, substituent effects were investigated. Substitutions on carbons 2, 3 and 4 were considered with six functional groups: methyl (CH 3 ), chlorine (Cl), cyano (CN), amino (NH 2 ), nitro (NO 2 ) and alcohol (OH). This sample of groups was chosen to be representative in the scale of chemical inductive effects as defined by their Taft sigma values σ* that ranged from 0.0 (for CH 3 ) to 4.0 (for NO 2 ) [START_REF] Dean | Lange's Handbook of Chemistry[END_REF] . These coefficients were introduced by Taft in the early 1950s on reactions of aliphatic systems [START_REF] Taft | [END_REF][54][55] and more recently used for aliphatic carboxylic acids and alcohols by Zhang et al. 56 and for aliphatic CH 3 -X compounds by Travnikova et al. 57 . The two main decomposition mechanisms highlighted for 1,1,1-trinitrobutane (reactions 1 and 3) were computed for these different compounds and Gibbs free energy profiles were summarized in table 5. The geometric structures of stationary points are summarized in Supporting Information (Tables S10-S21). An influence of the nature of the substituent was highlighted on carbons 2 and 3. Indeed, for substitution on carbon 2 the dissociation Gibbs free energy of the C-N bond varied from 18. Nevertheless, if an influence was exhibited, no linear correlation was directly observed with the Taft sigma values neither for carbon 2 nor for carbon 3 and neither for reaction 1 or for reaction 3. This can be explained by the fact that the substitutions in carbon 2 and 3 led to steric hindrance issued from the closeness with nitro groups. So, inductive effects as characterized by the σ* values revealed not to entirely govern the influence of substituent on carbons 2 and 3. Indeed, direct interaction between the nitro groups and the substituent would also have an effect, in particular on the carbon 2.

To quantify this effect, the different global reactivity descriptors for the different substituted on carbon 2 reactants are given in the table 6. Charges on C 1 and C 2 and the substituent groups were added up, calculated based on the Atomic Polar Tensors (APT) approach 58 . Among the 6 substituents groups, one is electron-donating group (+I), CH 3 , and five are electron withdrawing groups (-I), NO 2 , NH 2 , OH, Cl and CN. As expected, the partial charge of CH 3 group is slightly more positive than for the hydrogen atom whereas for the other groups the partial charges are negative. Concerning the electron withdrawing groups, the following sequence is found q(OH) < q(Cl) < q(CN) < q(NO 2 ) < q(NH 2 ) meaning that the stronger group electronegativity in our case was found for OH group. Looking at the group electronegativity scale obtained by De Proft et al. 59 at CISD/6-31++G(d,p) level, the same trend was observed exception given for OH and CN. However, these two groups belong, with OCH 3 , CCH, SCH 3 and SiH 3 , to the "problem" groups as defined by De Proft during their correlation study with other group electronegativity scales. The presence of the substituent groups affects not only the partial charges of the  carbon but also the trinitro-substituted carbon with, for this latter, a decrease of 50% as the substituent groups range from H to OH. For the others descriptors scales no clear tendency can be extracted.

However, if both mechanisms were affected by the substitution on carbons 2 and 3, the dissociation of the C-N bond remained preferred compared to the HONO elimination, as illustrated in Figure 6 (for carbon 2).

Concerning the substitution on carbon 4, no significant effect was observed neither on the C-N breaking, with ΔG diss ranging from 21.1 kcal/mol for NH 2 to 21.6 kcal/mol for OH, or the HONO elimination, with ΔG TS2 ranging from 37.5 kcal/mol for NH 2 to 38.2 kcal/mol for NO 2 . So, reaction 1 remained the favored mechanism. This is explained by the fact that inductive effects were of only low influence due to the distance between substituents and nitro groups. Moreover, steric hindrance of substituent did not affect the nitro groups.

Finally, it appeared that only substitution on carbons on the close surrounding of nitro groups (from a geometric point of view) led to significant changes in the energetic profiles of the decomposition channels of 1,1,1-trinitrobutane. The inductive effects of substituents, steric hindrance and their possible interaction with nitro groups led to conformational rearrangement. However, whatever the position of the substituents, the competition between the C-N bond breaking and the HONO elimination was unchanged and the breaking of C-NO 2 bond remained the favored mechanism.

Case of trinitroalkyl aromatic compounds

In previous works, nitroaromatic compounds were studied in view of deriving predictive models to evaluate their impact sensitivity based on their molecular structures 11 . In this work, compounds presenting -R-C(NO 2 ) 3 fragments were highlighted as a particular class of compounds among aromatic nitro compounds. In previous sections, the C-N bond breaking was presented as the primary channel of decomposition of trinitroalkanes with ΔG diss =17.3-23.4 kcal/mol whereas the loss of nitro groups were calculated to be ΔG diss =47.8-63.0 kcal/mol in nitrobenzene derivatives (through C-N dissociation) 17 and ΔG ortho =40.2-43.6 kcal/mol in o-nitrotoluene derivatives (initiated ortho interaction between NO 2 and CH 3 ) 16 . So, aromatic compounds presenting trinitroalkyl groups may be considered as decomposing in an aliphatic type process by C-N bond breaking on the terminal C(NO 2 ) 3 groups. To validate this statement, two additional molecules were studied, 4-(1,1,1-trinitrobutyl)-nitrobenzene (p-NO 2 -Bz) and 2-(1,1,1-trinitrobutyl)-nitrobenzene (o-NO 2 -Bz) (in Figure 7). For these two molecules, the aromatic and aliphatic type mechanisms were computed (see Table 7).

For p-NO 2 -Bz, the loss of the nitro groups linked to the aryl ring by dissociation of the C-N bond was calculated with ΔG diss =59.0 kcal/mol whereas similar dissociation for a carbon of the trinitroalkyl fragment was computed with only ΔG diss =21.9 kcal/mol. So, p-NO 2 -Bz decomposes in an aliphatic type decomposition process.

In the case of o-NO 2 -Bz, an additional aromatic decomposition channel was computed. Similarly to the case of o-nitrotoluene derivatives, the nitro group directly interacts with hydrogen on the first carbon of the alkyl chain in the ortho position. The rate limiting step of this mechanism corresponds to the initial hydrogen transfer from the carbon of the alkyl chain to an oxygen of the nitro group with an activation Gibbs free energy ΔG ortho =41.2 kcal/mol. Otherwise, ΔG diss was calculated at only 21.8 kcal/mol for the breaking of the C-N bond in the C(NO 2 ) 3 terminal fragment. So, the aliphatictype channel remained favo red in the decomposition of o-NO 2 -Bz.

Finally, the decomposition process of these two aromatic nitro compounds was initiated on the trinitroalkyl fragment. This confirmed the specificity of such compounds since their decomposition was evidenced to follow the aliphatic behaviour of the trinitroalkyl fragment. Besides, the activation ΔGs of both the C-N bond breaking and the HONO elimination were in the range of the classical trinitroalkanes considered in the previous section. So, such compounds would have to be considered within the same class of compounds in terms of decomposition initiation process.

If the aim of this paper was not to develop QSPR models, the linear correlation between the activation energy of this mechanism, i.e. the dissociation energy of the C-NO 2 bond on the trinitroalkyl fragment, and the impact sensitivities (log h 50% ) for such class of compounds gathered in previous works 11 was considered. But only low correlation between ΔG diss and logh 50% was observed with 0.58 in R², for the only seven available molecules (as shown in Figure 8). To go further in the finding of a QSPR model, a dedicated study could be realized in the future with a larger dataset.

CONCLUSION

In this study, the decomposition mechanisms of a series of nitro compounds presenting trinitroethyl fragments were studied by theoretical means based on density functional theory calculations. From the five studied decomposition channels, the breaking of the C-N bond (reaction 1) was highlighted as the primary decomposition path with ΔG diss =21.3 kcal/mol. HONO elimination (reaction 3) was the most competitive other channel with ΔG TS2 =37.7 kcal/mol.

The influences of various structural parameters on these reaction channels were also studied.

Considering the length of the carbon chain, the activation energies of both reactions were slightly lowered from trinitromethane to trinitrooctane. Identically, activation energies were influenced by the substitution of the 1,1,1-trinitrobutane in various sites (in particular on carbon 2 and 3) due to inductive and steric effects of the functional groups. Nevertheless, in both cases, no change in the competition between reactions 1 and 3 was highlighted. More influence was observed when considering the presence of double bonds in the carbon chain, explaining the particular character of nitroalkenes, but the C-N breaking remained favoured.

From these studies, the breaking of the C-N bond was evidenced as the primary channel for the decomposition of such trinitroalkanes. Moreover, it was highlighted, based on the analysis of 4-(1,1,1-trinitrobutyl)-nitrobenzene and 2-(1,1,1-trinitrobutyl)-nitrobenzene, that aromatic nitro compounds presenting trinitroalkyl groups are a particular class of aromatic nitro compounds. For these specific compounds, the breaking of the aliphatic C-NO 2 bonds were calculated to be favoured compared to the loss of nitro groups linked to the aromatic cycle. This last observation was in accordance with the identification of this particular family of aromatic nitro compounds in previous QSPR study 11 .

Figure 3 -

 3 Figure 3 -Relative Gibbs free energy diagram for the main decomposition paths of 1,1,1trinitrobutane calculated at PBE0/6-31+G(d,p) level (see Figure 4 for nomenclature of the different stationary points).

Figure 4 -

 4 Figure 4 -Structures of the different stationary points (minima and transition states) of the main decomposition channels for 1,1,1-trinitrobutane presented in Figure 3.

Figure 5 -Figure 6 -Figure 8 -

 568 Figure 5 -Influence on the presence of double bonds on the geometric structures of the radicals issued from the C-N dissociation (Rad1)

  2 kcal/mol for Cl to 23.4 kcal/mol for OH whereas activation energy of the HONO elimination ΔG TS2 changed from 35.5 kcal/mol for CN to 39.4 for OH. For substitution on carbon 3, ΔG diss ranged from 17.3 kcal/mol for Cl to 21.5 kcal/mol for OH and ΔG TS2 ranged from 33.7 for Cl to 37.7 kcal/mol for H. It has to be noticed that ΔG diss associated to the breaking of the C-N bond decreases when adding a chlorine atom on both carbons 2 and 3 (by -3.1 and 4.0 kcal/ mol, respectively). This observation is in agreement with the trend highlighted byKhrapkovskii 31 for a series of nitromethane derivatives that noticed a decrease of 5.8 kcal/mol and 11 kcal/mol from CH 3 NO 2 to CH 2 ClNO 2 and CHCl 2 NO 2 , respectively.

Table 1 -

 1 Bond dissociation energies (kcal/mol) for trinitroalkanes

		experimental	calculated
		Nazin data a	Lebedev data a	∆E [PBE0/6-31+G(d,p)]
	CH 3 -C(NO 2 ) 3	94.6	95.5	96.2
	C 2 H 5 -C(NO 2 ) 3	91.3	93.3	92.1
	C 4 H 9 -C(NO 2 ) 3	89.0	90.9	92.8
	a extracted from ref. 28

Table 2 -

 2 Relative energies (ΔE), enthalpies (ΔH) and Gibbs free energies (ΔG) in kcal/mol for the main decomposition channels of 1,1,1-trinitrobutane at PBE0/6-31+G(d,p) level. The decomposition channels are illustrated in Figure1.

		ΔE (kcal/mol)		ΔH (kcal/mol)		ΔG (kcal/mol)
		NO 2 (a)	NO 2 (b)	NO 2 (a)	NO 2 (b)	NO 2 (a)	NO 2 (b)
	S0	0.0	0.0	0.0	0.0	0.0	0.0
	Reaction 1
	Rad1 + NO 2	35.4	35.4	36.0	36.0	21.3	21.3
	Reaction 2
	TS1	60.0	58.5	60.0	58.5	60.0	58.1
	LM1	-10.7	-12.0	-10.7	-12.0	-11.1	-12.7
	Rad2 + NO	24.2	23.1	24.7	23.7	12.4	10.9
	Reaction 3
	TS2	62.4	38.0	62.4	38.0	62.3	37.7
	Prod1 + HONO	6.1	6.1	6.7	-6.7	-6.9	-6.9
	S0'	N.A.	1.7	N.A.	1.7	N.A.	2.5
	TS3	N.A.	83.2	N.A.	83.2	N.A.	86.0
	Prod2 + HONO	N.A.	0.5	N.A.	1.1	N.A.	-10.7
	Reaction 4
	Rad3 + O	155.8	156.2	156.4	156.8	147.1	147.3
	Reaction 5
	Rad4 + C  (NO 2 ) 3	92.9		93.5		74.7

Table 3 -

 3 Relative Gibbs Free Energies (ΔG) in kcal/mol for the main decomposition channels of 1,1,1trinitrobutane (n C =4) and its analogues with different lengths of carbon chain (n C =1-8) at PBE0/6-31+G(d,p) level

	n C	1	2	3	4	5	6	7	8
	S0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
	Reaction 1
	Rad1 + NO 2	26.6	23.5	21.8	21.3	21.5	21.5	21.2	21.1
	Reaction 3
	TS2	N.A.	40.7	38.0	37.7	37.6	37.6	37.7	37.5
	Prod1 + HONO N.A.	1.0	-7.0	-6.9	-7.1	-7.3	-7.4	-7.6

Table 4 -

 4 Relative Gibbs Free Energies (ΔG) in kcal/mol for the main decomposition channels of 1,1,1trinitrobutane and its analogues presenting double bonds along the carbon chain at PBE0/6-31+G(d,p) level

		1,1,1-trinitrobutane	2=3 a	3=4 a	2=3=4 a
	S0	0.0	0.0	0.0	0.0
	Reaction 1
	Rad1 + NO 2	21.3	9.2	22.6	11.4
	Reaction 3
	TS2	37.7	48.4	35.4	44.0
	Prod1 + HONO	-6.9	4.6	-10.7	2.0
	a i=j: double bonds between carbon atoms labelled i and j, respectively

Table 5 -

 5 Relative Gibbs Free Energies (ΔG) in kcal/mol for the main decomposition channels of 1,1,1trinitrobutane and its analogues by substitutions on carbons 2, 3 and 4 at PBE0/6-31+G(d,p) level

	on carbon 2	H	CH 3	Cl	CN	NH 2	NO 2	OH
	S0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
	Reaction 1
	Rad1 + NO 2	21.3	19.6	18.2	20.4	22.9	18.6	23.4
	Reaction 3
	TS2	37.7	38.2	38.8	35.5	38.1	37.0	39.4
	Prod1 + HONO	-6.9	-12.1	-9.8	-9.3	-26.0	-5.0	-23.1
	on carbon 3	H	CH 3	Cl	CN	NH 2	NO 2	OH
	S0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
	Reaction 1
	Rad1 + NO 2	21.3	18.2	17.3	18.6	19.6	21.1	21.5
	Reaction 3
	TS2	37.7	34.7	33.7	35.0	36.9	36.4	36.1
	Prod1 + HONO	-6.9	-10.5	-9.7	-8.1	-6.1	-5.6	-6.4
	on carbon 4	H	CH 3	Cl	CN	NH 2	NO 2	OH
	S0	0.0	0.0	0.0	0.0	0.0	0.0	0.0
	Reaction 1
	Rad1 + NO 2	21.3	21.5	21.3	21.2	21.1	21.3	21.6
	Reaction 3
	TS2	37.7	37.6	38.0	38.1	37.5	38.2	37.6
	Prod1 + HONO	-6.9	-7.1	-6.4	-5.9	-7.1	-6.2	-7.1
	σ* 52	0.49	0.0	2.96	3.30	0.62	4.0	1.34

Table 6 -

 6 Vertical Ionization Potential (I, in eV), Vertical Electron Affinity (A, in eV), Electronegativity (, in eV), Chemical Hardness (, in eV), Electrophilicity Index (, in eV), charges of the substituent groups (q(group)) and for the carbon atoms C 1 and C 2 (q(C 2 ) and q(C 1 ), respectively).

	Groups	I	A				q(group)	q(C 2 )	q(C 1 )
	H	9.6	3.4	6.5	3.1	6.8	0.04	-0.04	0.28
	CH 3	9.6	3.2	6.4	3.1	6.5	0.07	0	0.25
	Cl	9.3	3.5	6.4	2.9	7.2	-0.24	0.27	0.18
	CN	9.9	3.7	6.8	3.1	7.5	-0.18	0.08	0.22
	NH 2	8.1	3.2	5.7	2.4	6.7	-0.14	0.3	0.17
	NO 2	9.6	3.7	6.7	2.9	7.6	-0.16	0.04	0.18
	OH	9.4	3.5	6.5	2.9	7.2	-0.27	0.42	0.14

Table 7 -

 7 Activation energies in terms of Gibbs free energies (ΔG) in kcal/mol for the main decomposition channels of 4-(1,1,1-trinitrobutyl)-nitrobenzene (p-NO 2 -Bz) and 2-(1,1,1-trinitrobutyl)nitrobenzene (o-NO 2 -Bz) at PBE0/6-31+G(d,p) level

	p-NO 2 -Bz	o-NO 2 -Bz


Yves Le Jan
email: yves.lejan@math.u-psud.frjrosen30@optimum.net

Yves Le

Jan Jay Rosen

Perturbation of the loop measure

come

Introduction

Professor Fukushima's contribution to probabilistic potential theory was seminal. His book on Dirichlet spaces [START_REF] Fukushima | Dirichlet Forms and Symmetric Markov Processes[END_REF] was the first complete exposition in which the functional analytic, potential theoretic, and probabilistic aspects of the theory were considered jointly, as different aspects of the same mathematical object. In particular, transformations of the energy form, such as restriction to an open set, trace on a closed set, change of the equilibrium (or killing) measure, change of reference measure, etc... have probabilistic counterparts. In the second chapter, the possibility of superposing different Dirichlet forms was recognized. The present paper follows this line of research. But we will focus on Markov loops and bridges rather than Markov paths. We will present them in the next section. Let us stress however that the existence of loop or bridge measures seems to require more than the assumption of a Markov process defined up to polar sets, which is the basic assumption for a Markov process associated with a Dirichlet form. The existence of a Green function seems to be required. Our purpose is to compute the variation of the loop measure induced by an infinitesimal variation of the generator. This variation may in particular affect the killing rates or the jumping rates. In the case of symmetric continuous time Markov chains, the question has been addressed in chapter 6 of [START_REF] Jan | Markov paths, loops and fields[END_REF]. The results are formally close to formulas used in conformal field theory for operator insertions (see for example [8]). We try here to extend them to a more general situation, to show in particular that the bridge measures can be derived from the loop measure.

Background on loop measures

We begin by introducing loop measures for Borel right processes (such as Feller processes) on a rather general state space S, which we assume to be locally compact with a countable base. Let X = (Ω, F t , X t , θ t , P x ) be a transient Borel right process [START_REF] Sharpe | General theory of Markov processes[END_REF] with cadlag paths (such as a standard Markov process [START_REF] Blumenthal | Markov Processes and Potential Theory[END_REF]) with state space S and jointly measurable transition densities p t (x, y) with respect to some σ-finite measure m on S. We assume that the potential densities u(x, y) = ∞ 0 p t (x, y) dt are finite off the diagonal, but allow them to be infinite on the diagonal. We do not require that p t (x, y) is symmetric. We assume furthermore that 0 < p t (x, y) < ∞ for all 0 < t < ∞ and x, y ∈ S, and that there exists another Borel right process X in duality with X (Cf [START_REF] Blumenthal | Markov Processes and Potential Theory[END_REF]), relative to the measure m, so that its transition probabilities P t (x, dy) = p t (y, x) m(dy). These conditions allow us to use material on bridge measures in [6]. In particular, for all 0 < t < ∞ and x, y ∈ S, there exists a finite measure Q x,y t on F t -, of total mass p t (x, y), such that

Q x,y t 1 {ζ>s} F s = P x (F s p t-s (X s , y)) , (2.1)
for all F s ∈ F s with s < t. (We use the letter Q for measures which are not necessarily of mass 1, and reserve the letter P for probability measures.) Q x,y t should be thought of as a measure for paths which begin at x and end at y at time t. When normalized, this gives the bridge measure P x,y t of [6]. We use the canonical representation of X in which Ω is the set of cadlag paths ω in S ∆ = S ∪ ∆ with ∆ / ∈ S, and is such that ω(t) = ∆ for all t ≥ ζ = inf{t > 0 |ω(t) = ∆}. Set X t (ω) = ω(t). We define a σ-finite measure µ on (Ω, F) by

F dµ = ∞ 0 1 t Q x,x t (F • k t ) dm(x) dt (2.2)
for all F measurable functions F on Ω. Here k t is the killing operator defined by

k t ω(s) = ω(s) if s < t and k t ω(s) = ∆ if s ≥ t, so that k -1 t F ⊂ F t -.
As usual, if F is a function, we often write µ(F ) for F dµ. µ is σ-finite as any set of loops in a compact set with lifetime bounded away from zero and infinity has finite measure.

We call µ the loop measure of X because, when X has continuous paths, µ is concentrated on the set of continuous loops with a distinguished starting 2 point (since Q x,x t is carried by loops starting at x). Moreover, it is shift invariant. More precisely let ρ u denote the loop rotation defined by

ρ u ω(s) = ω(s + u mod ζ(ω)), if 0 ≤ s < ζ(ω) ∆, otherwise.
Here, for two positive numbers a, b we define a mod b = a -mb as the unique positive integer m such that 0 ≤ a -mb < b . µ is invariant under ρ u , for any u. We let F ρ denote the σ-algebra of F measurable functions F on Ω which are invariant under ρ, that is F • ρ u = F for all u ≥ 0. Loop functionals of interest are mostly F ρ -measurable. Recall that Poisson processes of intensity µ appear naturally as produced by loop erasure in the construction of random spanning trees through the Wilson algorithm (see chapter 8 in [START_REF] Jan | Markov paths, loops and fields[END_REF] ). Although the definition of µ in (2.2), especially the 1 t , may look forbidding, µ often has a nice form when applied to specific functions in F ρ . A particular function in F ρ is given by

φ(f ) = ∞ 0 f (X t ) dt, (2.3)
where f is any measurable function on S. If f j , j = 1, . . . , k ≥ 2 are nonnegative functions on S, then by [11, Lemma 2.1]

µ   k j=1 φ(f j )   (2.4) = π∈P ⊙ k u(y π(1) , y π(2) ) • • • u(y π(k-1) , y π(k) )u(y π(k) , y π(1) ) k j=1
f j (y j ) dy j where P ⊙ k denotes the set of permutations of [1, k] on the circle. (For example, (1, 3, 2), (2, 1, 3) and (3, 2, 1) are considered to be one permutation π ∈ P ⊙ 3 .) We note however that in general when u is infinite on the diagonal

µ (φ(f j )) = ∞.
For k ≥ 2, the integral (2.4) can be finite if the f i satisfy certain integrability conditions: see the beginning of section 3. Consider more generally the multiple integral

π∈T ⊙ k 0≤r 1 ≤•••≤r k ≤t f π(1) (X r 1 ) • • • f π(k) (X r k ).
(2.5)

where T ⊙ k denotes the set of translations π of [1, k] which are cyclic mod k, that is, for some i, π(j) = j + i, mod k, for all j = 1, . . . , k.

Finite sums of multiple integrals such as these form an algebra (see exercise 11, p.21 in [START_REF] Jan | Markov paths, loops and fields[END_REF]) which generates F ρ , [3].

Finally, let a(x) be a bounded, strictly positive function on S. Define the time changed process Y t = X τ (t) where τ (t) = t 0 a(Y s ) ds is the inverse of the CAF A t = t 0 1/a(X s ) ds. It satisfies the duality assumption relative to the measure a • m. It then follows as in [5, section 7.3] that if u X (x, y), u Y (x, y) denote the potential densities of X, Y respectively with respect to m, a • m respectively, then u Y (x, y) = u X (x, y)/a(y).

(2.6)

It follows that µ Y is the image of µ X by the time change.

3 Multiple CAF's and perturbation of loop measures

Let M(S) be the set of finite signed Radon measures on B(S). We say that a norm • on M(S) is a proper norm with respect to a kernel u if for all n ≥ 2 and ν 1 , . . . , ν n in M(S)

n j=1 u(x j , x j+1 ) n i=1 dν j (x i ) ≤ C n n j=1 ν j , (3.1)
(with x n+1 = x 1 ) for some universal constant C < ∞. In Section 6 of [START_REF] Le Jan | Permanental fields, loop soups and continuous additive functionals[END_REF] we present several examples of proper norms which depend upon various hypotheses about the kernel u.

In particular, the following norm is related to the square root of the generator of X, which defines the Dirichlet space in the m-symmetric case:

ν w := w(x, y)w(y, z) dν(y) 2 dm(x) dm(z) 1/2 , ( 3.2)
where

w(x, y) = ∞ 0 p s (x, y) √ πs ds. (3.3)
To see that ν w is a proper norm we first note that u(x, z) = w(x, y)w(y, z) dm(y) (3.4) (It is interesting to note that w is the potential density of the process X Tt where T t is the stable subordinator of index 1/2. In operator notation (3.4) says that W 2 = U where W and U are operators with kernels w and u respectively.) Using (3.4)

n j=1 u(z j , z j+1 ) = n j=1 w(z j , λ j )w(λ j , z j+1 ) dm(λ j ) (3.5) = n j=1 w(z j , λ j )w(λ j-1 , z j ) dm(λ j )
in which z n+1 = z 1 and λ 0 = λ n . It follows from this that

n j=1 u(z j , z j+1 ) n j=1 dν j (z j ) (3.6) = n j=1 w(z j , λ j )w(λ j-1 , z j ) dν j (z j ) n j=1 dm(λ j ) ≤ n j=1 w(z j , s)w(t, z j ) dν j (z j ) 2 dm(s) dm(t) 1/2 ≤ n j=1 ν j w ,
where, for the first inequality, we use repeatedly the Cauchy-Schwarz inequality.

Lastly, set M ν (x, z) = w(x, y)w(y, z) dν(y).

(3.7)

Since ν w is the L 2 norm of M ν , and M ν+ν ′ = M ν + M ν ′ , we see that ν w is a norm. (This can also be viewed as the Hilbert-Schmidt norm of the operator defined by the kernel M ν ).

We denote by R + the set of positive bounded Revuz measures ν on S that are associated with X. This is explained in detail in Section 2.1 of [START_REF] Le Jan | Permanental fields, loop soups and continuous additive functionals[END_REF]. We use L ν t to denote the CAF with Revuz measure ν. Let • be a proper norm on M(S) with respect to the kernel u. Set

M + • = {positive ν ∈ M(S) | ν < ∞}, (3.8) and R + • = R + ∩ M + • . (3.9) Let M • and R • denote the set of measures of the form ν = ν 1 -ν 2 with ν 1 , ν 2 ∈ M + • or R +
• respectively. We often omit saying that both R • and • depend on the kernel u.

Let • be a proper norm for u. For

ν j ∈ R • , j = 1, . . . , k, let M ν 1 ,...,ν k t = π∈T ⊙ k 0≤r 1 ≤•••≤r k ≤t dL ν π(1) r 1 • • • dL ν π(k) r k . (3.10)
We refer to M ν 1 ,...,ν k t as a multiple CAF. Let

Q x,y (F ) = ∞ 0 Q x,y t (F • k t ) dt. (3.11)
We have the following analogue of [9, Proposition 5] and [10, Lemma 2.1].

Lemma 3.1 For any measures

ν j ∈ R + • , j = 1, . . . , k ≥ 2, µ(M ν 1 ,...,ν k ∞ ) = u(y 1 , y 2 ) • • • u(y k-1 , y k )u(y k , y 1 ) k j=1 dν j (y j ), (3.12)
Q x,y (M ν 1 ,...,ν k ∞ ) (3.13) = π∈T ⊙ k u(x, y 1 )u(y 1 , y 2 ) • • • u(y k-1 , y k )u(y k , y) k j=1
dν π(j) (y j ),

and if ν ∈ R + • µ(M ν 1 ,...,ν k ∞ L ν ∞ ) (3.14) = k i=1   i-1 j=1 u(y j , y j+1 )   u(y i , x)u(x, y i+1 )   k j=i+1 u(y j , y j+1 )   k j=1 dν j (y j ) dν(x) = Q x,x (M ν 1 ,...,ν k ∞ ) dν(x),
with y k+1 = y 1 .

The proof of (3.12) follows that of [10, Lemma 2.1], noticing that the crucial step [10, (2.23)-(2.28)] used the fact that the set of permutations of [1, k] is invariant under translation mod k. Since T ⊙ k is invariant under translation mod k, the same proof will work here. The proof of (3.13), which is much easier, follows that of [START_REF] Le Jan | Permanental fields, loop soups and continuous additive functionals[END_REF]Lemma 4.2]. The first equality in (3.14) follows from (3.12) and the fact that

M ν 1 ,...,ν k ∞ L ν ∞ = k i=1 M ν 1 ,...,ν i-1 ,ν,ν i ,...,ν k ∞ , (3.15) which is easily verified using L ν ∞ = ∞ 0 dL ν t .
The second equality in (3.14) then follows by comparing with (3.13).

Let now X (ǫ) , ǫ ≥ 0, X (0) = X, be a family of Markov processes with potential densities u (ǫ) (x, y), and let µ (ǫ) denote the loop measure for X (ǫ) . Assume that we can use the same proper norm • for all u (ǫ) . Let

u ′ (0) (x, y) = du (ǫ) (x, y) dǫ | ǫ=0 , (3.16)
and assume that • is also a proper norm for u ′ (0) . Then using the last Lemma we have formally, that is, assuming we can justify interchanging derivative and integral in the second equality,

d dǫ µ (ǫ) (M ν 1 ,...,ν k ∞ ) ǫ=0 (3.17) = d dǫ u (ǫ) (y 1 , y 2 ) • • • u (ǫ) (y k-1 , y k )u (ǫ) (y k , y 1 ) k j=1 dν j (y j ) ǫ=0 = k i=1 u(y 1 , y 2 ) • • • u ′ (0) (y i , y i+1 ) • • • u(y k , y 1 ) k j=1 dν j (y j ) = π∈T ⊙ k u(y 1 , y 2 ) • • • u(y k-1 , y k )u ′ (0) (y k , y 1 ) k j=1 dν π(j) (y j )
where we have set y k+1 = y 1 . Assume now that for some distribution F on S × S we have

u ′ (0) (y k , y 1 ) = S×S u(y k , x)F (x, y) u(y, y 1 ) dm(x) dm(y). (3.18)
Let A • denote the space spanned by the multiple CAF's with ν j ∈ R • . Note it is an algebra. Then comparing (3.17) and (3.13) we would obtain

dµ (ǫ) (A) dǫ | ǫ=0 = S×S F (x, y)Q y,x (A) dm(x) dm(y), (3.19) for all A ∈ A • .
In the following sections we present specific examples where this heuristic approach is made rigorous.

Perturbation of Lévy processes

Let X be a transient Lévy process in R d with characteristic exponent ψ so that, as distributions

u(x, y) = e iλ(y-x) ψ(λ) dλ. (4.1)
In [START_REF] Le Jan | Permanental fields, loop soups and continuous additive functionals[END_REF] we showed that • ψ,2 is a proper norm for u where

ν 2 ψ,2 = 1 |ψ| * 1 |ψ| (λ) | ν(λ)| 2 dλ. (4.2)
Let κ be a Lévy characteristic exponent, so that the same is true for ψ +ǫκ, and let X (ǫ) be the Lévy process with characteristic exponent ψ + ǫκ. We let u (ǫ) (x, y) denote the potential of X (ǫ) so that, as distributions

u (ǫ) (x, y) = e iλ(y-x) ψ(λ) + ǫκ(λ) dλ. (4.3) If we assume that |κ(λ)| ≤ C|ψ(λ)| (4.4)
for some C < ∞, then for ǫ > 0 sufficiently small

ν ψ+ǫκ,2 ≤ C ′ ν ψ,2 , (4.5)
for some C ′ < ∞. Thus • ψ,2 is a proper norm for u (ǫ) .

Let F be the distribution given by F (x, y) = e iλ(y-x) κ(λ) dλ, (4.6) and let Q λ 1 ,λ 2 (A) denote the Fourier transform of Q x,y (A) in x, y.

Theorem 4.1 If (4.4) holds, then

dµ (ǫ) (A) dǫ | ǫ=0 = - R d ×R d Q y,x (A)F (x, y) dm(x) dm(y) (4.7) = -Q λ,-λ (A) κ(λ) dλ. for all A ∈ A • ψ,2 .
Proof of Theorem 4.1: It suffices to show that for any ν 1 , . . . ,

ν k ∈ R + • ψ,2 d dǫ k j=1 u (ǫ) (y j , y j+1 ) k j=1 dν j (y j ) ǫ=0 (4.8) = - k i=1   i-1 j=1 u(y j , y j+1 )   R d ×R d u(y i , x) F (x, y) u(y, y i+1 ) dm(x) dm(y)   k j=i+1 u(y j , y j+1 )   k j=1 dν j (y j ),
with y k+1 = y 1 . Using (4.3) we see that

I(ǫ) =: k j=1 u (ǫ) (y j , y j+1 ) k j=1 dν j (y j ) (4.9) = k j=1 e i(y j+1 -y j ) • λ j 1 ψ(λ j ) + ǫκ(λ j ) dλ j dν j (y j ) =   k j=1 e -i(λ j -λ j-1 ) • y j dν j (y j )   k j=1 1 ψ(λ j ) + ǫκ(λ j ) dλ j
where λ 0 = λ n . We take the Fourier transforms of the {ν j } to see that

k j=1 u (ǫ) (y j , y j+1 ) k j=1 dν j (y j ) (4.10) = k j=1 νj (λ j -λ j-1 ) 1 ψ(λ j ) + ǫκ(λ j ) dλ j .
We have

1 ψ(λ j ) + ǫκ(λ j ) = 1 ψ(λ j ) -ǫ κ(λ j ) ψ 2 (λ j ) + ǫ 2 κ 2 (λ j ) ψ 2 (λ j )(ψ(λ j ) + ǫκ(λ j )) . (4.11)
Substituting this into (4.11) we can write the result as

I(ǫ) = I(0) -ǫJ + K(ǫ), (4.12)
where

J = k i=1     k j=1 j =i 1 ψ(λ j )     κ(λ i ) ψ 2 (λ i ) k j=1 νj (λ j -λ j-1 ) dλ j , (4.13)
and K(ǫ) is the sum of all remaining terms. We now show that J is precisely the right hand side of (4.8), and that

|K(ǫ)| = O(ǫ 2 ), (4.14)
which will complete the proof of our Proposition.

For the first point, using the relation between Fourier transforms and convolutions we have

R d ×R d u(y i , x) F (x, y) u(y, y i+1 ) dm(x) dm(y) (4.15) = e iλ i (y i+1 -y i ) κ(λ i ) ψ 2 (λ i ) dλ i .
Using this in the right hand side of (4.8) and proceeding as in (4.9)-(4.10) we indeed obtain J. As for (4.14), K(ǫ) is the sum many terms each of which has a factor of ǫ m for some m ≥ 2. We need only show that the corresponding integrals are bounded uniformly in ǫ. For example, consider the term which arises when using the last term in (4.11) for all j. This term is ǫ 2k times

K(ǫ) =: k j=1 νj (λ j -λ j-1 ) κ 2 (λ j ) ψ 2 (λ j )(ψ(λ j ) + ǫκ(λ j )) dλ j . (4.16)
By our assumption (4.4), for sufficiently small ǫ

| K(ǫ)| ≤ C ′ k j=1 |ν j (λ j -λ j-1 )| 1 |ψ(λ j )| dλ j (4.17) ≤ C ′′ k j=1 ν j ψ,2 ,
by repeated use of the Cauchy-Schwarz inequality, as in our proof in [START_REF] Le Jan | Permanental fields, loop soups and continuous additive functionals[END_REF] that • ψ,2 is a proper norm for u.

Perturbation by multiplicative functionals

Let m t be a continuous decreasing multiplicative functional of X, with m t ≤ 1 for all t and m ζ = 0. By [12, Theorem 61.5], there is a right process X t with transition semigroup

P t f (x) = P x (f (X t )m t ) = Q x,y t (m t )f (y) dm(y), (5.1)
where Q x,y t are the bridge measures (2.1) for our original process X and the second equality is [6, (2.8)]. Thus X has transition densities

p t (x, y) =: Q x,y t (m t ), (5.2)
and using [6, (2.8)] once more we can verify that these satisfy the Chapman-Kolmogorov equations. It follows from the construction in [12, Theorem 61.5] that if X has cadlag paths so will X. Let r t (ω)(s) = ω(t -s) -, 0 ≤ s ≤ t be the time reversal mapping. If we set m t = m t • r t , then it is easy to check that m t is a multiplicative functional as above, see [2, p. 359]. If X is the dual process for X as described in Section 2, with bridge measures Q x,y t then as above there exists a process Y with transition densities Q x,y t ( m t ). By [6, Corollary 1]

Q x,y t ( m t ) = Q x,y t (m t • r t ) = Q y,x t (m t ), (5.3)
which shows that Y is dual to X. We now show that if Q x,y t are the bridge measures for X, then

Q x,y t (F ) = Q x,y t (F m t ), for F ∈ F s , s < t. (5.4)
To see this, using the fact that m t is continuous and decreasing, we have for

F ∈ F s , s < t Q x,y t (F ) = P x F Q Xs,y t-s (m t-s ) (5.5) = lim t * ↑t P x F Q Xs,y t-s (m t * -s ) = lim t * ↑t P x F P Xs (m t * -s p t-t * (X t * -s , y)) = lim t * ↑t P x F m s P Xs (m t * -s p t-t * (X t * -s , y)) = lim t * ↑t P x (F m s m t * -s • θ s p t-t * (X t * , y)) = lim t * ↑t P x (F m t * p t-t * (X t * , y)) = lim t * ↑t Q x,y t (F m t * ) = Q x,y t (F m t ),
which proves (5.4). If A t is a CAF, then m t = e -At is a continuous decreasing multiplicative functional of X. Let X (ǫ) denote the Markov process X with m t = e -ǫL ν t . If µ (ǫ) denotes the loop measure for X (ǫ) , and µ the loop measure for X, it now follows from (2.2) and (5.4) that

µ (ǫ) (A) = µ Ae -ǫL ν ∞ .
(5.6)

Theorem 5.1 If ν ∈ R + • , for some proper norm • , then dµ (ǫ) (A) dǫ | ǫ=0 = -µ(L ν ∞ A) = - S Q x,x (A) dν(x) (5.7)
for all A ∈ A • .

Proof: It suffices to prove this for A of the form

M ν 1 ,...,ν k ∞ with ν j ∈ R + • , 1 ≤ j ≤ k. Since 0 ≤ e -x -1 + x ≤ x 2 /2 for x ≥ 0, it follows from (5.6) that |µ (ǫ) (M ν 1 ,...,ν k ∞ ) -µ(M ν 1 ,...,ν k ∞ ) -ǫµ(L ν ∞ M ν 1 ,...,ν k ∞ )| ≤ ǫ 2 µ (L ν ∞ ) 2 M ν 1 ,...,ν k ∞ . (5.8) µ (L ν ∞ ) 2
A is bounded by our assumption about the proper norm • , so the first equality in (5.7) follows. The second equality is (3.14).

It follows from (5.2) that X (ǫ) has potential densities

u (ǫ) (x, y) = ∞ 0 Q x,y
t (e -ǫL ν t ) dt ≤ u(x, y).

(5.9)

Then, using [6, Lemma 1]

u ′ (0) (x, y) = - ∞ 0 Q x,y t (L ν t ) dt (5.10) = - ∞ 0 t 0 p s (x, z)p t-s (z, y) dν(z) ds dt = -u(x, z)u(z, y) dν(z).
In view of this, Theorem 5.1 is another example of our heuristic formula (3.19), where the distribution F on S × S of (3.18) is δ(x -y) dν(x).

For use in the next section we now recast Theorem 5.1. For

ν j ∈ R + • , 1 ≤ j ≤ k let I(ǫ) = k j=1 u (ǫ) (y j , y j+1 ) k j=1
dν j (y j ).

(5.11)

Since k j=1 u (ǫ) (y j , y j+1 ) k j=1 dν j (y j ) = µ (ǫ) (M ν 1 ,...,ν k ∞
) by (3.12), it follows from Theorem 5. In the next section we will also use the following observation. If m t = e -ǫ t 0 c(s) ds then m t = m t • r t = m t . Hence if X is the dual process for X as described in Section 2, it follows from (5.3) that Y , the dual process to X (ǫ) , is the process ( X) (ǫ) obtained by perturbing X by the same multiplicative functional m t . In particular

u (ǫ) (x, y) = u (ǫ) (y, x).
(5.13)

6 Perturbation by addition of jumps Let j(x, y) be a nonnegative m ⊗ m-integrable function on S × S. We will assume that c(x) =: j(x, y)m(dy) = j(y, x)m(dy). (6.1) c(x) is integrable and we will assume moreover that it is bounded and strictly positive. Then G(x, dy) =: 1 c(x) j(x, y)m(dy) (6.2) is a probability kernel on S × B(S). c(x) will govern the rate at which we will add jumps to the process, which may depend on the position x of the process, and G(x, dy) will describe the distribution of the jumps from position x.

In more detail, define the CAF

A t = t 0 c(X s ) ds, (6.3)
and let τ t be the right continuous inverse of A t . Let λ be an independent mean 1 exponential. We define a new process Y t to be equal to X t for t < τ λ , and then re-birthed at a random point independent of λ, distributed according to G(X τ λ , dy), with this process being iterated. We use U c,G to denote the potential operator of Y . Let

ν u 2 ,∞ := |ν|(S) ∨ sup x u 2 (x, y) + u 2 (y, x) d|ν|(y), (6.4)
where |ν| is the total variation of the measure ν. This is a proper norm for u, see [10, (3.25)].

We use µ (ǫ) to denote the loop measure associated to Y , where we have replaced c by ǫc. Then µ (ǫ) is well defined for ǫ small enough and

dµ (ǫ) (A) dǫ | ǫ=0 = S×S (Q y,x (A) -Q x,x (A)) c(x)G(x, dy) dm(x), (6.6) for all A ∈ A • u 2 ,∞ .
Before proving this theorem, we first show that U ǫc,G has densities u ǫc,G (x, y) for ǫ sufficiently small. Note first that by (6.1)-(6.2) c(z)G(z, dy) dm(z) = j(z, y) dm(y) dm(z) = c(y) dm(y), (6.7)

and since we assumed that c is bounded it follows from (6.5) that sup x c(z)G(z, dy)u(y, x) dm(z) = sup

x c(y)u(y, x) dm(y) < ∞. (6.8)

Let λ 1 , λ 2 , . . . be a sequence of independent mean 1 exponentials, and set T j = j i=1 λ j . Using the fact that τ t+u = τ t + τ u • θ τt we see that

W c,G,n f (x) =: P x τ (T n+1 ) τ (Tn)
f (Y t ) dt (6.9)

= P x τ (λ n+1 ) 0 f (Y t ) dt • θ τ (Tn) = P x P Y τ (Tn) τ (λ n+1 ) 0 f (X t ) dt = P x G(Y - τ (Tn) , dz)P z τ (λ n+1 ) 0 f (X t ) dt .
We have

P x λ τ λ 0 f (X t ) dt (6.10) = P x λ ∞ 0 1 {t<τ λ } f (X t ) dt = P x λ ∞ 0 1 {At<λ} f (X t ) dt = P x ∞ 0 e -At f (X t ) dt .
Hence setting

V c f (x) = P x ∞ 0 e -At f (X t ) dt , (6.11)
and writing Gh(x) = S G(z, dy)h(y) for any nonnegative or bounded function h, we have shown that

W c,G,n f (x) = P x G(Y - τ (Tn) , dz)V c f (z) (6.12) = P x GV c f (Y - τ (Tn) ) .
Using once again the fact that τ t+u = τ t +τ u •θ τt and the Markov property, we see that for any h

P x h(Y - τ (Tn) ) = P x h(Y - τ (λn) ) • θ τ (T n-1 ) (6.13) = P x P Y τ (T n-1 ) h(Y - τ (λn) ) = P x G(Y - τ (T n-1 ) , dz)P z h(X - τ (λn) ) .
Using the change of variables formula, [4, Chapter 6, (55.1)], and the fact that X - t = X t for a.e. t, we see that

P z λ h(X - τ λ ) = P z λ (h(X τ λ )) = E x ∞ 0 e -t h (X τt ) dt (6.14) = P z A ζ 0 e -t h (X τt ) dt = P z ∞ 0 e -As h (X s ) dA s = P z ∞ 0 e -As h (X s ) c(X s ) ds = V c (ch)(z).
Thus we can write (6.13) as

P x h(Y - τ (Tn) ) = P x GV c c h(Y - τ (T n-1 ) ) . (6.15)
Iterating this we obtain

P x h(Y - τ (Tn) ) = P x GV c c h(Y - τ (T n-1 ) ) (6.16) = P x (GV c c) 2 h(Y - τ (T n-2 ) ) = • • • = P x (GV c c) n-1 h(Y - τ (λ) ) = V c c(GV c c) n-1 h(x),
where the last step used (6.14). Applying this to (6.12) we have that

W c,G,n f (x) = V c c(GV c c) n-1 GV c f (x) = V c (cGV c ) n f (x).
(6.17) It follows from (5.9) that V c has a density which we write as v c (x, y), and therefore W c,G,n has the density

w c,G,n (x, y) = v c (x, z 1 )G(z 1 , dz 2 )v c (z 2 , z 3 ) • • • (6.18) • • • G(z 2n-3 , dz 2n-2 )v c (z 2n-2 , z 2n-1 )G(z 2n-1 , dz 2n )v c (z 2n , y) n j=1 c(z 2j-1 ) dm(z 2j-1 )
.

By (6.5) it follows from (5.9) that sup

z v c (z, y)dm(y) ≤ M, (6.19)
and thus by (6.18) that

sup z w c,G,n (z, y)dm(y) ≤ C n M n+1 . (6.20)
Replacing c by ǫc we have shown that for ǫ sufficiently small

U ǫc,G f (x) = ∞ n=0 ǫ n V ǫc (cGV ǫc ) n f (x) (6.21)
for all bounded measurable f . Hence U ǫc,G has a density

u ǫc,G (x, y) = ∞ n=0 ǫ n v ǫc (x, z 1 )G(z 1 , dz 2 )v ǫc (z 2 , z 3 ) • • • (6.22) • • • G(z 2n-3 , dz 2n-2 )v ǫc (z 2n-2 , z 2n-1 )G(z 2n-1 , dz 2n )v ǫc (z 2n , y) n j=1 c(z 2j-1 ) dm(z 2j-1 ), with sup z u ǫc,G (z, y)dm(y) < ∞. (6.23)
We can also write this as

u ǫc,G (x, y) = ∞ n=0 ǫ n v ǫc (x, z 1 )j(z 1 , z 2 )v ǫc (z 2 , z 3 ) • • • (6.24) • • • j(z 2n-3 , z 2n-2 )v ǫc (z 2n-2 , z 2n-1 )j(z 2n-1 , z 2n )v ǫc (z 2n , y) 2n j=1 dm(z j ).
Similar expansions can be given for the semigroup which has therefore a density:

p ǫc,G,t (x, y) = ∞ n=0 ǫ n 0≤t 1 ≤•••≤tn≤t q t 1 (x, z 1 )c(z 1 )G(z 1 , dz 2 )q t 2 -t 1 (z 2 , z 3 ) • • • • • • c(z 2n-1 )G(z 2n-1 , dz 2n )q t-tn (z 2n , y) n j=1
dm(z 2j-1 ) dt j , (6.25) where q t denotes the kernel of the semigroup associated with the process killed at rate ǫc.

To see this let

W c,G,n,t f (x) =: P x (f (Y t ) ; τ (T n ) < t < τ (T n+1 )) . (6.26)
Using the fact that τ t+u = τ t + τ u • θ τt we have

= P x f (Y t ) ; τ (T n ) < t < τ (T n+1 ) τ (T n ) (6.27) = P x f Y t-τ (Tn) ; t -τ (T n ) < τ (λ n+1 ) • θ τ (Tn) τ (T n ) = P x P Y τ (Tn) f X t-τ (Tn) ; t -τ (T n ) < τ (λ n+1 ) τ (T n ) = P x G(Y - τ (Tn) , dz) P z λ n+1 f X t-τ (Tn) ; t -τ (T n ) < τ (λ n+1 ) τ (T n ) .
Conditional on τ (T n ) we have

P x λ f X t-τ (Tn) ; t -τ (T n ) < τ (λ) (6.28) = P x λ f X t-τ (Tn) ; A t-τ (Tn) < λ = P x e -A t-τ (Tn) f X t-τ (Tn) .
Hence setting P c,t f (x) = P x e -At f (X t ) , (6.29)

we have shown that

W c,G,n,t f (x) = P x G(Y - τ (Tn) , dz)P c,t-τ (Tn) f (z) (6.30) = P x GP c,t-τ (Tn) f (Y - τ (Tn) ) .
Using once again the fact that τ t+u = τ t + τ u • θ τt and the strong Markov property, we see that for any h

P x h(Y - τ (Tn) , t -τ (T n )) τ (T n-1 ) (6.31) = P x h(Y - τ (λn) , t -τ (T n-1 ) -τ (λ n )) • θ τ (T n-1 ) τ (T n-1 ) = P x P Y τ (T n-1 ) h(Y - τ (λn) , t -τ (T n-1 ) -τ (λ n )) τ (T n-1 ) = P x G(Y - τ (T n-1 ) , dz)P z h(X - τ (λn) , t -τ (T n-1 ) -τ (λ n )) τ (T n-1 ) .
Using the change of variables formula, [4, Chapter 6, (55.1)], and the fact that X - t = X t for a.e. t, we see that conditionally on τ (T n-1 )

P z λn h(X - τ λn , t -τ (T n-1 ) -τ (λ n )) (6.32) = P z λn h(X τ λn , t -τ (T n-1 ) -τ (λ n )) = E x ∞ 0 e -r h (X τr , t -τ (T n-1 ) -τ r ) dr = P z A ζ 0 e -r h (X τr , t -τ (T n-1 ) -τ r ) dr = P z ∞ 0 e -As h (X s , t -τ (T n-1 ) -s) dA s = ∞ 0 P z e -As h (X s , t -τ (T n-1 ) -s) c(X s ) ds = ∞ 0 P c,s (ch(•, t -s -τ (T n-1 )))(z) ds.
Combining this with (6.31) we have

P x h(Y - τ (Tn) , t -τ (T n )) (6.33) = P x GP c,sn ch(Y - τ (T n-1 ) , t -s n -τ (T n-1 )) ds n .
Iterating this we obtain, conditionally on τ (T n-1 ), then τ (T n-2 ), • • •

P x h(Y - τ (Tn) ), t -τ (T n )) (6.34) = P x GP c,sn ch(Y - τ (T n-1 ) , t -s n -τ (T n-1 )) ds n = P x GP c,s n-1 cGP c,sn ch(Y - τ (T n-2 ) , t -s n -s n-1 -τ (T n-2 )) ds n-1 ds n = • • • = P x   GP c,s 2 c • • • GP c,sn ch(Y - τ (λ) , t - n j=2 s j -τ (λ))   n j=2 ds j = P c,s 1 cGP c,s 2 c • • • GP c,sn ch(x, t - n j=1 s j ) n j=1 ds j
where the last step used (6.32). Applying this to (6.30) we have that

W c,G,n,t f (x) = P c,s 1 cGP c,s 2 c • • • GP c,sn cGP c,t-n j=1 s j f (x) n j=1
ds j . (6.35) (6.25) then follows as in the proof of (6.22). Let Ĝ(x, dy) = 1 c(x) j(y, x)m(dy). (6.36) By (6.1), Ĝ(x, dy) is a probability kernel on S × B(S). We now add jumps to the dual process X, where c(x) will again govern the rate of jumps but we use Ĝ(x, dy) to describe the distribution of the jumps from position x. Performing the same calculation as before, but with the dual process and using (6.24) and (5.13), we see that the two processes obtained by adding jumps have dual potential kernels: ûǫc, Ĝ(x, y) = u ǫc,G (y, x) (6.37)

The same will be true for the associated resolvents and semigroups. The duality assumptions are verified and we can therefore define a loop measure associated with these processes. We use µ (ǫ) to denote the loop measure associated to Y , where we have replaced c by ǫc.

The next Lemma is needed for the proof of Theorem 6.1.

Lemma 6.1 Assume (6.5) (which implies (6.8)). Then for any positive measure ν and ǫ sufficiently small.

and ν u 2 ǫc,G ,∞ ≤ C ν u 2 ,∞ , (6.39)
Proof of Lemma 6.1: (6.38) follows immediately from (6.22) and (5.9). It also follows from (6.22) that for a positive measure ν

u 2 ǫc,G (x, y) dν(y) (6.40) = ∞ m,n=0 ǫ m+n V ǫc (cGV ǫc ) m (x, y)V ǫc (cGV ǫc ) n (x, y) dν(y) = ∞ m,n=0 ǫ m+n V ǫc (cGV ǫc ) m-1 cG(x, dz 1 )V ǫc (cGV ǫc ) n-1 cG(x, dz 2 ) v ǫc (z 1 , y)v ǫc (z 2 , y) dν(y) .
Hence for ǫ small enough sup x u 2 ǫc,G (x, y) dν(y) (6.41)

≤ ∞ m,n=0 ǫ m+n sup x V ǫc (cGV ǫc ) m-1 cG(x, dz 1 )V ǫc (cGV ǫc ) n-1 cG(x, dz 2 ) sup z 1 ,z 2 v ǫc (z 1 , y)v ǫc (z 2 , y) dν(y) ≤ C ν u 2 ,∞ . Similarly u 2 ǫc,G (y, x) dν(y) (6.42) = ∞ m,n=0 ǫ m+n V ǫc (cGV ǫc ) m (y, x)V ǫc (cGV ǫc ) n (y, x) dν(y) = ∞ m,n=0 ǫ m+n (cGV ǫc ) m (z 1 , x) (cGV ǫc ) n (z 2 , x) dm(z 1 ) dm(z 2 )
v ǫc (y, z 1 )v ǫc (y, z 2 ) dν(y) .

Using (6.8) it follows that for ǫ small enough sup x u 2 ǫc,G (x, y) dν(y) (6.43)

≤ ∞ m,n=0 ǫ m+n sup x (cGV ǫc ) m (z 1 , x) (cGV ǫc ) n (z 2 , x) dm(z 1 ) dm(z 2 ) sup z 1 ,z 2 v ǫc (y, z 1 )v ǫc (y, z 2 ) dν(y) ≤ C ν u 2 ,∞ .
It follows from [10, Lemma 3.3] that ν u 2 ,∞ is a proper norm for u ǫc,G,n . Theorem 6.1 follows from the next Lemma.

Lemma 6.2 Under the assumptions of Lemma 6.1, for any ν

1 , . . . , ν k ∈ R + • u 2 ,∞ d dǫ k j=1 u ǫc,G (y j , y j+1 ) k j=1 dν j (y j ) ǫ=0 (6.44) = k i=1   i-1 j=1 u(y j , y j+1 )   -u(y i , x)u(x, y i+1 )c(x) dm(x) + S×S u(y i , z 1 )c(z 1 )G(z 1 , dz 2 )u(z 2 , y i+1 ) dm(z 1 )   k j=i+1 u(y j , y j+1 )   k j=1 dν j (y j ),
with y k+1 = y 1 .

Proof of Lemma 6.2: Set

I(ǫ) = k j=1 u ǫc,G (y j , y j+1 ) k j=1
dν j (y j ) (6.45) with y k+1 = y 1 .

We can write (6.22) as

u ǫc,G (x, y) (6.46) = v ǫc (x, y) + ǫ v ǫc (x, z 1 )c(z 1 )G(z 1 , dz 2 )v ǫc (z 2 , y) dm(y 1 ) +ǫ 2 v ǫc (x, z 1 )c(z 1 )G(z 1 , dz 2 )v ǫc (z 2 , y)c(z 2 )G(z 2 , dz 3 )u ǫc,G (z 3 , y) dm(y 1 ) dm(y 3 ) = v ǫc (x, y) + ǫV ǫc cGV ǫc (x, y) + ǫ 2 V ǫc cGV ǫc cGU ǫc,G (x, y),
with operator notation. We substitute this in (6.45) and collect terms to obtain

I(ǫ) = I(ǫ) + ǫ k i=1 J i (ǫ) + K(ǫ), (6.47)
where

I(ǫ) = k j=1 v ǫc (y j , y j+1 ) k j=1
dν j (y j ) (6.48) However every term in K(ǫ) comes with a pre-factor of ǫ m for some m ≥ 2, so we need only bound the integrals uniformly in ǫ. For this we will use Lemma 6.1. We illustrate this with the most complicated term, which has the pre-factor ǫ 2k : Then by Lemma 6.1, and the fact that v ǫc ≤ u cGV ǫc cGU ǫc,G (x 1 , y 2 )v ǫc (y 2 , x 2 ) dν 2 (y 2 ) (6.56) = c(x 1 )G(x 1 , dz 1 )v ǫc (z 1 , z 2 )c(z 2 )G(z 2 , dz 3 )u ǫc,G (z 3 , y 2 )v ǫc (y 2 , x 2 ) dν 2 (y 2 ) dm(z 2 )

J i (ǫ) (6.49)
≤ c(x 1 )G(x 1 , dz 1 )v ǫc (z 1 , z 2 )c(z 2 )G(z 2 , dz 3 ) dm(z 2 ) sup z 3 ,x 2 u ǫc,G (z 3 , y 2 )v ǫc (y 2 , x 2 ) dν 2 (y 2 ) ≤ C ν 2 u 2 ,∞ c(x 1 )G(x 1 , dz 1 )v ǫc (z 1 , z 2 )c(z 2 )G(z 2 , dz 3 ) dm(z 2 ).

(6.57) Using (6.19), the fact that c is bounded and the fact that G(•, dz) is a probability density c(x 1 )G(x 1 , dz 1 )v ǫc (z 1 , z 2 )c(z 2 )G(z 2 , dz 3 ) dm(z 2 ) ≤ Cc(x 1 ). (6.58) Thus (6.55) is bounded independently of ǫ by C k j=1 c(x j ) dm(x j ) < ∞ (6.59) since c(x) is integrable. The other terms of K(ǫ) can be bounded similarly.

  j , y j+1 )   u(y i , x)u(x, y i+1 ) y j ) dν(x).

Theorem 6. 1

 1 Assume that sup z u(z, y)dm(y) < ∞, sup z u(y, z)dm(y) < ∞.(6.5)

sup zu

  ǫc,G (z, y)dν(y) ≤ C sup z u(z, y)dν(y), (6.38)

vK

  ǫc (y j , y j+1 )   V ǫc cGV ǫc (y i , y i+1 ) y j ), and K(ǫ) represents all the remaining terms. Noting that I(0) = I(0), we can write (6.47) asI(ǫ) -I(0) = I(ǫ) -I(0) + ǫ k i=1 J i (ǫ) + K(ǫ).(6.50)Note also that I(ǫ) of (6.48) is a special case of the I(ǫ) of (5.11) with ν(dx) = c(x) dm(x). Hence by (5i , x)u(x, y i+1 )c(x) dm(x) y j ).Since v ǫc (x, y) ↑ u(x, y) as ǫ ↓ 0, it follows by the Monotone Convergence Theoremi , z 1 )c(z 1 )G(z 1 , dz 2 )u(z 2 , y i+1 ) dm(z 1 )   k j=i+1 u(y j , y j+1 )   k j=1 dν j (y j ).To complete the proof of our Lemma it remains to show that lim ǫ→0

V

  ǫc cGV ǫc cGU ǫc,G (y j , y j+1 ) k j=1 dν j (y j ) (6.54)= v ǫc (y 1 , x 1 )cGV ǫc cGU ǫc,G (x 1 , y 2 ) • • • v ǫc (y k-1 , x k-1 )cGV ǫc cGU ǫc,G (x k-1 , y k ) v ǫc (y k , x k )cGV ǫc cGU ǫc,G (x k , y 1 ) k j=1 dν j (y j ) dm(x j ) = cGV ǫc cGU ǫc,G (x 1 , y 2 )v ǫc (y 2 , x 2 ) • • • cGV ǫc cGU ǫc,G (x k-1 , y k ) v ǫc (y k , x k ) cGV ǫc cGU ǫc,G (x k , y 1 )v ǫc (y 1 , x 1 ) k j=1 dν j (y j ) dm(x j ),where the last step is just a rearrangement. We can rewrite this ascGV ǫc cGU ǫc,G (x 1 , y 2 )v ǫc (y 2 , x 2 ) dν 2 (y 2 ) • • • (6.55) cGV ǫc cGU ǫc,G (x k-1 , y k ) v ǫc (y k , x k ) dν k (y k )cGV ǫc cGU ǫc,G (x k , y 1 )v ǫc (y 1 , x 1 ) dν 1 (y 1 ) k j=1 dm(x j ).