Molecular modelling of graphene nanoribbons on the effect of porosity and oxidation on the mechanical and thermal properties

Graphene is considered as the most promising nanomaterial of the recent decades given the huge amount of studies that have been performed to characterize its outstanding properties and in searching of novel applications. Following this tendency, this study covers the modelling of graphene nanoribbons (GNRs) with the aim of analyzing the effect of porosity and oxidation on the tensile mechanical properties and in-plane thermal conductivity through molecular dynamics (MD). Using quasi-static simulations the mechanical properties were evaluated in first place. A ‘hardening’ mechanism was observed for GNRs at porosities below 1%, i.e. perfect or near-perfect GNRs, by which the GNRs could withstand higher loading levels. This hardening effect was manifested in the carbon network by the generation of dislocation lines formed by pentagon-heptagon pairs (5–7 defects), which acted as a stress reliever. The failure of GNRs was produced as a tearing mechanism with cracks growing along the armchair or zigzag directions. The porosity affected all the analysed tensile mechanical properties (i.e., Young’s modulus, Poisson’s ratio, tensile strength and deformation at break), but with different tendency in the fracture properties due to the presence or absence of hardening behaviour in the GNRs. Nevertheless, the oxidation affected only the tensile modulus and Poisson’s ratio but not to the tensile strength and deformation at break. The thermal conductivity of the GNRs was affected either by the porosity and oxidation. Pores and oxidation groups acted as phonon scatterers since they disrupted the carbon network by the generation of vacancies or out-of-plane carbons, respectively, which decreased the phonon mean free path and thus the thermal conductivity. In conclusion, the porosity and oxidation of GNRs greatly determine the tensile mechanical properties and in-plane thermal conductivity of such materials and must be considered when tuning the synthetic pathways. The effect of porosity and oxidation on the tensile mechanical and thermal conductivity properties of graphene nanoribbons are evaluated through molecular dynamics simulations.


