Thermal Conductivity and Mechanical Properties of Nitrogenated Holey Graphene

TL;DR

This study investigates the mechanical and thermal properties of nitrogenated holey graphene using first principles calculations and molecular dynamics, revealing high elastic modulus, tensile strength, and moderate thermal conductivity that decrease with temperature.

Contribution

It provides the first detailed computational analysis of nitrogenated holey graphene's mechanical and heat transport properties, guiding future experimental research.

Findings

Elastic modulus of 335 GPa at room temperature

Tensile strength of 60 GPa at room temperature

Thermal conductivity of 64.8 W/m-K at 300 K

Abstract

Nitrogenated holey graphene (NHG), a two-dimensional graphene-derived material with a C2N stoichiometry and evenly distributed holes and nitrogen atoms in its basal plane, has recently been synthesized. We performed first principles calculations and molecular dynamics simulations to investigate mechanical and heat transport properties of this novel two-dimensional material at various temperatures. First principles calculations based on density functional theory yield an elastic modulus of 400 +/- 5 GPa at 0 K, 10% larger than predicted by molecular dynamics simulations at low temperatures. We observed an overall decreasing trend in elastic modulus and tensile strength as temperature increases. At room temperature, we found that NHG can present a remarkable elastic modulus of 335 +/- 5 GPa and tensile strength of 60 GPa. We also investigated the thermal conductivity of NHG via non-equilibrium molecular dynamics simulations. At 300 K an intrinsic thermal conductivity of 64.8 W/m-K was found, with an effective phonon mean free path of 34.0 nm, both of which are smaller than respective values for graphene, and decrease with temperature. Our modeling-based predictions should serve as guide to experiments concerning physical properties of this novel material.

Tables

Table 1. Table 1: Intrinsic thermal conductivity and effective phonon mean free path for NHG.

T (K)	${\kappa }_{\infty }$ (W/m-K)	$\Lambda$ (nm)
300	64.8	34.0
450	44.7	23.8
600	39.4	21.8

Equations

$$T_i=\frac{2}{3N_i{k}_B}\sum _j\frac{p_j^2}{2m_j},$$

$${\kappa }_L=\frac{⟨J_x⟩}{⟨dT/dx⟩}.$$

$$\nu =\frac{{\sigma }_t}{{\sigma }_l}.$$

$$E=\frac{{\sigma }_l-\nu {\sigma }_t}{ϵ},$$

$$\frac{1}{{\kappa }_L}=\frac{1}{{\kappa }_{\infty }}\left(1+\frac{\Lambda }{L}\right),$$

Full text

††journal: Carbon

1 Introduction

Since graphene’s first reported observation [1, 2], the family of two-dimensional (2D), atomically thin materials has grown considerably [3]. Materials such as silicene[4], germanene [5], phosphorene [6], stanene [7] and borophene [8], to name just a few, have either been produced or isolated in single-layer or few-layer form. Together with graphene, these materials present a complete range of electronic properties, including metals, semiconductors and insulators. Regarding mechanical properties, most of these materials are quite resistant to tensile strain while presenting almost no resistance to compressive strain, due to their sheet-like structure. Meanwhile, concerning heat transport properties, these 2D wonder-materials, also cover a wide range of values from low to very high.

Recently, yet another 2D material has been synthesized via direct functionalization of graphene sheets  [9]. This novel graphene-derived material presents a crystal structure with evenly distributed holes and nitrogen atoms with a C2N stoichiometry in its basal plane. The atomic structure of nitrogenated holey graphene (NHG) consists of rings of carbon atoms terminated by nitrogen atoms, as show in Fig. 1. In this work we employed extensive first principles calculations and classical molecular dynamics simulations in order to predict mechanical and heat transport properties of this novel material. To the best of our knowledge there are no experimental nor theoretical investigations of such properties yet, therefore predictions based on our modeling results can serve as a guide to groups interested in the physical properties of this novel 2D material.

2 Methods

