Impact of oxidation morphology on reduced graphene oxides upon thermal annealing

Thermal reduction of graphene oxide (GO) is an essential technique to produce low-cost and higher quality graphene-based materials and composites used today in a plethora of applications. However, despite a demonstrated efficiency of high-temperature annealing in reducing the oxygen content of GO, the impact of the morphology of the initially oxidized samples on the restored sp2 graphene plane versus remaining sp3 imperfections remains unclear and out-of-control. Here using classical molecular dynamics, we simulate the process of thermal reduction on several GO samples for a variety of initial conditions and elucidate how both the concentration of oxygen functional groups and their spatial distribution jeopardize the reduction process efficiency. Our simulations suggest thermal annealing strategies to further optimize the crystallinity of reduced GO, enhancing their transport properties and hence making the resulting composites even more performant for electronic applications.

Common fabrication methods of graphene include mechanical exfoliation [1], chemical vapour deposition (CVD) [17,18], epitaxial growth [19], electro-chemical and solution exfoliation [20,21]. These methods, however, are known to yield only small-scale production. Alternatively, large-scale production (at a comparatively lower cost) can be obtained via the reduction of graphene oxide (GO) [22], which is the current form of graphene-based materials used in today composites for industrial applications. Reduced GO-based composites are intensively used as additives to plastics, incorporated in glass-reinforced polymers or in concrete to enhance strength and thermal conductivity performances. Besides, such graphene-related materials are also suitable for coatings and printing applications [23][24][25].
A GO sample typically results from the oxidation process of graphite (under strong acid/base treatments). It is made by oxygen covalently-functionalized, atomically thin carbon sheets, showing O:C ratios between 0.3 and 0.5. The presence of oxidizing species is responsible for a dramatic degradation of the properties with respect to pristine graphene. In particular, GO is electrically insulating (in opposition to graphene) since the presence of sp 3 carbon atoms disrupts the flow of charge carriers through sp 2 bond networks [26][27][28]. Similarly, a strong reduction in the thermal conductivity (of 90%) is also observed in GO for degrees of oxygen functionalization as small as 5% [29].
This poses the need for a reduction process transforming GO into the parent state of graphene, yielding the so-called reduced graphene oxide (RGO). For reducing GO, mostly post-synthesis chemical and thermal treatments are employed, favoring the loss of oxygen content, and thus the partial restoration of graphene-like properties. In particular, within chemical treatments, reducing agents are often used such as hydrazine and Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. hydroiodic acid. Their employment, however, generates products with low C:O ratio and a comparatively poor quality with respect to graphene produced via mechanical/physical exfoliation or CVD. On the contrary, thermal reduction is generally performed through high-temperature annealing (in vacuum, inert of reducing atmosphere) which usually results in much more effective reduction, efficiently restoring sp 2 carbon domains and improving the electrical and thermal properties of GO sheets [14,30].
In this respect, theoretical computations based on the density functional theory and Monte Carlo techniques have explored the basic reduction mechanisms and the thermodynamic stability of the reduced samples [31][32][33][34][35][36]. Today there is general consensus on the fact that: (i) GO samples contain oxygen mainly in the form of epoxy (i.e. the bridge site oxygen) as well as hydroxyl (i.e. -OH) groups; however, the hydroxyl groups are less stable, and so, after annealing, most of the hydroxyl groups are removed. (ii) The key factor underlying the GO stability is the tendency of oxygen to diffuse, to agglomerate and to form highly oxidized domains surrounded by regions of pristine graphene. (iii) Within these agglomerates, all the binary decomposition reactions producing O 2 or H 2 O are hindered by both geometrical and energetic factors; moreover, the oxygen binding energy in an oxidized island increases with the size of the functionalized region. (iv) The formation of lattice damage in GO is related to the surface density of epoxy species: at low coverage, these undergo reversible desorption, while they create the precursors for the release of CO/CO 2 mixtures upon increasing coverage.
An overall picture has emerged where RGO derives from the interplay among several chemical reaction processes whose effect is to alter the sample morphology by removing functionalizing species (and carbon atoms) from the basal plane. In this respect, the morphology of the GO configuration and the initial density of oxidizing species represent pivotal factors that dictate the predominant chemical mechanisms at play and hence the final recovery of the hexagonal lattice cristallinity. The occurrence of the mentioned reactions is triggered by thermal energy, which represents the ultimate force driving the reduction process. Such a factor, however, is unavoidably responsible also for structural modifications of the basal plane and for its deviations from the planar configuration. Thermal energy relaxation processes, indeed, are expected to determine rearrangements of the atomic positions on large spatial scales possibly leading to strong in-and out-of-plane displacements in RGO. As an extreme case, strong thermal stresses could produce extended disruptions in the original plane, jeopardizing the overall reduction process. The impact of these rearranging mechanisms on the overall cristallinity of RGO is still unclear. In particular, while the nature of the basic chemical reduction mechanisms has been largely investigated, a systematic analysis of the interplay between such chemical reactions and the thermal relaxation processes (as a function of the annealing temperature) still needs to be addressed. A major factor hindering such an achievement is represented by the multi-scale nature of the problem in both the spatial and time domain. Thermal reduction is, indeed, a dynamical process taking place in systems with macroscopic dimensions: on the one side, it involves chemical reactions occurring within a reduced number of atoms; on the other side, the equilibration of thermal energy causes rearrangements of the atoms on a much larger scale. Also the time scales involved are different, with chemical reactions occurring in few femtoseconds while purely diffusive phenomena could take from seconds to hours. A fully quantum-mechanical description of the dynamical evolution of large atomistic systems would clearly satisfy the previous requirements but it is practically unfeasible. In this respect, classical molecular dynamics with reactive force fields represents a valid computational tool to investigate such a problem [37,38]. Here we systematically analyze the thermally reduced GO on both the chemical and structural level. In particular, by simulating the thermal reduction of several GO samples, we investigate the impact of the annealing temperature on the chemical content and structural properties of RGO. By monitoring the number of sp 2 carbon atoms and the oxygen-to-carbon ratio as a function of the annealing temperature, we study the effect of the density of oxidizing agents adsorbed and of the spatial distribution of the functionalized areas onto the efficacy of the thermal treatment. Moreover, in order to estimate the mechanical integrity of the systems, we calculate the strain-field and the atomistic out-of-plane displacements. The knowledge of these indicators represents a fundamental step towards a correct interpretation of the experimental results, which generally relate such parameters to transport coefficients (e.g. electron and thermal conductivities). We show how large initial concentrations of oxygen-based functionalizing groups and large oxidized areas do not lead to a sizeable recovery of the plane crystallinity. At the same time, strong local rearrangements of the atomic configuration are obtained with noticeable deformations of the basal plane whose regularity is consequently lost.

