High-temperature decomposition of amorphous and crystalline cellulose: reactive molecular simulations

We study the thermal decomposition of cellulose using molecular simulations based on the ReaxFF reactive force field. Our analysis focuses on the mechanism and kinetics of chain scission, and their sensitivity on the condensed phase environment. For this purpose, we simulate the thermal decomposition of amorphous and partially crystalline cellulose at various heating rates. We find that thermal degradation begins with depolymerization via glycosidic bond cleavage, and that the order of events corresponds to a randomly initiated chain reaction. Depolymerization is followed by ring fragmentation reactions that lead to the formation of a number of light oxygenates. Water is formed mainly in intermolecular dehydration reactions at a later stage. The reaction rate of glycosidic bond cleavage follows a sigmoidal reaction model, with an apparent activation energy of 166 ± 4 kJ/mol. Neither the condensed phase environment nor the heating programme have appreciable effects on the reactions. We make several observations that are compatible with mechanisms proposed for cellulose fast pyrolysis. However, due to the absence of anhydrosugar forming reactions, the simulations offer limited insight for conditions of industrial interest. It remains unclear whether this is a natural consequence of the reaction conditions, or a shortcoming of the force field or its parameter set.


Introduction
Cellulose pyrolysis underlies the burning of wood and plant-based textiles, and is an essential step in the thermochemical conversion of plant biomass into solid, liquid and gaseous fuels. The thermal degradation process is remarkably complex, and is known to involve a large number of elementary reactions and an interplay of solid, liquid and gas-phase mechanisms. (Pecha et al. 2019) Describing the elementary reactions remains a fundamental challenge for a thorough understanding of cellulose pyrolysis. (Mettler et al. 2012b) During the past decade, computer simulation methods have found increasing use in this pursuit. (Ciesielski et al. 2020).
Recently, a number of studies have explored reactive force field methods (Harrison et al. 2018) as a means to investigate the fundamentals of cellulose pyrolysis. Molecular dynamics (MD) simulations based on reactive force fields are of particular interest, as they can predict both reaction mechanisms and rates without pre-determined chemistry. Moreover, the simulations capture the effects of a changing chemical environment. The computational cost of force field methods is substantially lower than that of quantum chemistry methods. This means that the structural hierarchy of native cellulose (Terrett et al. 2019) can be considered up to the level of microfibril bundles.
On the other hand, reactive molecular simulations are subject to the common limitations of atomistic simulation methods. Most importantly, time scale limitations dictate that pyrolysis reactions have to be studied at temperatures and heating rates that are significantly higher than experimental ones.
At the moment, there are two reactive force field formulations applicable to carbohydrates. These are the ReaxFF reactive force field (van Duin et al. 2001) and the adaptive intermolecular reactive bond-order potential (AIREBO) (Stuart et al. 2000). AIREBO-MD simulations have been used to study the mechanisms of irradiation damage in isolated cellulose chains (Polvi et al. 2012;Polvi and Nordlund 2014). Pyrolysis chemistry has been studied solely using ReaxFF-MD simulations. One factor that differentiates the earlier pyrolysis-related studies is their description of the condensed phase environment. This includes periodic crystals (Qiao et al. 2020), amorphous melts (Zheng et al. 2016;Chen et al. 2017;Rismiller et al. 2018) and individual chains . Some studies have used cellobiose as a model compound (Hao et al. 2020;Zhang et al. 2020a, b). The organization of cellulose into microfibrils, which is the native state in plant cell walls, has not been considered. Another question that has received limited attention is that of reaction kinetics.
Recently, we carried out stochastic ReaxFF-MD simulations on the thermal decomposition of isolated cellulose chains. (Paajanen and Vaari 2017) We found that the decomposition begins with depolymerization via glycosidic bond cleavage, and that the apparent activation energy is comparable to values reported in global mass loss kinetics studies. The findings were encouraging for further use of the modelling approach. These earlier simulations were limited in two important aspects. Firstly, we studied isolated chains, which are representative of gas phase conditions. This eliminates possible effects due to a condensed phase chemical environment (Seshadri and Westmoreland 2012;Hosoya and Sakaki 2013). Secondly, we only followed the simulations up to the first reactive event, which means that we could not observe reaction sequences. For these reasons, the simulations could only provide a partial description of the degradation process.
In the present work, we use ReaxFF-MD simulations to study the thermal decomposition of amorphous and partially crystalline cellulose (i.e. microfibrils). We focus on the mechanism and kinetics of chain scission, and their sensitivity on the condensed phase environment. We follow the complete transformation of cellulose into low molecular weight products (LMWPs). The current simulations reproduce many of our previous observations on isolated chains. Most importantly, we find that the decomposition begins with glycosidic bond cleavage, and that the sequence of events corresponds to randomly initiated chain depolymerization. We find that crystallinity has no appreciable effects on the mechanism or kinetics of chain scission, the evolution of the molecular weight distribution, or the LMWPs. Curiously, we do not observe levoglucosan (LGA), or other anhydrosugars, among the reaction products. We discuss this and other open questions regarding the validity of the results, and their relevance for conditions of practical interest.