We have performed atomistic modeling of NHG via classical molecular dynamics simulations (MD) where the Tersoff bond-order potential was employed to describe interatomic forces [10, 11]. Although the Tersoff potential has not been parameterized to model NHG specifically, it has already been parameterized for graphene [12], hexagonal boron nitride (hBN) [13], and graphene–hBN structures [14]. Therefore, in the present work we employed potential parameters from Lindsay and Broido [12, 13] for C-C interaction and obtained the parameters for C-N via the usual mixing rules for the Tersoff potential [11].

In order to verify the accuracy of the chosen potential parameters, we first obtained phonon dispersion relations of NHG by diagonalizing the dynamical matrix via the General Utility Lattice Program (GULP) [15, 16]. The phonon dispersion presented in Fig. 2 was obtained for an 18–atom unit cell (shown in Fig. 1). From the dispersion curves we see that NHG presents three acoustic modes: two of the modes present linear dispersion around the $\Gamma$-point, while the third mode has a quadratic dispersion, similar to graphene and hBN [12, 13]. The absence of imaginary frequencies in Fig.  2 indicates that the crystal structure is stable when described by the chosen potential parameters, at least at [math] K. Furthermore, the features observed in Fig. 2 are, at least qualitatively, in agreement with recent ab-initio calculations based on density functional theory (DFT), which also confirmed the stability of NHG’s predicted crystal structure [17].

We also investigated the mechanical response of NHG sheets by MD simulations and first principles DFT calculations. MD simulations were performed with the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)[18] package, with periodic boundary conditions along in-plane directions and free boundary conditions in the perpendicular direction. Mechanical properties were evaluated by performing uniaxial tensile strain along in-plane directions. Each simulation cell was equilibrated at the desired temperature via a Nosé-Hoover thermostat[19, 20], while the equations of motion were integrated with a $0.25$ fs timestep. Each cell was relaxed to zero stress along planar directions at the respective temperature via coupling to a thermostat and a barostat. Uniaxial strain was applied by increasing the length of the simulation cell along the loading direction at a constant engineering strain rate of $1\times 10^{8}$ s*-1*. Even though the atomic positions were not explicitly rescaled in the simulation, the small strain rate prevented the formation of voids, and the atoms moved following their natural dynamics. The width of the simulation cell, i.e. the length perpendicular to the loading direction, was allowed to relax in order to ensure zero average stress along that direction. Virial stresses were calculated at each strain level to obtain engineering stress-strain response.

DFT calculations were performed as implemented in the Vienna ab initio simulation package (VASP)[21] using the Perdew-Burke-Enzerhof (PBE) generalized gradient approximation exchange-correlation functional [22]. The projector augmented wave method was employed with an energy cutoff of 400 eV [23], and a k-point mesh size of $7\times 7\times 1$ in the Monkhorst-Pack mesh was employed for Brillouin zone sampling [24]. A supercell consisting of 72 atoms was fully relaxed with geometry optimization and periodic boundary conditions by the conjugate gradient method. For smearing, the tetrahedron method was applied with a width of $0.1$ eV [25]. After obtaining the optimized structure, we increased the periodic simulation box size in multiple steps along the loading direction with a small engineering strain of 0.0025. In this case to avoid an instant void formation, the atomic positions were rescaled with respect to the box size such that a gap is not formed in the simulation box. A full energy minimization was performed at each step. We studied the mechanical properties under biaxial and uniaxial loading conditions. In the case of biaxial loading, the box size along the perpendicular direction remained constant which resulted in the biaxial stress condition. On the other hand, we allowed the change of simulation box size in perpendicular direction to reach negligible stress to ensure uniaxial stress conditions. The stresses on the minimized structures were measured and plotted against the strain in order to obtain the elastic properties. For the density of state calculation the smearing width was set to $0.01$ eV and a denser k-point mesh size of $11\times 11\times 1$ was used. For the ab initio molecular dynamics simulations (AIMD), a Langevin thermostat was used for maintaining the temperature and the integration time step was set at $1.0$ fs.