Methods
Molecular simulations performed in this work were carried out within the LAMMPS package [39] to study the formation of RGO materials and investigating the effect of the thermal reduction process.
To properly describe the bond association/dissociation events during the reduction process, the reactive force-field ReaxFF [40] was adopted. ReaxFF is a general bond-order dependent potential providing energies, transition states, reaction pathways, and reactivity trends in a number of hydrocarbon-oxygen systems in good agreement with quantum mechanical calculations and experiments [41,42]. Furthermore, it has been previously fruitfully employed to study the structural evolution of RGO/GO materials in several experimental conditions [37,38,43].
In our simulations we used a 160 Å×160 Å graphene sheet with epoxy and hydroxyl functional groups attached to both sides of the sheet. The total number of Carbon atoms in the systems is as large as 10488, with periodic boundary conditions applied along the graphene basal plane. The initial concentrations of epoxy and hydroxyl groups adsorbed are in a ratio of 1:1; different values of such a parameter have proved not to modify the main conclusions of our simulations to a significant extent.
In order to disentangle the effects of the spatial aggregation of the functionalizing groups and their spatial density, two different sets of initial GO configurations have been built: (a) for a fixed C:O concentration, four systems have been generated with epoxy and hydroxyl groups agglomerated to form oxidized areas with increasing average size; (b) for a given morphological distribution of oxidized regions (and total defected area), we considered four configurations with an increasing density of epoxy and hydroxyl groups filling the defected areas. In all the samples the actual position of the oxidizing groups within the defected clusters is chosen in a random fashion.
Once the initial groups were distributed on the plane, the GO sheet was subjected to a thermal annealing process consisting of a multi-step procedure. After a relaxation at 5 K for 12.5 ps, the system temperature was gradually varied by iteratively repeating the following steps: (1) the temperature is linearly increased of ΔT=150 K in 25 ps; (2) the temperature is kept fixed in the following 200 ps. The above iterative scheme is repeated until a maximum temperature larger than 2000 K has been reached. The length of the different time intervals is chosen to be sufficiently long to host the reduction process: indeed, after 200 ps no significant changes in the morphology of the sheet have been observed in most of the simulations. To perform the structural analysis on the systems, the structures obtained at the end of each iterative step have been further relaxed at 300 K for 1.25 ps. The simulations were performed in a constant-temperature (NVT) ensemble with a Berendsen thermostat for temperature control and a time step of 0.25 fs.
While the present simulation setup and protocol are likely limited in their quantitative accuracy, all previous successful applications of the same force field here adopted make us confident that our results are reliable enough to identify a fundamental trend of the most relevant quantities of interest as a function of the input parameters. In particular, we believe that the conclusions of this work still hold when the applied heating rate is much larger than the one considered here, like in real annealing situations. In such cases the diffusive mechanisms are indeed expected to have a stronger impact on the evolution of the systems. In order to consider this aspect, the effect of impurity diffusion (previous and ongoing) is encoded in the initially different spatial distributions of oxidized areas.