Introduction
Graphene is an one-atom thick layer of covalently bonded sp 2 -hybridized carbon atoms that form ideal hexagonal crystal lattices as honeycomb network [1]. This molecular structure originates a long-range pbonding and is responsible of the outstanding mechanical [2,3], electric [4,5] and thermal properties [6] of graphene. The applications and properties of graphene has been intensively studied in the last two decades, after the work published by Geim and Novoselov [7], and covers biomedical engineering [8], electronics [9,10], aerospace [11], among others. Nevertheless the oxidized forms of graphene, graphene oxide (GO) and reduced graphene oxide (rGO), are usually preferred over the neat structure since they are more active towards chemical modification and have higher interlayer gap than pristine graphene [12] due to the presence of polar groups (i.e., oxygen-containing groups). Improved electrochromic kinetics [13], lower resistance to heat transference through interface [14] and higher hydrophilicity [15] are reported for the oxidized varieties of graphene compared to the pristine form.
The synthesis of graphene oxide is usually carried out by treating graphene or graphite with strong acids and oxidizers, such as permanganate and sulfuric acid (Hummers method) [16,17], which delivers hydroxyl (C-OH) and oxirane (C-O-C) groups in the surface of the graphene sheet, together with carboxyl (C-COOH) and carbonyl (C-C = O) groups in the edge positions, as indicated the 13 C NMR technique [18,19]. The amount and dispersion of these oxygencontaining groups determines the extension of the unoxidized benzene ring and oxidized aliphatic sixmembered ring regions [20], which is identified as L-K structure [12]. Another important structural characteristic of the oxidized forms of graphene is the presence of structural defects in form of pores or holes due to overoxidation and lamellar peeling [21,22], which also affects the structure of the sheet. Li et al. [23] demonstrated through MD simulations that the presence of hole defects decreased the interfacial strength between polymer and graphene in nanocomposite systems. Thus, the properties and structure of the oxidized graphene derivatives depend on the amount and distribution of the polar groups and the hole defects. While pristine graphene is an hydrophobic material [24][25][26], the GO and rGO varieties are hydrophilic [27].
The modulation of properties through the amount of such oxygen-containing groups and pores in their structure is a way of tuning the graphene derivatives into the required application. Several methods have been described to synthesize graphene with controlled porosity such as the template approach [28,29] and the chemical/thermal activation [30,31]. Lin and Grossman [21] found that there existed a relationship between the size of nanopores and synthesis parameters, such as oxygen and epoxy concentration in the GO, as a potential application for the preparation of selective reverse osmosis membranes. The production of reduced graphene oxide (rGO) through the oxidative exfoliation-reduction of graphite yields a product with a certain degree of defects [1,32], in contrast with the chemical vapor deposition (CVD) method where the graphene is produced with low defects [33]. Achawi et al. [34] found that the structure of the graphene derivatives directly affected their surface reactivity and cytotoxicity, which revealed the importance of controlling the structure of the graphene materials depending on the desired final properties.
The experimental characterization of nanoparticles commonly is limited to the determination of structural parameters which includes composition, morphology, specific surface and size, among others [35]. Nevertheless, the determination of some properties such as the mechanical (e.g., Young's modulus and strength) or the thermal ones of individual nanoparticles is difficult to carry out experimentally and often requires techniques complex techniques such as atomic force microscopy (AFM) [2,[36][37][38] or Raman spectroscopy [39,40], but they generally lack of sufficient certainty and the results are very dependent on the experiment configuration and technician skills. Molecular simulations offer the best alternative to study the properties of nanoparticles [41] and they can also provide a powerful insight into the mechanism that underly complex phenomena such as nucleation of metallic nanoparticles [42,43] and biologic interactions [44][45][46]. In the case of graphene nanoparticles molecular dynamics (MD) simulations have been used to study a wide variety of properties such as the mechanical properties [47,48], thermal stability [49] and conductivity [50][51][52], interface properties [53,54], synthesis mechanisms [55] and reinforcing characteristics [56,57], among others.
In our previous study [14] MD simulations were conducted to compute the thermal conductivity across the interface, i.e. Kapitza resistance, between graphene and different polymeric matrices. Continuing in this line of work in the present study a series of molecular-based models of pristine graphene (G) and graphene oxide (GO) nanoribbons (GNR) with different percentage of pores (i.e. vacancy or hole defects) were developed focusing the characterization on the properties of isolated graphene. With this aim the effect of porosity and oxidation on the tensile mechanical properties and in-plane thermal conductivity were assessed by means of molecular dynamics (MD) simulations. Specifically the Young's modulus, Poisson's ratio, tensile strength and deformation at break were evaluated from quasistatic MD simulations, using the energy minimization approach, from which the failure mechanism was also assessed. Moreover, the thermal conductivity was calculated using the Mü ller-Plathe method [58], a non-equilibrium approach. Given that the thermal conductivity is fully dependant on the phonon transmission mechanism, the density of states (DOS) functions together with the mean free path (MFP) were also determined from MD simulations to help understanding the results found of thermal conductivity. The effect of porosity and oxidation on the graphene properties were assessed from a statistical point of view in terms of the F parameter (i.e., Fisher factor). This parameter is based on the analysis of the associated errors (i.e., variance) found in the results, in the context of an ANOVA analysis, and is a simple way to check for possible significative effect of factors in a series of results by comparing the variance of each series [59,60].

Modelling of graphene nanoribbons
A series of pristine graphene (G) and graphene oxide (GO) nanoribbon models having different percentage of pores were prepared. The carbon structure of the graphene sheets was generated using the VMD software [61], using dimensions of 250 Å length and 75 Å width, a C-C distance of 1.4579 Å which constituted 7200 carbon atoms and armchair structure in the edges. Then the oxide functionalities, i.e. hydroxyl and oxirane groups, were randomly bonded to the carbon atoms on both sides of the sheet. A total of 500 hydroxyl groups and 50 oxirane groups were added, which supposed a degree of oxidation (i.e. percentage of oxidized carbon atoms) of 8.3% within the range of experimental values [62]. The pores were generated by taking out a certain number of carbon atoms in the structures in random positions along the length of the nanoribbons. The morphology of the pore was irregular and its sizes varied between 10 and 25 Å of internal length. The carbon atoms at the edges (either at the sides of the nanoribbons or in the pore) were not passivated with hydrogen atoms, which according to Hu and Maroudas [63] doesn't have a significant effect on the thermal conductivity of graphene nanoribbons. The pore percentage was considered in a carbon atom basis, i.e. considering the number of carbon atoms removed in the defective structure respect to the total atoms in the original structure. Models having 0.5, 1, 2, 3, 4 and 6% of pores were prepared both for the pristine (G) and the oxide sheet (G). The carbon structure of the graphene sheets is shown in Fig. 1. The Tersoff potential [64] was employed to model covalent bonds between carbon atoms in the graphene sheet. As a bond order potential, the Tersoff potential allows for C-C bond formation or dissociation during the course of a simulation, which allows for example to study fracture properties. On the other hand the covalent bonds of the graphene with the oxide functionalities (i.e. hydroxyl and oxirane groups) were modelled through a Class II force field, whose parameters were based on the COMPASS version developed by H. Sun [65]. This type of force fields includes bond, angles and dihedrals interactions, but also cross-terms which increase the accuracy of the force field on the description of the properties of the system [66]. The non-bonded interactions were modelled by the Lennard-Jones 9-6 and Coulomb potentials for distances below the cutoff of 10.5 Å , while the Particle Mesh Ewald (PME) method was used to account for long electrostatic interactions. The partial atomic charges were derived according to bond increments as stated in COMPASS force field (see Fig. 2).
The LAMMPS simulation code [67] was used to perform the MD simulations. Periodic boundary conditions (PBC) were only applied in the x direction of the graphene sheets (i.e. in the 250 Å dimension of the sheet), while the y and z directions were kept non-periodic during the simulations. The graphene models (either G or GO) were equilibrated under the NP x T ensemble at 1 atm of pressure and 300 K of temperature during at least 0.5 ns until getting zerostress structures and reaching the plateau in potential energy.

Characterization of mechanical properties
The tensile mechanical properties of the graphene models were calculated from equilibrium MD simulations, using a quasi-static approach. The equilibrated graphene models, with zero stress, were initially minimized using the FIRE method [68] during 10,000 iterations. The minimized structures were then subjected to subsequent increments of deformation in the x direction of the cell box, controlled by the corresponding strain vector, while keeping the other cell parameters fixed. PBC were only applied to the x direction (direction of loading). After each deformation the structure were re-minimized under the FIRE method without any changes in the cell parameters. The strain was applied in two types of increments: the first 20 increments were of amplitude of 0.002, and the next increments were of 0.008. The response of the system to these increments was monitored as the virial stress calculated with the next expression: where V is the total volume of the graphene sheet, m k and u k denote the mass and velocity of the kth particle, respectively, r kl stands for the distance between kth and lth particles, and f lk is the force exerted on lth particle by kth particle. The Young modulus was calculated from the initial slope (i.e. the first derivative of the virial stress respect to the strain) of the stress-strain curve while the tensile strength was calculated as the highest stress registered in the simulation. The thickness of the graphene sheets was considered as 3.35 Å to calculate their corresponding volume, V. The in-plane Poisson's ratio, m xy , was calculated by monitoring the transversal length of the graphene sheet according to next equation: where e x and e y refers to the strain of the graphene in x (longitudinal) and y (transverse) direction, respectively. Given that the Poisson's ratio is strain dependant, an average value was compute from the plateau region between the strain range of [0, e C ], being e C the critical longitudinal strain where the Poisson's ratio changes from negative to positive [69].
In such cases where no e C was observed, due to premature fracture of the graphene, the Poisson's ratio was calculated up to fracture.

Characterization of thermal conductivity
The thermal conductivity j of the graphene models was calculated using the Mü ller-Plathe method [58], which is based on reverse non-equilibrium MD simulations (RNEMD). In this method the simulation box  is divided into slabs along the direction in which the thermal conductivity is being calculated, and a heat flux is imposed between the cold and hot slabs by exchanging the velocity vectors of the hottest atom of the cool slab with the coldest atom of the same mass in the hot slab. This process induces a gradient of temperature between the two slabs, from which thermal conductivity can be calculated as stated in the next formula: where J is the heat flux and dT/dx is the gradient of temperature in the system along x direction. The heat flux was calculated according to next formula: where t is the simulation time, A is the cross-sectional area of the graphene model, m is equal to the mass of the atoms transferred, and v hot and v cold are the velocities of the atoms exchanged. The graphene sheets were divided into a total of 20 slabs of about 12.6 Å , the cold slabs were placed in the ends while the hot region was placed in the middle of the graphene sheet. Initially the systems were equilibrated in the canonical ensemble (NVT) during 100 ps, and then the Mü ller-Plathe algorithm was applied during 1 ns in the microcanonical ensemble (NVE) after getting a steady temperature gradient.

Characterization of in-plane phonon density of states (DOS)
The thermal conductivity of a system is highly related with the modes and frequencies of vibration of its constitutive atoms, i.e. the phonons. The phonon density of states (DOS), which describes the proportion of states that are occupied by the system at each energy in the range of modes of vibrations, was calculated for each graphene nanoribbon. The DOS was computed by applying the Fourier transformation to the velocity autocorrelation functions (VACF) as stated next: where v(0) and v(t) are the velocities of each atom at initial time and at time t, respectively. The VACF functions were acquired along 1 ps trajectories with five repetitions for statistical average. Only the x and y components of the velocity of atoms were considered for the calculation to get in-plane DOS, since the z component doesn't contribute to thermal conductivity of graphene nanoribbons due to their planar or sheet geometry.

Statistical analysis
With the aim of analysing the results from a statistical point of view, the statistical F-value was calculated according to next formula: where s series is the standard deviation associated to a series of results in which an associated factor could be acting (i.e. the effect of oxidation or the effect of pores) and s error is the standard deviation associated to the MD method itself due to random errors. The Fvalue compares if two standard deviations (i.e., variances) are statistically equal or not, and thus the existence of a factor which produces variations in the results in addition to random errors. If the null hypothesis H 0 was true, then the standard deviation in the results would only be associated to random errors (s series = s error ). On the other hand, if the alternative hypothesis H 1 was true then there would be an additional factor which produced the differences found in the results in addition to random errors (s series = s factor ? s error ). The F-value was compared to a critical F-value, considering a statistical significance of 95%, in order to accept or reject the null hypothesis H 0 .

Results and discussion
Mechanism of failure under tensile stress of graphene nanoribbons The stress-strain curves obtained during the tensile simulations of the graphene models are plotted in Fig. 3 and Fig. 4 for pristine graphene (G) and graphene oxide (GO), respectively, with a comparison between the different pore percentages. It was observed that the stress-strain curves of some of the models, depending on the percentage of pores, exhibited a complex behaviour with two steps of deformation in which after the initial linear elastic response and plastic yield there was a sharp increase of stress at deformations around 0.28. This 'hardening' behaviour [70] associated to high strain was produced in the pristine graphene (G) models of 0 and 0.5% pores (Fig. 3), whereas in the graphene oxide (GO) it was also observed for the 1% pore model (Fig. 4). When analysing the snapshots from MD simulations in these models it was observed that, around 0.28 of deformation, it was produced the reorganization of some atoms in form of dislocations. In Fig. 5 a series of snapshots from the tensile simulation of the GO model with 0.5% pores have been selected to show a graphic schematic. The dislocations appeared in the structure as pentagonheptagon pairs due to the breakage of one C-C bond and later displacement of the atoms. This pentagonheptagon pairs have also been identified in the literature as '5-7' defects [71] or Stone-Wales defects [23]. It was observed that these dislocations started preferably at the borders of the pores, which acted someway as nucleation points of this type of defects. While increasing the deformation of the sheet, new pentagon-heptagon pairs were created next to the previous ones and the dislocation grew in a straight line which followed the armchair or zigzag directions of the crystal lattice of the graphene sheet. In Fig. 6 a snapshot showing the formation of a dislocation line with detail is plotted, while in Fig. 7 the stress map of snapshots before dislocation of atoms, during the generation of a dislocation line and after cracking are represented. As it can be observed from the stress maps, the dislocation of atoms decreased the stress of the atoms just besides the dislocation line, i.e. the dislocation acted as a stress reliever, and helped the GNR sheet to accommodate higher degrees of deformation before failure. Other authors such as Mao et al. [72] and Liang and Huan [73] also reported this 'hardening' behaviour, while the presence of such dislocations has also been demonstrated experimentally using sophisticated TEM microscopy [74,75]. When the pore percentage was higher than the one indicated previously for each type of graphene (i.e., higher than 0.5% for pristine GNRs and higher than 1.0% for GO), it was observed that the graphene sheet broke at lower deformations without the previous formation of such dislocation defects. Thus, no 'hardening' behavior was observed in the stressstrain curves for these GNRs. At a certain stress, the structure finally collapsed and tore with even the production of monoatomic carbon chains at certain points of the crack [76]. It was also observed that the crack growth was predominantly produced along straight lines aligned with the armchair or zigzag directions of the hexagonal crystal lattice of graphene (Fig. 8), which was in agreement with other studies [77,78]. These cracks were produced in form of tears (i.e. tearing mechanism) which occurred in the direction perpendicular to the direction of loading. In some cases, the direction of the tear occurred with some tilt respect to the direction of loading, and it      was observed that even it occasionally changed the direction from armchair to zigzag (or vice versa). This mechanism of deformation and failure of graphene sheets, through the formation of dislocations and ripping, has also been observed for carbon nanotubes [71].

Tensile mechanical properties of GNRs
A summary of the main tensile mechanical properties of the pristine graphene (G) and graphene oxide (GO) models in terms of Young's modulus, tensile strength, deformation at break and Poisson's ratio is shown in Table 1 and Table 2. Graphic representations of these properties versus the porosity are plotted in Fig. 9 and Fig. 10. The Young modulus of pristine graphene (G) without pores (0% model) was calculated as 1.17 TPa, which was near the value of 1 TPa reported experimentally [2] and by ab initio simulations [79]. Nevertheless, the intrinsic strength (i.e. tensile strength) was overestimated for this model (238 GPa) in comparison with that reported in the literature of 130 GPa [2] or even lower [80], which would be closer to the values of tensile strength obtained for the models with pores. In the same way, the values of deformation at break were higher than those reported in the literature [79], which indicated that the MD method tended to overestimate the value of properties at high levels of deformation. Liang et al. [81], which conducted a similar study on pristine graphene sheets using the AIREBO potential, obtain equivalent results either in terms of the stressstrain curve and tensile strength. In the case of the graphene oxide (GO), the reported values for Young's modulus are dependant on the degree and type of oxidation (i.e., oxirane and/or hydroxyl). Kang et al. [82] obtained from experimental measurements a Young's modulus of about 700 GPa for a graphene sheet with 25% oxidation, while Suk et al. reported a modulus of 207 GPa for a graphene sheet of 40% oxidation. On the other hand, Khoei et al. [83], who also studied the mechanical properties of graphene oxide through MD simulations, found a Young's modulus of around 950 GPa for a graphene sheet with 10% degree of oxidation. In the present study a value of 903 GPa was obtained for the 0% porosity model and 8.3% degree of oxidation, which was consistent with the results found in the literature.
According to previous studies, the Poisson's ratio depends on the size and chirality of the sheet [84], but also on the applied strain [69]. The values of Poisson's ratio reported in Table 1 and Table 2 were calculated in the range of deformation [0, e C ] as explained before, and were negative in all the cases, either G or GO, indicating thus an auxetic behaviour for these nanoribbons [85,86]. According to the study by Jiang and Park [69], the negative Poisson's ratio phenomenon can be observed in graphene sheets with widths up to about 10 nm, which agrees with the results found in the present study. The 'hardening' mechanism explained previously also affected the transversal deformation of the graphene nanorribons and thus to the value of the associated Poisson's ratio. It was observed that at deformation levels below 0.28 the Poisson's ratio had negative and approximately constant values but, at higher deformations (i.e., in those cases where no failure was produced before) the Poisson's ratio tended to increase linearly and even turned to positive values (see Fig. 10). Besides, a graphic comparative of the snapshots between the longitudinal section of the graphene nanoribbon at different deformation states have been plotted in Fig. 11 to illustrate this behaviour. The values of the parameter F were calculated from the standard deviation associated to the results to analyse, from a statistical point of view, the effect of porosity and oxidation of the graphene nanoribbons on its tensile mechanical properties. The results are summarized in Table 3, in which bold values indicate that the value of F was higher than the critical value, and thus the alternative hypothesis H 1 was accepted, i.e., there was an effect of the factor. The first analysed factor, i.e. porosity, had a significative effect in all the assessed tensile mechanical properties and either the Young modulus, tensile strength, deformation at break or Poisson's ratio were statistically affected due to the presence of pores in the graphene sheet. The tendencies versus porosity found on these properties (see Fig. 9) were similar in both types of graphene nanoribbons (G and GO), although different for each analysed property. In the case of the Young's modulus it decreased approximately in a linear fashion when the porosity increased, up to a decrease of around 16-17% when the porosity was 6%. The tendency found in Poisson's ratio also   followed a linear trend, although for this property the porosity produced an increase in its value, reaching up to 70% for a 6% porosity respect to the non porosity nanoribbon. For the tensile strength and deformation at break, both properties related with fracture parameters, the tendency followed a exponential decrease. At porosities below 2% there was produced a dramatic decrease either in the tensile strength and deformation at break (with maximum decreases of about 54-51%), but at higher porosities it was observed an asymptotic behaviour where the porosity had a lower impact on the values of such properties. This different tendency versus porosity found in one hand in the Young's modulus and Poisson's ratio, and on the other hand on the tensile strength and deformation at break could be explained by the hardening mechanism described before. The hardening behaviour, by which the graphene nanoribbons could afford higher stress and deformations (see Fig. 3 and Fig. 4), was only produced in perfect (i.e., 0% porosity) or practically perfect (i.e., below 1% porosity) graphene nanoribbons, but at porosities higher than 1% this mechanism was not produced. The absence of this hardening mechanism in the graphene nanoribbons produced the sharp decrease observed in the properties related with fracture, i.e. tensile strength and deformation at break. The Young's modulus and Poisson's ratio were not affected by this mechanism since these properties are calculated in the elastic region.
Regarding the second analysed factor, i.e. the oxidation of the graphene nanoribbons, the F values indicated that only the Young's modulus and Poisson's ratio were significatively affected, whereas the tensile strength and deformation at break were not modified due to this factor. As can be seen in Fig. 9, the oxidation decreased the values of Young's Figure 11 Snapshots showing a comparative of the longitudinal section of the 0% pores GO nanoribbon during the tensile deformation at different deformation states. In bold: H 0 is refused (significative effect of factor)/In italics: H 0 is accepted (no effect of factor) modulus by around 240 GPa (around 20-24%) respect to the G nanoribbons, which was in agreement with a 8.3% degree of oxidation or functionalization [87]. This decrease in the rigidity of the graphene sheet due to the oxidation can be explained considering the disruption of the sp 2 network due to the presence of hydroxyls and oxirane groups [88]. As was described by Zheng et al. [87], the sp 3 hybridization of carbon atoms breaks the local p bond of the graphene sheet and also makes the sp 3 with an off-plane structure easier to be bent by tension, which weaken the graphene structure. The fact that the oxidation has no significant effect either in the tensile strength and deformation at break, i.e. tensile fracture properties, indicated that at high deformations the effect of oxidation was negligible. Zandiatashbar et al. [89], who also studied the effect of defects on the mechanical properties of graphene, found that the porosity govern the mechanical properties over the sp 3 defects (i.e. oxidation) in the vacancy-defect regime. In the present study it was observed that this statement applied essentially to the fracture properties, but in the case of tensile mechanical properties in the elastic region (i.e. at low strains) the effect of oxidation was comparable to the effect of porosity and both factors contributed significatively to the Young's modulus and Poisson's ratio.

Thermal conductivity of graphene nanoribbons
A summary of the results of thermal conductivity (j) measured in the pristine (G) and graphene oxide (GO) nanoribbons can be found in Tables 4 and 5, respectively. The thermal conductivity of non-porous pristine graphene (G) nanoribbon (0% model) was calculated as 160.5 W/mÁK, according to the RNEMD simulations performed in this study. This value was in agreement with the values reported for armchair edge graphene nanoribbons in other MD studies [90,91] as well with the experimental range found for suspended graphene of 100-300 W/mÁK [92]. The thermal conductivity of graphene nanoribbons is usually reported an order of magnitude lower than that of graphene [93], which is a consequence of the edge and defect scattering due to its narrow width (7.5 nm width in the present study) [94]. Besides, Li and Zhang [95] demonstrated that the Tersoff potential, i.e., the potential used in the present study, produced similar heat flux distribution along graphene nanoribbons than the AIREBO potential.
The statistical analysis results made in terms of parameter F are summarized in Table 6. It was found that either the porosity or oxidation had a significative effect on the calculated thermal conductivity, since the determined parameter F was higher than the critical value in each case. The curves of thermal conductivity versus porosity plotted in Fig. 12 showed a linear trend in both types of nanoribbons, although two regions with different slope were observed in both types of graphene nanoribbons. Below 1% porosity, the thermal conductivity exhibited a higher slope versus porosity, but at higher porosities, the curve had a lower slope, i.e. a lower Table 4 Thermal conductivity (300 K) for pristine graphene (G) nanoribbons (variation percentages are in brackets)

Porosity (%)
Thermal conductivity j (W/m K) 0 160.5 ± 11.8 0.5 142.6 ± 2.2 (-11%) 1 125.4 ± 1.4 (-22%) 2 111.1 ± 1.6 (-31%) 3 103.6 ± 2.2 (-35%) 4 95.0 ± 2.3 (-41%) 6 77.0 ± 0.9 (-52%) Table 5 Thermal conductivity (300 K) for graphene oxide (GO) nanoribbons (variation percentages are in brackets) Porosity (%) Thermal conductivity j (W/mÁK) 0 43.6 ± 2.8 0. 5 40.2 ± 0.8 (-8%) 1 38.8 ± 0.5 (-11%) 2 36.6 ± 0.8 (-16%) 3 34.1 ± 0.9 (-22%) 4 32.3 ± 1.5 (-26%) 6 28.9 ± 0.2 (-34%) is accepted (no effect of factor) dependence of porosity. This double-side linear behaviour of thermal conductivity versus porosity was also obtained by other authors [96,97] for graphene sheets. In addition, other studies found an exponential decay of thermal conductivity when increasing porosity in pristine GNR [98,99]. Given that the thermal conductivity is based on the transmission of phonons through the net of carbon atoms, the presence of pores acted as phonon scattering points where the transmission of vibrations was cancelled. Nevertheless, the scattering effect is more significant in small pores (i.e., lower porosity) than in longer pores (i.e., higher porosity) [97,100]. The analysis of the phonons and its connection with thermal conductivity is discussed in next section. The oxidation of GNRs approximately decreased their thermal conductivity for about 60-70% respect to that found for the pristine GNR, indicating thus that GO nanoribbons were worst heat conductors than their non-oxidized version (i.e., are better insulators). This sharp decrease of thermal conductivity was in the order as calculated by other studies for a 8.3% of oxidation of the graphene sheet [91,101]. The oxidation of GNR through the introduction of hydroxyl and oxirane groups produced a change in the hybridization of the implied carbon atoms from sp 2 to sp 3 , which generated structural or morphological changes in the graphene sheet. These structural changes were due to variations found in the C-C bond distances and C-C-C-C dihedral angles produced in the environment of the oxidation groups, as can be seen from the representations shown in Figs. 13 and 14. The C-C bond distance in a non-oxidated environment (i.e., in pristine GNR) was of 1.46 Å (Fig. 13a), but this distance increased in the first neighbour atoms of the C atoms involved in the oxidation groups. When a hydroxyl or oxirane group was attached to the C atoms, the C-C bond distance of the first neighbours increased to around 1.51 Å (Fig. 13b) and 1.48 Å (Fig. 13c), respectively, or even to 1.52 Å in the C atoms directly bonded to the oxirane group. Consequently, the C atoms bonded to the oxidative groups had to get out of the plane to accommodate these bond increments, as were manifested in the dihedral angle values (Fig. 14). The alteration of the structure of the GNR sheet, either by an increase of the bond distance or by a change in the dihedral angles, generated points were a phonondefect scattering mechanism was produced, which finally manifested in a lower thermal conductivity. Mu et al. [101] also indicated that the change in mass that produced the oxidation also perturbed the vibration of the C atoms, but the change produced in the structure affected more significatively to the thermal conductivity than the change in mass. On the other hand, Lin and Buehler [91] explained that the decrease found in the in-plane thermal conductivity after the oxidation of GNR was also due by a reduction of the phonon mean free path, as it is discussed in the next section.

In-plane phonon density of states (DOS)
According to Balandin et al. [102], acoustic phonons are the main heat carriers in carbon materials such as graphene, thanks to the exceptionally strong sp 2 covalent C-C bonds of its lattice structure. In this way, the results of thermal conductivity can be explained from the analysis of in-plane phonon density of states (DOS), which were obtained from the Fourier transformation of the atomic velocity autocorrelation functions (VACF). According to the kinetic theory [91,101], the thermal conductivity can be related with phonons through the next expression: where C i is the volumetric specific heat of the vibration mode, v i the group velocity, l i the mean free path of the phonon mode i (with an angular frequency x), and l eff the effective phonon mean free path. The previous relation establishes that the thermal conductivity is proportional to the phonon mean free path, the heat capacity and the phonon velocity. On the other hand, according to the Debye model, the specific heat capacity C i is related with the phonons as: where h is the Planck constant, D(x) is the DOS for phonon mode i with frequency x, T is the absolute temperature and f(T,x) is the Bose-Einstein distribution function. Given that the derivative term is constant for a fixed phonon mode i with frequency x at constant T [91], Eq. 8 states that the specific heat capacity C i of a phonon mode is mainly dependant on the product of frequency and DOS function.
Considering thus the previous relationships, the DOS multiplied by the frequency spectra and the phonon mean free path were calculated to get a better understanding of the effect of the porosity and oxidation on these parameters and therefore in the conductivity. The phonon DOS spectra are plotted in Fig. S15 (see Supplementary Material) whereas the product of frequency and DOS function spectra are plotted in Fig. 15. In first place, the porosity had negligible effect on the DOS spectra of either G and GO nanoribbons, which indicated that the presence of pores did not affected the vibrational modes of the graphene. Nevertheless, when comparing the DOS spectra of G and GO (see Fig. 16), the oxidation produced the apparition of a new in-plane vibrational mode at high frequencies, around 120 THz, associated to the C-O bonds. According to Eq. 8,  these high frequency phonons contributed increasing the in-plane thermal conductivity of the GO nanoribbons, nevertheless their contribution was not enough since the thermal conductivity of the GO structures was lower than the G structures as was described in previous section. At this point, the modification of the mean free path of the phonons due to the oxidation played a major role in the thermal conductivity of the nanoribbons than the vibration frequencies. The effective phonon mean free path l eff can be expressed through the Matthiessen's rule as: where l ph-ph is the intrinsic mean free path due to the phonon-phonon scattering and l ph-b is the phonon Figure 15 In-plane phonon DOS multiplied by the frequency of a pristine graphene (G) and b graphene oxide (GO) nanoribbons.
free path due to boundary scattering, which is the nanoribbon length, L (i.e., l ph-b = L). Combining Eqs. (7) and (9) the thermal conductivity can be related with the nanoribbon length as: According to Eq. (10), from the linear plot of 1/j against 1/L the thermal conductivity for an infinitely long nanoribbon can be obtained from the intercept of the extrapolation whereas l ph-ph from the slope of the fitting. With the aim of calculating the intrinsic mean free path l ph-ph at each porosity, a series of models were generated by replicating twice or three times the models presented in Fig. 1 along the longitudinal direction to obtain doubled (2L) and tripled length (3L) graphene nanoribbons. After their equilibration, the thermal conductivity of each model was calculated according to the Mü ller-Plathe method. The individual results of thermal conductivity for these 2L and 3L models can be found in the Supplementary Material. From the linear plots of 1/j against 1/L the intrinsic mean free path l ph-ph was calculated for each porosity level on each type of graphene nanoribbon, G and GO, and the results are summarized in Fig. 17. When analysing in first place the effect of porosity, in the case of G nanoribbons it was observed that l ph-ph decay rapidly with the porosity from 2000 Å at 0% porosity to around 240 Å at 3% porosity, from which l ph-ph maintained constant with no further decrease. For the GO nanoribbons there were observed two regions of l ph-ph versus porosity: between 0 and 1% the value of l ph-ph was around 180 Å , while at porosities higher than 2% the l ph-ph decreased and took values around 115 Å . In comparison with G nanoribbons, the effect of porosity had less impact on GO nanoribbons. On the other hand, it can also be observed that l ph-ph was an order of magnitude higher in G nanorribons compared to GO nanoribbons at porosities below 1%, and approximately twofold higher in G nanorribons at porosities over 2%. This higher values of l ph-ph for G nanoribbons would explained their higher thermal conductivity and the effect of oxidation was mainly manifested in a sharp decrease of the l ph-ph parameter. Although was not the purpose of the present study to analyse different oxidation degree levels, Chen and Li [103] determined that the degree of oxidation was directly proportional to the inverse of the phonon mean free path. In summary, either the porosity and the oxidation affected mainly the phonon mean free path, and either the pores and the oxidized groups (i.e., hydroxyl and oxirane) acted as scattering points for the transmission of the phonons throughout the graphene sheet.

Conclusions
A series of molecular simulation models of graphene nanoribbons (GNRs) were constructed with different degrees of porosity and with or without oxidation, i.e. neat (G) or oxide (GO), to assess the effect of such parameters in the tensile mechanical and thermal conductivity properties. In first place, regarding the failure mechanism of graphene nanoribbons under tensile loadings, it was observed that at low porosity levels and around 0.28 deformations a 'hardening' behaviour was manifested in the structures as dislocations lines formed by pentagon-heptagon pairs (i.e., 5-7 defects) which helped the nanoribbons to accommodate higher tensile stress and acted as stress reliever. Such GNRs exhibited a complex stressstrain curve with two steps of deformation. This 'hardening' behaviour required of near perfect or non-defective graphene nanoribbons, so that at porosities higher than around 1% the nanoribbons broke at deformations below 0.28 with no hardening effect. The breaking of the nanoribbons was produced as a tearing mechanism in which the crack growth predominantly along straight lines perpendicular to the direction of loading and aligned with the armchair Figure 17 Plots of intrinsic mean free path l ph-ph versus porosity.
or zigzag directions of the hexagonal crystal lattice of graphene. The porosity affected all the analysed tensile mechanical properties (i.e., Young's modulus, tensile strength, deformation at break and Poisson's ratio), but with different tendency in the fracture properties (tensile strength and deformation at break) which was explained by the presence or absence of hardening behaviour at high deformations in the graphene nanoribbons. Nevertheless, the oxidation affected only the tensile and Poisson's ratio but not to the tensile strength and deformation at break. On the other hand, the thermal conductivity of the GNRs was affected either by the porosity and oxidation according to the F analysis. Pores and oxidation groups (i.e., hydroxyl and oxirane) acted as phonon scatterers since they disrupted the carbon network by the generation of vacancies or out-of-plane carbons, respectively, which decreased the phonon mean free path and thus the thermal conductivity. In conclusion, the porosity and oxidation of GNRs must be considered when tuning the tensile mechanical properties and thermal conductivity of such materials through the synthetic pathways. Moreover, it can also be established that the MD methodology is a valid technique to compute the mechanical and thermal properties of graphene derivatives.

Funding
Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.

Data and code availability
The raw data is available upon request to the authors.

Declarations
Conflicts of interests The authors have no relevant financial or non-financial interests to disclose.
Ethical approval Not Applicable.

Open Access
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://c reativecommons.org/licenses/by/4.0/.