The thermal conductivity of NHG was calculated via non-equilibrium molecular dynamics (NEMD) simulations with periodic boundary conditions, in order to investigate size (length) and temperature dependence of heat conduction in this material. Within this methodology, a temperature gradient is imposed in the system and the resulting heat flux can be measured, such that the thermal conductivity is directly calculated from Fourier law [26]. During NEMD simulations, the entire system was initially equilibrated at the desired temperature via a Nosé-Hoover thermostat for $100$ ps, while the equations of motion were integrated with a $0.25$ fs timestep, excluding fixed atoms at the two ends, as shown in Fig. 3. Then the thermostat was turned off, and the simulation cell was divided in $21$ slabs along the heat transport direction. The first and last slabs were assigned to be the hot and cold reservoirs, respectively, and were coupled to independent thermostats, set at different temperatures, while the intermediate slabs were not directly coupled to any thermostats. The temperature of each slab can be computed from the equipartition theorem as

[TABLE]

where $T_{i}$ is the temperature of $i$-th slab, $N_{i}$ is the number of atoms in $i$-th slab, $k_{B}$ is Boltzmann’s constant, $m_{j}$ and $p_{j}$ are atomic mass and momentum of atom $j$, respectively. At each simulation step a given amount of energy was added to the atoms in the hot reservoir by the thermostat and an equivalent amount of energy was removed from the atoms in the cold reservoir by the other thermostat, such that on average the total energy exchanged between the system and the thermostats is zero, and the average system temperature is stable. The heat flux is given by the difference between the energy added and the energy removed from the system. After an initial transient time of approximately $200$ ps the heat flux reaches a stationary state and a stable temperature gradient is established along the simulation cell. At this stage, the simulation was carried out for at least another $2.0$ ns to calculate the thermal conductivity from the average heat flux and the average temperature gradient.

Fig. 4(a) shows the energy added to the hot slab and removed from the cold slab as a function of time. The net amount of energy exchanged with the reservoirs averages to zero over time, and thus the total energy is conserved and the system is in a stationary state. In addition, the linear dependence with time, shows that energy is exchanged at constant rate. The heat flux along the $x$-direction is given by $J_{x}=\frac{dE/dt}{A}$, where $A$ is the cross-sectional area of the sample given by its width multiplied by thickness. We assumed a thickness of $0.33$ nm for nitrogenated holey graphene which is based on the thickness of graphitic carbon nitride  [27, 28]. Fig. 4(b) shows the time-averaged temperature at each slab. In this case, by neglecting the initial jumps adjacent to the hot and cold slabs, a linear temperature gradient ($dT/dx$) is established along the sample. Finally, the thermal conductivity of a sample with size $L$, $\kappa_{L}$ was computed from the ratio between the average heat flux in a given direction $\langle J_{x}\rangle$ and the average temperature gradient along the same direction

[TABLE]

3 Results and discussion

3.1 Mechanical properties

In Fig. 5, the DFT results for stress-strain response of NHG under uniaxial and biaxial loading is shown. As described above, unidirectional strain was applied along the longitudinal direction while on the transverse direction the box size was fixed. In this case, the stress along the longitudinal ($\sigma_{l}$) and transverse ($\sigma_{t}$) directions were calculated at each strain level in order to evaluate elastic properties. As shown in Fig. 5, the stress presents a linear dependence with the strain, which confirms that the specimen is stretched within its elastic limit. The Poisson ratio $\nu$, can be obtained using the following ratio:

[TABLE]

Finally, the elastic modulus $E$, can be obtained by:

[TABLE]

where $\epsilon$ is the applied strain. Based on our DFT calculations, the elastic modulus under biaxial loading is $398\pm 5$ GPa, and under uniaxial loading it is $401.5\pm 5$ GPa. Therefore, the elastic modulus and Poisson ratio of NHG at [math] K are found to be $400\pm 5$ GPa and $0.27$, respectively.