Results
In this section we present the results of our simulations of the thermal reduction on some significant samples of GO. We discuss the morphological properties of the resulting RGO structures as a function of the temperature and analyze the relevance of the involved parameters. For the sake of clarity, we consider separately the effects of the different density of defects in the functionalized areas and of the initial GO morphology. This structure represents a convenient scheme to provide some meaningful physical insight.

Density of functionalizing groups
We begin investigating the impact of the density of oxygen groups in GO on the morphological features of RGO obtained via thermal annealing. The presence of functionalized areas with different concentrations of oxidizing groups has emerged as a characteristic feature of GO systems. More specifically, they can result from local aggressive chemical oxidation steps or from the post-oxidation diffusive motion of oxygens on the graphene surface [33]. It is clear that the number of possible configurations of defected areas and functionalizing groups experimentally obtainable is remarkably large: they can differ in shape, dimensions and the way oxygen groups spatially distribute on their surface. Performing simulations of thermal annealing on all of them is impossible and we need an alternative strategy which could be nevertheless result of some help at a reasonable computational cost. Therefore, we focus on an individual GO system characterized by a unique spatial distribution of functionalized areas covering 50% of the total area. We take it as a representative of a general GO sample undergoing thermal annealing. For the given spatial distribution of defected areas, we considered four samples of GO differing for the density of the distributed oxidizing agents. In detail: by keeping fixed the total area of defected regions in the sample and their spatial distribution, we oxidized them with a different amount of epoxy and hydroxyl groups. The configurations of GO sheets that we analized are shown in We, at first, monitored the amount of oxygen adsorbed onto the surface at different annealing temperatures. The corresponding evolution is presented in figure 2. A meaningful variation of the oxygen content is observed in all the cases for temperatures larger than 400 K. Below this value, thermal energy, which can be considered the ultimate cause triggering all the possible basic mechanisms in the reduction process (e.g. desorption of oxygen compounds from the surface, complex chemical reactions among different oxidizing agents forming other products, release of CO 2 molecules in the environment and surface diffusion of impurities) is not sufficient to determine the desorption of oxygens from the graphene plane. On the contrary, at higher temperatures a monotone decrease of oxygen is found. Importantly, we note a stronger reduction of the oxygen content for higher functionalizing densities: for a given annealing temperature, the higher the initial density of oxygencontaining species, the higher the amount of desorbed oxygen. In this respect, the evolution of the quantity of sp 2 -hybridized carbon atoms offers a different perspective on the mechanisms taking place during the reduction process. Such a quantity is usually taken as representative parameter of the purity of the RGO sample: it is reported in figure 3 as a function of temperature normalized to its initial value. For all the systems, a similar trend of sp 2 carbons can be observed. After an initial transient over which reduction is not at work (in agreement with the same range behaviour of figure 2), sp 2 carbons generally increase for any T>800 K. Such an increase is more evident in low-density functionalized samples, where up to 10% of the original amount of carbons is regained  for sufficiently high temperatures. Systems with a higher concentration of oxygen groups undergo a weaker reduction process and only a small percentage of carbons goes back to the unoxidized state.
These results suggest a twofold conclusion. First of all, we remark that for sufficiently high annealing temperatures, no improvement in the cristallinity of the basal plane is obtained. Although a decrease of the oxygen content is simultaneously gained, the amount of sp 2 atoms remains almost unchanged. Since it has been pointed out elsewhere that a strong connection between some RGO properties (e.g. electrical conductivity) and the sp 2 carbon atom can be observed [27], our analysis suggests the existence of saturating trend of such observables for increasing annealing temperatures. Next, we observe that for a given annealing temperature, the amount of sp 2 -hybridized carbons of the reduced system decreases for larger values of the initial concentration of functionalizing groups. This fact is evident in the two low-density samples for T>1500 K and can be better explained by looking at the typology of chemical desorption processes occurring within the samples. Figure 4 shows the number of reactions involving a carbon atom removal in the four systems considered. At each temperature, the number of occurrences of such processes increases with initial functionalization density. Differently stated, the most probable reactions taking place in highly-oxidized samples involve the removal of a carbon from the sheet, consequently damaging the plane crystallinity. It is convenient to visualize the atomic sites of the underlying graphene plane from which carbon atoms have been removed through damaging chemical processes. The inset of figure 4 shows a map for the four systems at temperature T=1500 K where lighter points correspond to the removed atoms. As expected, most of the vacancies are formed in the originally oxidized areas, with more removals in the highly-oxidized samples. The original oxidized areas are strongly subject to damaging processes, which in the high-density cases can be completely emptied of carbon atoms. In those regions, the increasing number of vacancy formations could even give rise to the removal of the entire oxidized area leaving a hole. Such situation corresponds to a limiting case in which all the originally oxidized atoms are removed, i.e. with a minimum increase of sp 2 -hybridized carbons obtained through reduction. An image of the samples after annealing at temperature 1500 K is given in figure 5.
We proceed with the structural analysis of the reduced samples. By exploiting the atomistic picture developed in this work, we mapped the 'strain field' after thermal treatment. Specifically, we calculated the displacement of each carbon atom at the end of the reduction process with respect to its ideal position in an unoxidized graphene plane. In detail, the positions used are those extracted from the low-temperature system configurations obtained minimizing the total energy, after the process and in the pristine samples, respectively. Figure 6 shows the map of the 'strain field' for the four systems after reduction at T=1000 K: for each atom, a circle is depicted with a radius proportional to the defined displacement. It is clear how the strain field varies with the density of defects in the oxidized areas: the deviation from a pure graphene plane is indeed stronger in highly functionalized areas (which are characterized by larger regions with a darker colour). The regions involved correspond to the originally functionalized areas and the closer regions: they have lost oxygen groups as an effect of the thermal process but, as a consequence of the occurred carbon removals, a corresponding deformation of the system from the graphene plane is also gained. As expected, the chemical and structural modifications play a pivotal role in the transport performance of the samples: non-sp 2 domains and local changes in the graphene lattice constants represent detrimental scattering sources both for thermal and electronic transport.   We conclude now by discussing the dependence of such mechanical deformations as a function of the annealing temperature. Figure 7 shows the mean value of the atomic displacements in all the systems after the treatment at different temperatures. All the samples undergo a stronger mechanical deformation as the annealing temperature is increased with generally stronger effects in strongly oxidized samples. The larger energy which is transferred to the systems at high temperatures is then sufficient to desorb atoms from the plane but at the same time is sufficient to modify the local atomic configuration.