Methods
MD simulations of thermal decomposition were carried out on two types of systems: amorphous and partially crystalline. Five structural variants and seven heating rates were considered for both system types. The simulated trajectories were analyzed for (1) the number of molecules and molecular species; (2) the distribution of molecular weight; (3) the state of glycosidic, pyranose ring and alcohol group bonds; (4) the chemical identity and origin of LMWPs; and (5) the presence of cyclic compounds. Lastly, the kinetics of glycosidic bond cleavage was studied by fitting reaction models and kinetic parameters against the temperature evolution of a conversion fraction. The molecular models, the simulation set-up and the postprocessing steps are detailed in what follows.

Molecular models
The studied systems consisted solely of cellulose. The starting point was a molecular model of an individual chain in the twofold helical screw conformation found in native cellulose. The atomic coordinates were generated using crystallographic unit cell data for the cellulose I b allomorph (Nishiyama et al. 2002). Models of individual chains with degrees of polymerization (DP) 10 and 20 were used to construct the amorphous and partially crystalline models. Native cellulose exists in an ordered state, in which the chains form fibrous aggregates known as microfibrils. (Jarvis 2018) The partially crystalline models are thus our primary representation of native cellulose. The amorphous models serve as a baseline for studying the effects that crystallinity might have on the pyrolysis reactions.

Amorphous systems
The amorphous systems were constructed as follows (see Fig. 1a for reference). Firstly, 16 chains of DP 20 were distributed in a periodic cubic simulation domain of 30 9 30 9 30 nm 3 . The molecules were distributed randomly both with respect to their location and orientation. The initial density was low, roughly 3 mg/cm 3 . The atomic coordinates were brought into a local minimum of potential energy using the conjugate gradient method, after which the system was compressed in a dynamic simulation of 1 ns in the isothermal-isobaric ensemble, at 300 K temperature and 1 atm pressure. During the simulation, the density would converge to a value of roughly 1.2 g/cm 3 (cf. typical cell wall density of wood, 1.50-1.53 g/cm 3 , which includes the contributions of hemicelluloses, lignins and bound water). The OPLS-AA force field for carbohydrates (Damm et al. 1997), adapted for cellulose by Paavilainen et al. (Paavilainen et al. 2011), was used in the preparatory phase to speed up the simulations.
After this, the interaction model was changed from OPLS-AA to the ReaxFF reactive force field. The ReaxFF parameter set was that of Chenoweth et al. (2008), originally developed to study hydrocarbon oxidation. This choice follows our previous work on the thermal decomposition of isolated chains (Paajanen and Vaari 2017). Moreover, a set of preliminary simulations showed that the Mattsson et al. (2010) and Rahaman et al. (2011) parameter sets, which have been previously applied to cellulose (Dri et al. 2015;Zheng et al. 2016;Chen et al. 2017) yield unrealistically low values for the apparent activation energy of depolymerization. Again, the system was brought into a local minimum of potential energy, and then equilibrated in a dynamic simulation of 100 ps in the isothermalisobaric ensemble, at 300 K temperature and 1 atm pressure. The ReaxFF bond information output was used to confirm that no chemical reactions had occurred during the preparatory phase. Five variants of the amorphous system were created, each starting from a different spatial distribution of the cellulose chains.