In the next step, we performed MD simulations of uniaxial tension at different temperatures. To compare the elastic modulus predicted by MD with that obtained by DFT, we performed the MD modeling of tension at $2$ K. In this case we found an elastic modulus of around $360\pm 5$ GPa for NHG, which corresponds to a $10\%$ difference from our prediction based on DFT. This agreement further validates our choice of parameters for the Tersoff potential and confirms the validity of our MD modeling even though the parameterizations were not trained for NHG. In addition, it should be noted that the Tersoff potentials proposed by Lindsay and Broido were developed to simulate heat conduction along graphene and hexagonal boron-nitride layers and not to exactly reproduce the elastic constants [12, 13]. Furthermore, it is worthy to note that Tersoff potential developed by Lindsay and Broido slightly underestimate (around $3\%$) the elastic modulus of graphene [29]. Fig. 6, present stress-strain relations for NHG at various temperatures from $300$ K to $800$ K. At room temperature, our MD simulations predict an elastic modulus of $335\pm 5$ GPa and tensile strength of $60$ GPa at corresponding failure strain of $0.206$. These values are comparable to single-layer graphitic carbon nitride structures [28]. Our results show overall linear decreasing trends in elastic modulus and tensile strength as the temperature increases. This way, at $800$ K the elastic modulus and tensile strength are found to be around $300$ GPa and $32$ GPa, respectively.

The deformation process of a NHG membrane at $800$ K is depicted in Fig. 7. It is found that at strain levels close to the tensile strength, several defects form along the structure. These defects are mainly initiated by breaking of C-C bonds perpendicular to the loading direction. This initial bond breakage instantly cause breaking of adjacent C-C bonds resulting in the formation of large triangular voids consisting of monoatomic chains of carbon and nitrogen atoms. By increasing the strain, more triangular defects form mostly close to the existing defects throughout the sheet ( Fig. 7(b) ). The tensile strength is then found to be a point in which the coalescence of existing defects occur ( Fig. 7(c) ) leading to the sample rupture (Fig. 7(d) ). At high temperatures, the formation of initial defects are slightly earlier than the tensile strength point. In this case, as the initial defects form, the load bearing ability of the structure is reduced. This crack initiation could be realized from the stress-strain response which is the point that the stress start to present non-linear response. Nevertheless, we found that by decreasing the loading temperature the crack initiation and tensile strength occur at closer strain levels. In this regard, at room temperature the crack initiation happens at the tensile strength point and the stress-strain relation follow almost linear pattern up to the tensile strength.

3.2 Electronic properties

Fig. 8, shows the partial electron charge density distribution of NHG. Electrons with energies ranging from $-1$ eV up to the Fermi energy were included in the estimation of charge density. Our results illustrate higher concentration of electron charge on the more electronegative nitrogen atoms, especially the lone pairs of electrons are clearly visible. It is found that not all C-C bonds are equivalent as seen from the different distributions of electron charges on them. For half of the C-C bonds, the electron charge is concentrated at the center of the bonds. While for the other half, the electron charge is concentrated on both carbon atoms. These two types are regularly distributed depending on the geometry of the neighboring atoms.

NHG is a semiconductor with a finite band gap of $1.96$ eV as determined optically or $1.7$ eV as determined by DFT [9]. Fig. 9 shows the evolution of the electronic density of states and the energy band gap with uniaxial strain as calculated in this work. At completely relaxed state (strain=$0.00$) the band gap was found to be $1.71$ eV, in agreement with the previous DFT study mentioned above. Our first principles calculations show that the band gap decreases monotonically with increasing strain. These results suggest that the band gap of NHG can be manipulated via mechanical strain.

3.3 Point defects and thermal stability