Morphology of functionalizing areas
We now focus on the spatial distribution of oxidizing agents on the pristine graphene plane. This is critical to understand the demonstrated action of the agglomerated and highly concentrated oxygen functionalities to improve the overall stability of the material.
The possible spatial configurations of defected areas within a GO sample is extremely large: in particular, an exceedingly large number of possibilities exist for the shape and dimension of the oxidized areas as well as for the local configuration of epoxy and hydroxyl groups within the same areas. This makes also the treatment of this problem unpractical via a systematic simulative approach. Hence, we adopt an alternative approach which allows to extract information on the effect of annealing on the morphology of GO. We built a set of GO systems which are characterized by different dimensions of the oxidized clusters and by different shapes of the  agglomerates. In order to focus on the effect of the spatial distribution of oxidizing impurities, they have been endowed with the same total initial carbon-oxygen ratio (about 0.26 in our case). In particular, the similar oxygen concentration will allow an almost fair comparison of the selected samples. The systems considered in this section are shown in figure 8. This set of GO samples contemplates different possibilities of agglomerated oxygen functionalizations: from a random distribution of oxygen groups in system 1 (with several and extremely small oxidized areas) to system 4 which represents the most clusterized system whose oxidized areas have linear dimensions of tenths of nanometers.
We begin calculating the amount of oxygen adsorbed onto the basal plane as a function of the annealing temperature: results are given in figure 9. After a weak decrease in the oxygen content at low temperatures, all the systems undergo a stronger reduction process for T>500 K, resulting in an efficient decrease of functionalizing   Figure 10 shows the evolution of this quantity normalized with respect to its initial value as a function of the temperature. A clear trend is now observed: for sufficiently high temperatures, less clusterized samples show a larger increase of the sp 2 carbons and hence a larger recovered cristallinity. These findings highlight the effect of GO morphology on the annealing process, which is found almost twice more effective in samples with diluted impurities. The images of the systems after annealing at temperature 1500 K offer an immediate picture of the different degree of purity of the system in the different cases. Figure 11 provides strong evidence that the release of oxygen groups is associated with two competing mechanisms: (i) the release of oxygen from the surface with a complete restoration of the sp 2 character of the carbon atoms and (ii) the release of oxygen from the surface with a corresponding removal of carbon atoms from the plane. The latter irreversibly affects the plane cristallinity by introducing topological defects like vacancies. These mechanisms generally occur in all the systems resulting in the similar O:C ratios of figure 9; the corresponding occurrencies are, instead, strongly dependent on the spatial distribution of the functionalized areas. Figure 12 clarifies this point by showing the occurrencies of carbon removals from the plane (normalized to the total number of carbons) in the different cases considered. A correspondigly larger number of carbon release is observed for more clusterized samples at any given  temperature, revealing how the presence of close oxygens favours reactions involving the extraction of carbons from the plane. In the inset of figure 12 the spatial distribution of the carbon removals is shown at T=1500 K in the four systems.
We conclude this section by analizing the structural features of the systems after thermal reduction. Figure 13 shows the strain field of the samples after a thermal treatment at T=1000 K. The deviation from the planar character of pristine graphene plane is indeed more pronounced for higher impurity concentrations. The regions involved correspond to the originally functionalized areas, which upon removal of oxygen groups are then subjected to relaxation process leading to a local rearrangement of the atomic configuration. The corresponding evolution of the average deviation as a function of the annealing temperature is shown in figure 14. We note how more clusterized samples yield a less planar final configuration after the treatment, with a deformation increasing as a function of the temperature. The trend observed is indeed common to all the systems. After an initial transient over which the atomic displacements change to a reduced extent, a sudden increase is found for a sufficiently high temperature. The threshold temperature depends on the initial morphology of the oxidized areas and shows a clear correlation with the onset of removal of carbon atoms from the basal plane ( figure 12). This points to the vacancy formation mechanism as the responsible force leading to a spatial rearrangement of the atoms in GO. Figure 13. Map of the strain field after thermal treatment at T=1000 K for the systems of figure 8. For each atom, a circle is depicted with a radius proportional to its displacement with respect to its ideal position in an unoxidized graphene plane.