Partially crystalline systems
The partially crystalline systems were created as follows (see Fig. 1b for reference). Firstly, 24 cellulose chains of DP 10 were arranged in the 3-4-5-5-4-3 stacking (Yang and Kubicki 2020) to form a microfibril segment with a hexagonal cross-section. The fibril segment was placed in a periodic simulation domain of 12 9 15 9 5 nm 3 . The chains were set to be periodic along their length to mimic an infinitely long fibril. Additionally, four chains of DP 20 were arranged in a loose grid around the fibril, perpendicular to its longitudinal axis and parallel to the edges of the simulation domain. After this, the preparatory phase was similar to that of the amorphous systems. During the compression simulation, the DP 20 chains would adhere to the surface of the fibril segment and form a disordered layer between the fibril and its periodic images. Again, the compression resulted in a system with a density of roughly 1.2 g/ cm 3 . Lastly, it was confirmed that no chemical Fig. 1 Condensed-phase models of cellulose: a an amorphous system that consists of entangled chains; and b a partially crystalline system that consists of a periodic microfibril segment and surrounding disordered chains. The aligned chains of the microfibril segment are colored green and the disordered chains grey. The molecules are visualized using the stick-representation of covalent bonds. The red rectangle shows the periodic simulation domain. Simulations of thermal decomposition were carried out on five variants of both system types reactions had occurred during the preparatory phase. Five variants of the partially crystalline system were created, each using a different initial velocity vector, and resulting in a different distribution of loose chains around the fibril.
The partially crystalline systems approximate the internal structure of a cellulose microfibril bundle in spruce wood (Terrett et al. 2019)-with the exception that the hemicelluloses galactoglucomannan and glucuronoarabinoxylan are substituted by the loose cellulose chains. The mass ratio of disordered to ordered chains, 1:4, corresponds to that found for hemicelluloses and cellulose in conifer secondary cell walls (Scheller and Ulvskov 2010). The size and crosssectional shape of the microfibril is also based on studies on spruce wood (Fernandes et al. 2011;Paajanen et al. 2019). The use of a periodic model is motivated by the large aspect ratios of the microfibrils (Jarvis 2018).

Decomposition simulations
The amorphous and partially crystalline systems were subjected to thermal decomposition in simulations with a linear heating scheme. The simulations were carried out at constant heating rates of 0.6, 1, 2, 3, 4, 5 and 6 K/ps, and at 1 atm external pressure. The initial and final temperature were 300 and 3300 K, respectively, and the simulated time was adjusted according to the heating rate. A decomposition simulation was carried out for each combination of condensed-phase environment, structural variant and heating rate, which resulted in altogether 70 simulations.

Details of the simulation set-up
Both the preparatory and the decomposition simulations were carried out using the MD software LAMMPS (Plimpton 1995). As mentioned earlier, the OPLS-AA force field (Jorgensen et al. 1984) and the ReaxFF reactive force field (van Duin et al. 2001) were used to describe interatomic bonding in the preparatory phase, while only the latter was used in the decomposition phase. In ReaxFF, the electronegativity equalization method (Mortier et al. 1986) was used for charge equilibration with a precision parameter of 10 À6 .
The initial velocities of the atoms were assigned randomly to achieve the desired kinetic temperature. The Nosé-Hoover thermostat and barostat were used for temperature and (isotropic) pressure control, respectively. (Tobias et al. 1993;Shinoda et al. 2004) Time constants of 10 fs and 100 fs were used for temperature and pressure control, respectively. The total linear and angular momentum of the system was set to zero at 1 ps intervals. The equations of motion were integrated using the velocity-Verlet algorithm with a constant time step of 0.1 fs. The coordinates of the atoms and information on chemical bonding was output at 1 ps intervals.