In addition to insights from classical MD simulations, we employed DFT to further investigate the stability of NHG against point defects. A supercell of 72 atoms was used for this calculation. We investigated four point defects which included two cases of monovacancy by removing a C atom (case a) or a N atom (case b), and two cases of divacancies where we removed a pair of C atoms (case c) or one C and one N atom (case d). After introducing the defect, the geometry of each structure was relaxed via energy minimization. Fig. 10 presents the energy-minimized configurations, and shows that the introduction of defects does not change the global structure of NHG. During minimization, new bonds were formed with neighboring atoms in order to compensate for missing atoms. The overall distortions caused by the defects are local by nature, which confirms the robustness and stability of NHG’s structure against defects.

Next, the thermal stability of the defective structures was also tested via AIMD simulations at $300$ K. Each structure was simulated for a total time of $5.0$ ps, and the defective structures remained intact at the end of each simulation. Comparing the potential energies of each structure during AIMD simulations we concluded that divacant structures are less stable than monovacant structures. We also noticed that the monovacant C structure is less stable than the monovacant N structure, and that the divacant CC structure is less stable than the divacant CN structure.

To further elaborate on the role of defects in the thermal stability of NHG we conducted AIMD simulations at very high temperatures, namely $T=3000$ K and $T=4000$ K for pristine and defective NHG. These simulations were also conducted for $5.0$ ps. At the end of each simulation pristine NHG and all four defective structures were intact at $T=3000$ K. This suggests high thermal stability of NHG in spite of local defects. At $T=4000$ K the pristine NHG structure was partially disintegrated at the end of the simulation, and all of the four defective structures were completely disintegrated at this temperature. This shows that point defects can play a role in initiating thermal disintegration of NHG at very high temperatures. Thus, we conclude that the point defects considered in this study only affect NHG locally. Nonetheless, at a very high temperature these defects can play a role in disturbing its structural integrity and thermal stability.

3.4 Thermal conductivity

Thermal conductivity calculations with the NEMD method present strong size dependence effects [30]. In order to investigate this size dependence in NHG we performed simulations with periodic cells of increasing length and a constant width of $23$ nm, as shown by the data points in Fig. 11 for three temperatures. For very small simulation cells $\kappa_{L}$ increases linearly with sample length, which corresponds to the ballistic transport limit. In NHG, for large enough cells the conductivity converges to a finite value, corresponding to the diffusive regime. At intermediate lengths, the system is in a cross-over region from ballistic to diffusive heat transport. In general, the dependence of $\kappa$ on $L$ can be described by [26]

[TABLE]

where $\kappa_{\infty}$ is the size-independent thermal conductivity of the material (i.e. the conductivity in the diffusive regime), and $\Lambda$ is an effective phonon mean free path (MFP) for the material. From the equation we see that $L=\Lambda$ yields $\kappa_{L}=\kappa_{\infty}/2$, such that this effective MFP corresponds to the system length at which the system reaches $50\%$ of its diffusive thermal conductivity [31, 32, 33]. Therefore, we can simultaneously estimate the intrinsic diffusive thermal conductivity of a given material as well as its effective phonon MFP via Eq.5. The solid lines in Fig. 11 represent the best fit to Eq. 5 for each temperature. The values of $\kappa_{\infty}$ and $\Lambda$ for each temperature are summarized in Table 1.

We predict a thermal conductivity of $64.8$ W/m-K for NHG at room temperature, which is two orders of magnitude smaller than the thermal conductivity of graphene at the same temperature [34, 30, 35]. The effective phonon MFP of isotopically pure graphene at $300$ K can be estimated to be $\approx 5$ $\mu$m [35], also two orders of magnitude longer than the calculated $\Lambda=34.0$ nm for NHG. We attribute such difference mostly to the lower density of NHG relative to graphene, as well as the presence of N atoms in the structure which soften the phonon modes of pristine graphene. Nonetheless, it is interesting to notice that the problem of collective excitations with very long MFP does not materialize in NHG [35], even near room temperature. Indeed, from Fig. 11 we observe that $\kappa$ seems to converge at a sample length of $\approx 150$ nm even at room temperature. It is also interesting to compare the effective phonon MFP with the characteristic length scales of this crystal structure, namely the diameter of the holes and the distance between them. At room temperature, the average hole diameter equals $0.62$ nm while the average distance between them is $0.88$ nm. Therefore, at room temperature, both length scales are much shorter than the effective phonon MFP.