Conclusions
In this work, we investigate, by means of classical molecular dynamics, the effect of thermal annealing on reducing different samples of GO. By monitoring the number of sp 2 carbon atoms and the oxygen-to-carbon ratio as a function of the annealing temperature, we study the effect of the density of oxidizing agents adsorbed and of the spatial distribution of the functionalized areas onto the efficacy of the thermal treatment.
We show how the initial concentration of oxygen-based functionalizing groups impacts the efficacy of the reduction process: for O:C densities larger than 35% the thermal reduction does not yield a sizeable recovery of the plane crystallinity. In these cases, the main mechanisms activated by thermal energy are reactions involving the removal of carbon atoms from the plane. As a result, vacancies are formed. A strong local rearrangement of the atomic configuration is obtained with a sizeable deformation of the basal plane whose regularity is consequently lost. Similarly, clusterization of the oxidized areas largely affects the quality of the reduced GO upon application of a thermal process. The amount of carbons going back to the unoxidized state after the treatment decreases for spatially concentrated functionalizations. At the same time, more vacancies are generally generated in such systems. Moreover, the thermal energy within the systems engenders a stronger deformation of the atomic configuration. In both cases, the spatial deformation of the basal plane following the thermal annealing increases with the annealing temperature. These findings suggest the possibility to improve the cristallinity of RGO through optimal choice of the parameters of the thermal reduction process, i.e. annealing temperature and initial GO morphology. In particular, for a given spatial distribution of oxidizing agents on the sample, an optimal value of the annealing temperature could be found, yielding a minimal distortion of the plane and a sizeable recovery of the sp 2 carbon domains. Finally, we mention that our simulation data could actually be further integrated into materials databases [44,45], for further use of machine learning algorithms [46][47][48][49] to extrapolate on morphologies and physical properties (electronic and thermal) of a large spectrum of reduced GOs morphologies. This could enable faster access to important information for designing composites with improved performances.