Post-processing
The bond information files were used to determine the chemical composition of the systems as a function of time. This was done using in-house post-processing tools that were based on the OpenBabel Toolkit (O'Boyle et al. 2011) and the Ovito Open Visualization Tool (Stukowski 2010). The software were used to identify product species based on their elemental composition and bond topology, and to report their frequency of incidence. The canonical form of the simplified molecular-input line-entry system (SMILES) was used to produce human-readable reports. Products of intra and intermolecular reactions were identified based on the locations of their constituent atoms within the cellulose chains. The dissociation of glycosidic, pyranose ring and alcohol group bonds was identified from changes in the bond topology. Lastly, the smallest set of smallest rings (SSSR) algorithm (Figueras 1996) was used to identify cyclic compounds, and to count the number of pyranose, furanose and carbocyclic rings. Certain quantities of interest, such as the number of molecules and the system volume, were obtained directly from the output of the simulation software.

Reaction kinetics
The kinetics of glycosidic bond cleavage was quantified using the Shuffled Complex Evolution (SCE) algorithm (Duan et al. 1993) as implemented in the Gpyro software (version 0.8186) (Lautenberger and Fernandez-Pello 2009). Reaction models and the Arrhenius kinetic parameters were fitted against the temperature evolution of the fraction of intact glycosidic bonds. In other words, depolymerization was considered to be analogous to a thermogravimetric analysis (TGA) experiment, with the fraction of intact bonds corresponding to the residual mass. To comply with the usual operating range of the Gpyro software, the heating rate of the simulations was scaled from units of K/ps to K/min, as one would expect in practical experiments. After the kinetic parameters were estimated using the objective data in K/min, the parameters were scaled to the values corresponding to heating rates in K/ps. The Flynn-Wall-Ozawa (FWO) isoconversional method (Khawam and Flanagan 2006a) was used to support the estimation of kinetic parameters. Similarly, the method of generalized master curves (Sánchez-Jiménez et al. 2013) was used to support the screening of reaction models.

Overview
We performed 70 simulations of thermal degradation, which includes five parallel simulations for each combination of condensed phase environment and heating rate. On a general level, the degradation sequence is similar regardless of the crystallinity and the heating rate (see Fig. 2). Thermal motion induces reactions that break cellulose down into progressively smaller molecules and molecular fragments. This leads to an increase in internal pressure, which drives the expansion of the simulation domain and ultimately results in vaporization. The degradation reactions continue in the gas phase, and ultimately lead to a system that consists solely of LMWPs. The parallel simulations display repeatability for all of the recorded quantities. In the following sections, we discuss in detail the condensed phase reactions that break cellulose down to the oligo and monosaccharide level. We also briefly discuss the LMWPs that are formed in the condensed phase.
It should be noted, that the vaporization step is possible due to the simulation protocol that allows expansion against atmospheric pressure. This is one of the ways to deal with pressure build-up, which would otherwise lead to extreme and unwanted conditions. An alternative approach is to perform the simulations in conditions of constant volume, and to regularly remove species that are deemed volatile. In this case, the user controls the residence time of the volatiles, and needs to provide the criteria for identifying them, as well as their vaporization rates. These parameters necessarily affect the frequency of secondary reactions, and, by extension, the overall sequence of reactive events. While the chosen approach has a similar limitation with respect to the choice of external pressure, we consider it less ambiguous and preferable for the purpose of studying primary reactions.

Depolymerization reactions
The evolution of the molecular weight distribution shows that chain segments of all intermediate lengths are present at the onset of vaporization (Fig. 2c). The following analysis shows that chain scission events are almost exclusively depolymerization reactions that occur via glycosidic bond cleavage. This is demonstrated using bond dissociation statistics for the amorphous system at the lowest heating rate ( Fig. 3a-d). The other heating rates and the partially crystalline systems lead to almost identical data sets. Figure 3a shows the fraction of cleaved glycosidic bonds as a function of temperature, and the distribution of the dissociation events between the C1'-O4 and C4-O4 bonds. Two key observations can be made. Firstly, glycosidic bond cleavage does take place, and all of the bonds are cleaved over the course of the simulations. Secondly, the dissociation takes place at the C1'-O4 bond in roughly 60% of the cases, which means that the glycosidic oxygen tends to remain bonded towards the reducing end of the chain. It is also worth noting that the transition into gas phase occurs when roughly 40% of the bonds have been cleaved. Glycosidic bond cleavage leads to chain scission, and thus to the formation of progressively shorter chains. However, chain scission can also occur via ring fragmentation reactions. Figure 3b shows the fraction of anhydroglucose units (AGUs) in which one or more of the ring bonds have been cleaved. Chain scission via ring fragmentation requires the cleavage of at least two of the pyranose ring bonds (n [ 1). At the onset of vaporization, this has occurred in less than 15% of the AGUs, which implies that glycosidic bond cleavage is the primary chain scission pathway. As further evidence, Fig. 3c shows the timing of ring fragmentation reactions relative to glycosidic bond cleavage.
The timing is considered from the viewpoint of an AGU and its two bonds with the neighboring units (see inset of Fig. 3c). In this setting, glycosidic bond cleavage precedes ring fragmentation for more than 90% of the AGUs. This confirms that glycosidic bond cleavage dominates chain scission and, consequently, the molecular weight distribution.
The timing of ring opening and fragmentation reactions relative to glycosidic bond cleavage can be used to deduce further details. Firstly, glycosidic bond cleavage precedes ring opening for roughly 85% of the AGUs (Fig. 3c). This means that ring opening reactions typically occur at chain ends. Secondly, both of the glycosidic bonds are cleaved before ring opening for roughly 25% of the AGUs, and before ring fragmentation for roughly 70% of the AGUs (not shown). This means that individual repeat units are typically released in their acyclic form. Lastly, the timing can also be considered from the viewpoint of a glycosidic bond and its two neighboring AGUs. In this setting, ring opening precedes glycosidic bond cleavage for roughly 85% of the glycosidic bonds (not shown). This is the same number as for the reverse situation discussed above. It is only possible if the depolymerization occurs via a chain-end mechanism in which glycosidic bond cleavage and ring opening reactions alternate along the chain.
To be thorough, we also looked at the cleavage of the alcohol group bonds. Figure 3d shows the fraction of AGUs in which one or more of the alcohol group bonds have been cleaved (i.e. the C2-O2, C3-O3, C5-C6 or C6-O6 bonds). At the onset of vaporization this has occurred in roughly 10% of the AGUs, which implies that reactions involving the alcohol groups are not central to the depolymerization stage. In summary, we reach the following conclusions: (1) depolymerization occurs predominantly via glycosidic bond cleavage, (2) it frequently leads to ring opening reactions of the neighboring AGUs, (3) there is a propagating chain-end mechanism that involves alternating glycosidic bond cleavage and ring opening, (4) ring fragmentation occurs mostly in the monomeric species released in the process, and (5) reactions of the alcohol groups only become important at a later stage.
As mentioned earlier, we found that the bond dissociation statistics of the amorphous and crystalline systems are similar. We made the same observation with regard to the depolymerization sequence, the depolymerization kinetics and the distribution of LMWPs. This indicates that the degradation process is insensitive to crystallinity. Several experimental studies have reached the opposite conclusion regarding cellulose slow and fast pyrolysis. Sample crystallinity has been reported to affect both the onset temperature of volatilization (Kim et al. 2010;Wang et al. 2013) and the yields of various volatile products. (Wang et al. 2013(Wang et al. , 2014Liu et al. 2013;Mukarakate Fig. 3 Distribution and timing of bond dissociation events: temperature evolution of a the fraction of cleaved glycosidic bonds; b the fraction of AGUs with one or more cleaved pyranose ring bonds; c the fraction of AGUs for which glycosidic bond cleavage (t GBC ) occurs before or after ring opening (t RO ) and ring fragmentation (t RF ) reactions; and d the fraction of AGUs with one or more cleaved alcohol group bonds. The bonds are highlighted in the insets in green and yellow. The data is averaged over the simulations on the amorphous system at 0.6 K/ps heating rate. The vertical dashed lines indicate the onset of vaporization et al. 2016; Leng et al. 2018) However, some of these reports have been challenged on the basis of possible heat and mass transfer limitations, as well as sample purity issues. (Zhang et al. 2014). The validity of the simulated behavior thus remains unclear. To simplify the remaining discussion, we only present the results for the amorphous systems (with the exception of the kinetic analysis, where additional plots are not required).

Chain depolymerization
We then studied how the depolymerization proceeds within the individual chains. For this purpose, we determined the order of glycosidic bond cleavage within each chain. We first looked at a number of example cases, and observed that series of three or more bonds were frequently cleaved in consecutive order (see Fig. 4). The cleavage would proceed either towards the reducing or the non-reducing end, with seemingly no preference. We then determined the relative frequency of intra and end-chain events during the first two depolymerization steps. We found that each glycosidic bond is an equally probable site for the first cleavage, but the second one occurs more probably at the newly formed chain ends. Specifically, the frequency of chain-end cleavage is roughly ten times higher in the second step. These observations are indicative of a chain depolymerization mechanism.
We then looked at how the molecular weight distribution evolves as a function of the fraction of cleaved glycosidic bonds (i.e. conversion) (Fig. 5a), and compared the distributions with those predicted by a stochastic chain depolymerization model (Fig. 5b) (Simha and Wall 1952). In the stochastic model, the depolymerization sequence is controlled by four parameters: the probabilities of (1) intra-chain initiation (p i ), (2) end-chain initiation (p e ), (3) propagation (p p ) and (4) termination (p t ). Contrary to the model of Simha and Wall, we assumed the probability of termination to be constant, and thus independent of the chemical environment. Based on the earlier observations, we assumed that the process is randomly initiated (p i % p e ) and that the probability of propagation is significantly higher than that of initation (p p [ p i;e ). With these assumptions, the molecular weight distribution is controlled solely by the probabilities of propagation and termination (p p , p t ).
Through manual fitting, we found a good correspondence with the parameters p i =p e ¼ 1; p p =p i ¼ 13, and p t ¼ 0:43. Figure 5a-d compare the molecular weight distribution observed in the simulations and that produced by the chain depolymerization model. Figure 5a, b show the distributions at a glance, and Fig. 5c, d provide a visually easier comparison of the low-DP range. The ReaxFF-MD results are an average over the simulations on the five amorphous systems at the seven different heating rates. To make the comparison meaningful, we removed the effects of fragmentation reactions from the results. The broad agreement supports the interpretation that, in the ReaxFF-MD simulations, cellulose thermal decomposition is driven by randomly initiated chain depolymerization. It should be noted, that the frequent ring opening reactions lead to the absence of cyclic monomeric species among the reaction products.
A random initiation step is in agreement with our previous work on isolated chains, where we observed a uniform scission probability across the glycosidic bonds. (Paajanen and Vaari 2017) It is also compatible with the recent findings of Dauenhauer and co-workers (Krumm et al. 2016;Zhu et al. 2017). Their thin film pyrolysis experiments show that cellulose undergoes a transition from end-chain to intra-chain depolymerization at 467°C (referred to as the reactive melting point), with the intra-chain mechanism dominant at higher temperatures. On the other hand, several studies suggest that cellulose fast pyrolysis involves a chain depolymerization mechanism, in which LGA-the major volatile product-is formed during the depropagation step. (Mayes and Broadbelt 2012;Vinu and Broadbelt 2012;Hosoya and Sakaki 2013). While we do observe chain depolymerization, we do not find LGA or other anhydrosugars among the products. We assign this to the high temperatures used for accelerating the rate of the primary reaction of interest in this paper, glycosidic bond cleavage. However, at the same time, all other homolytic bond cleavage reactions are accelerated, most of which lead to ring opening. In addition, these bond cleavages produce radicals that may interfere with the rest of the decomposition sequence. High temperatures also quickly degrade crystallinity in the partly crystalline systems, and contribute to the expansion of the simulation box. These effects tend to suppress inter-chain interactions thought to be important for LGA formation (Hosoya and Sakaki 2013).

Depolymerization kinetics
We then studied the kinetics of depolymerization using the fraction of cleaved glycosidic bonds as the conversion variable. Figure 6a shows the temperature evolution of conversion for the amorphous and partially crystalline systems. The conversion curves at each heating rate are closely similar, which implies that the condensed phase environment does not significantly affect the reactions. We then determined the apparent activation energy of depolymerization using the FWO isoconversional method (Khawam and Flanagan 2006a). The FWO method is based on the approximately linear dependence between the logarithm of the heating rate and the inverse temperature at a given conversion: ln b a % ln A a E a;a g a ð ÞR À 2:315 À 0:457 where b a is the heating rate, A a is the pre-exponential factor, E a;a is the apparent activation energy, g a ð Þ is the (integral) reaction model, T a is the temperature, and R is the universal gas constant. The subscript a refers to a constant value of conversion. Importantly, the activation energy can be determined from a linear fit of ln b a versus T À1 a without assuming a reaction model. Figure 6b shows an example fit at a conversion of 0.4. A linear fit was appropriate at all conversion levels. Figure 6c shows the estimates for the activation energy at several conversions between 0.1 and 0.9. We find that the activation energy is roughly constant up to the conversion of 0.8, which is well after the transition into gas phase. The average activation energy is 166 ± 4 kJ/mol at conversions below 0.8, and then increases to 191 ± 4 kJ/mol at 0.9 conversion. The values are remarkably close to what we observed previously for coiled and elongated chains in vacuum (168 ± 2 and 188 ± 2 kJ/mol, respectively; DP-16) (Paajanen and Vaari 2017). At the time, we attributed AGUs are shown in blue, yellow and red, respectively. Fragmented-ring AGUs that are not part of a chain segment are hidden the lower activation energy in coiled chains to the presence of self-interactions. The current results support this conclusion, as we observe significant clustering of chain segments after the transition into gas phase. Both experimental and theoretical studies have suggested an activation energy for cellulose thermal depolymerization in the 190-200 kJ/mol range. (Mamleev et al. 2007a, b;Sánchez-Jiménez et al. 2011;Agarwal et al. 2012) The predicted average activation energy is roughly 13% below this range. We used the value 166 ± 4 kJ/mol as a fixed parameter in the subsequent SCE search for the pre-exponential factor and the reaction model.
During the SCE screening, we found that common reaction models could not reproduce the simulated conversion curves. The closest correspondence was achieved using an Avrami-Erofeyev (i.e. nucleation and growth) type model (Khawam and Flanagan 2006b), which produced a satisfactory match for conversions below 0.5, but failed at higher values. Our  (c) and (d) for monomers and oligomers, respectively. In (c), the dashed line shows the frequency of monomers due to random depolymerization. In (d), the darker coloring indicates the standard error of the mean, and the lighter one the standard deviation, as predicted by the stochastic model. The ReaxFF-MD results are averaged over all simulations on the amorphous systems earlier observations support a consistent mechanism throughout the simulations, and we thus continued the screening using the method of generalized master curves (Sánchez-Jiménez et al. 2013). The master curve method is based on expressing the reaction rate in terms of a reduced time variable, s: where t is time. The reduced reaction rate becomes: where f a ð Þ is the (differential) reaction model. Consequently, In other words, the generalized reaction rate can be used to reproduce the reaction model up to a scaling factor. Fig. 6 a Fraction of intact glycosidic bonds (1 À a) as a function of temperature in the partially crystalline and amorphous systems. b Example of a linear fit of ln b a versus T À1 a for the FWO method. c Apparent activation energy of depolymerization determined using the FWO method. d Generalized master curves calculated from the simulations and a number of analytical reaction models. The error bars in (c) show the standard error of the linear least squares fit. The grey color in (c) shows a range of activation energies proposed for cellulose thermal depolymerization (see text for references). In (d), the ReaxFF-MD data points are averaged over all simulations Figure 6d shows the generalized master curves determined from the ReaxFF-MD simulations, the stochastic chain depolymerization model, and a number of commonly used analytical reaction models. We made the following assumptions in deriving the reaction model from the chain depolymerization scheme: (1) the initiation and propagation steps have the same activation energy, (2) their rates are directly proportional to the number of susceptible bonds, and (3) their relative rate is directly proportional to their relative probability. Thus, the reaction model can be expressed as a function of the number of intra-chain bonds, end-chain bonds, and end-chain bonds due to initiation: where N i and p i are the number and scission probability of bonds of type i, respectively, and N ref is the total number of bonds at zero conversion. Due to the random initiation step, the master curve produced by the chain depolymerization model is insensitive to the initial molecular weight distribution.
Two key observations can be made from Fig. 6d. Firstly, the reaction model derived from the chain depolymerization scheme reproduces the master curve of the ReaxFF-MD simulations satisfactorily. This implies that the reaction model formulated in Eq. (5) is a reasonable approximation of the depolymerization kinetics. Secondly, the figure shows why the common analytical models fail to reproduce the simulated conversion curves. Acceleratory, deceleratory and linear models inevitably fail, as the reaction clearly has an acceleratory (a\0:3) and a deceleratory phase (a [ 0:3), and thus a peak reaction rate at an intermediate conversion (a % 0:3). This group of models includes reaction-order (F), geometrical contraction (R), diffusion (D) and power-law nucleation models (P), among others. On the other hand, an Avrami-Erofeyev type model (A) can produce a satisfactory match in the acceleratory phase, but at the cost of deviations in the deceleratory phase. Figure 6d shows examples of the analytical models with given parameters (e.g. F1 refers to a reaction-order model with an order parameter of one).

Low molecular weight products
Lastly, we looked at the formation of LMWPs before the transition into gas phase. We determined the chemical identity and frequency of species with mass below a given threshold, and classified them into products of intra and intermolecular reactions. Among the products of intramolecular reactions, we further identified products of fragmentation reactions. The classification was based on tracing the atoms to their original locations within the molecular topology of cellulose. With the exception of dehydration reactions, we did not consider hydrogen transfer to constitute an intermolecular reaction. We chose m 2 3 m AGU for the mass threshold to exclude partially dehydrated monomeric species. Figure 7 shows the number and mass fraction of LMWPs as a function of conversion. At the onset of vaporization, the mass fraction of LMWPs has reached roughly 10%, and the number fraction roughly 40%. With the exception of water, the LMWPs are almost exclusively products of intramolecular reactions. After the transition into gas phase, encounters between the molecules become infrequent and intermolecular reactions cease. The LMWPs formed in the condensed phase are almost exclusively products of fragmentation reactions. We identified more than one hundred different species. However, a relatively small group of frequent LMWPs constitute the majority both by number and by mass. Frequently observed species Fig. 7 Number and mass fraction of low molecular weight products as a function of conversion. The data represents an average over the simulations on the amorphous systems include, in the order of increasing molecular weight, water, formaldehyde, carbon dioxide, ethanedial, hydroxyacetaldehyde, 1,2-ethenediol and hydroxypropanedial. Together, they account for roughly 75% of the LMWPs by number, and 60% by weight. With the exception of water, the frequent species are formed from specific fragments of the AGUs, as summarized in Fig. 8 and Table 1, and briefly discussed below.
Formaldehyde is formed exclusively from the hydroxymethyl group, i.e., the C6-O6 fragment. Its formation thus proceeds via the cleavage of the C5-C6 bond. Carbon dioxide is formed almost exclusively ([ 95%) from the O4-C1-O5 fragment, which contains both the glycosidic and the pyranose ring oxygens. This requires both glycosidic bond cleavage and ring fragmentation. The formation of carbon dioxide has been linked to that of formic acid, both derived from the C1 carbon. (Mettler et al. 2012a) Ethanedial (i.e. glyoxal) is formed exclusively from the C1-C2, C2-C3, C3-C4 and C4-C5 segments. Its formation thus requires ring fragmentation, either alone (C2-C3 segment) or accompanied by glycosidic bond cleavage (C1-C2, C3-C4 and C4-C5 segments). The C1-C2 segment accounts for the majority of the cases ([ 75%), which indicates that the latter route is dominant. Hydroxyacetaldehyde (i.e. glycolaldehyde) is formed exclusively from the C5-C6 segment, and its formation thus proceeds via ring fragmentation. Banyasz et al. (2001a, b) and Piskorz et al. (1986) have proposed that hydroxyacetaldehyde is a product of ring fragmentation, and that the reaction path that leads to its formation competes with that of LGA. Moreover, the yield of HAA has been observed to increase with temperature at the expense of LGA in cellulose fast pyrolysis. (Luo et al. 2004) These observations are compatible with ours, and provide an explanation for the lack of LGA and other anhydrosugars. 1,2-ethenediol is formed from all possible carbon-carbon segments, including the C1-C2, C2-C3, C3-C4, C4-C5 and C5-C6 segments. The C2-C3 and C5-C6 segments account for roughly half of the cases. These routes can proceed via ring fragmentation alone, while the others also require glycosidic bond cleavage. Hydroxypropanedial (i.e. glucic acid) is formed exclusively from the C1-C2-C3, C2-C3-C4 and C3-C4-C5 segments. Its formation thus requires both glycosidic bond cleavage and ring fragmentation. Lastly, water is formed mainly in intermolecular dehydration reactions involving the C2, C3 and C6 alcohol groups. Water formation in the condensed phase is minor, and amounts to less than 0.3% of the dry mass. The yield increases significantly after the transition into gas phase, and reaches values above 5% of the dry mass. It should be noted, that our models only account for chemical loss of water, as absorbed water is not present. This means that possible catalytic effects due to the absorbed water are not captured.
The experimental reference yields given in Table 1 are from the cellulose powder and thin film pyrolysis experiments of Mettler et al. (2012a). We chose this work as a reference since the observed product distribution had the most overlap with that of our simulations. Due to the different reaction conditions, the experimental yields should be considered a point of reference, and not an attempt at validation. The Authors also report ab initio MD simulations on the thermal decomposition of a-cyclodextrin, which they use as a small-molecule surrogate for cellulose. Based on the simulations, they propose reaction pathways for the formation of formaldehyde, carbon dioxide, ethanedial and hydroxyacetaldehyde. The predicted segments of origin are in agreement with those given in Table 1, with the exception that ethanedial formation only involves the C3-C4 segment. The Authors Table 1 Frequent low molecular weight products and their origins within the anhydroglucose unit (see Fig. 8   do not report 1,2-ethenediol or hydroxypropanedial among the reaction products.

Conclusions
We carried out ReaxFF-MD simulations on the thermal decomposition of amorphous and crystalline cellulose. Our aim was to obtain predictions on the elementary mechanisms of cellulose pyrolysis, and to evaluate the capability of ReaxFF-MD for this purpose. The simulations reproduce many of our previous observations on the thermal decomposition of isolated chains. Most importantly, we find that the process begins with glycosidic bond cleavage, and that the apparent activation energy is comparable to values reported in global mass loss kinetics studies. In the condensed phase simulations, we follow the reactions further, and find that the depolymerization proceeds as a randomly initiated chain reaction. We find that crystallinity has no appreciable effects on the mechanism or kinetics of chain scission, the evolution of the molecular weight distribution, or the identity of LMWPs. Frequently observed LMWPs are formed in fragmentation reactions, and can be traced to specific locations within the AGUs. They include formaldehyde, carbon dioxide, ethanedial, hydroxyacetaldehyde, 1,2-ethenediol and hydroxypropanedial. We do not observe the formation LGA or other anhydrosugars, which raises the question of whether ReaxFF-MD can adequately represent conditions of slow or fast pyrolysis. It is not clear whether the absence of anhydrosugars is a natural consequence of the reaction conditions, or a shortcoming of the ReaxFF formalism or the used parameter set. The notable formation of hydroxyacetaldehyde, which has been shown to compete with that of LGA, in fact suggests the former. The comparability between high temperature simulations and pyrolysis in natural and industrial processes remains a central question for future studies. Systematic force field evaluation should, in our opinion, begin from the pyrolysis chemistry of simple sugars, especially that of glucose. At the same time, the question of cellulose pyrolysis alone warrants the development of a ReaxFF parameter set optimized for carbohydrate chemistry. In this connection, it should be evaluated whether an explicit description of electron transfer is required. In fact, such extension to the ReaxFF formalism already exists (Islam et al. 2016). Lastly, to overcome the limitation of artificially high temperatures and heating rates, the possibilities of acceleration methods, such as the accelerated ReaxFF reactive dynamics (aARRDyn) (Cheng et al. 2014) and reactive parallel replica dynamics (RPRD) (Joshi et al. 2013), should be explored thoroughly.