From Table 1 we observe a decrease of both $\kappa$ and $\Lambda$ as temperature increases, as qualitatively expected. In Fig. 12 we present the thermal conductivity for samples with a fixed length $L=150$ nm, as a function of simulation temperature, for temperatures above $300$ K. In this temperature range, and in the absence of structural defects and impurities, the thermal conductivity of a bulk crystalline material is expected to be inversely proportional to temperature [36, 37]. In other words, when the only source of thermal resistance is phonon-phonon scattering then $\kappa\sim T^{-1}$. However, in the presence of defects and impurities, $\kappa$ becomes less dependent on $T$, and the exponent assumes a smaller (absolute) value. When phonon-defect scattering dominates over phonon-phonon scattering, the conductivity is expected to be temperature-independent, such that one can write $\kappa\sim T^{0}$ [36, 38]. Fig. 12 shows that our data agrees with an exponent smaller than unit, a behavior we attribute to the dimensionality of NHG, as well as to the presence of holes and two different atomic species in the material. Based on our AIMD simulations, we believe that no structural defects are formed in our MD simulations. Nonetheless, the scattering rate of phonons by lattice holes increases with temperature, leading to higher thermal resistance and lower thermal conductivity. We have recently performed a study of mechanical properties and thermal conductivity of amorphized graphene samples, in which we found that certain types of defects produce a larger decrease in thermal conductivity  [39]. On the same study we observed that, for samples with high defect concentration, thermal conductivity is almost temperature-independent. A full characterization of structural defects and their individual effect on the thermal conductivity of NHG shall be addressed in future work.

4 Conclusions

In conclusion, we performed first principles calculations and extensive atomistic MD simulations to investigate mechanical properties and thermal conductivity of NHG, a novel two-dimensional material, at various temperatures. MD simulations were based on the Tersoff potential with parameters obtained from graphene and hBN parameterizations. The phonon dispersion of NHG presents no instabilities with the chosen empirical potential. Mechanical properties were evaluated via uniaxial tensile strain to obtain elastic constants. First principles DFT calculations predict an elastic modulus of $400\pm 5$ GPa at [math] K, which is $10\%$ larger than the value obtained by atomistic MD simulations at very low temperatures. Results show an overall decreasing trend in elastic modulus and tensile strength as temperature increases. At room temperature, our MD simulations predict an elastic modulus of $335\pm 5$ GPa and a tensile strength of $60$ GPa at a corresponding failure strain of $0.206$. During mechanical failure under strain we observed defect initiation by breakage of C-C bonds mainly perpendicular to the loading direction, formation of large triangular voids consisting of mono-atomic chains and finally void coalescence leading to the final rupture of the sample. We also investigate the thermal conductivity of NHG via NEMD simulations, in particular its length dependence. Fitting simulation results to an empirical expression for $\kappa$ we predict an intrinsic thermal conductivity of $64.8$ W/m-K and a corresponding effective phonon mean free path of $34.0$ nm at $300$ K. Both limiting thermal conductivity and effective phonon mean free path decay with temperature for temperatures above $300$ K. Finally, due to the current lack of experimental studies on the physical properties of NHG, our modeling-based predictions should serve as a guide to future experiments concerning the physical properties of this novel 2D material.

Acknowledgments

B.M. and T.R. greatly acknowledge the financial support by European Research Council for COMBAT project. L.F.C.P. would like to thank Hasan Sahin and Mehmet Yagmurcukardes for providing DFT-optimized unit cell structure for comparison of phonon dispersions. L.F.C.P. acknowledges financial support from Brazilian government agency CAPES for project “Physical properties of nanostructured materials” via its Science Without Borders program.