Thermodynamic properties of pure and doped (B, N) graphene

Abstract

Ab initio density functional perturbation theory (DFPT) has been employed to study the thermodynamical properties of pure and doped graphene sheet and the results have been compared with available theoretical and experimental data. The concentration of B and N has been varied upto 50% of the carbon atoms in graphene. Phonon frequencies are essential ingredients entering into such a calculation, which have been computed by using the dynamical matrix provided by VASP software in combination with phonopy code in the harmonic approximation. This easily provides us with the Helmholtz free energy and leads us to numerical estimates of various thermodynamical properties. The results for specific heat are in good agreement with various theoretical and experimental studies obtained earlier for pure graphene. Interesting new results have been reported for B and N substituted structures. It has been observed that the specific heat decreases with the increase in concentration of doping while the entropy increases. Furthermore, large doping concentrations result in unstable sheets leading to imaginary frequencies in the transverse directions. The instability needs to be compensated by external strains and that has been assumed while carrying out Brilluoin summations. These results will be useful for calculations of thermal conductivity of doped graphene and thus feasibility of using these for device applications. Our preliminary results for thermal expansion have indicated negative thermal expansion behavior of pure graphene at low temperatures which needs investigation on doped graphene as well.

---

1. Introduction

The synthesis of graphene back in 2004,¹ as a two dimensional thermodynamically stable crystal, led to a lot of research interest by virtue of its unique band structure and salient features.² It is basically a planar sheet of sp² bonded carbon atoms arranged in a honeycomb lattice structure. Its various properties like ballistic transport at room temperature, high intrinsic mobility, quantum Hall effect, combined with mechanical strength make it a potential candidate for a lot of applications.^3,4 But its use in electronics is limited by its unique band structure, which leads to the absence of a band gap. So a lot of experimental and theoretical studies have been done to modify its electronic and optical properties. We in our previous work^5–7 have done a detailed study to engineer band gap of graphene and optical properties by substituting the carbon atoms with boron and nitrogen atoms, i.e. forming a hetero structure.

In addition to its unique electronic properties, graphene also has very high thermal conductivity⁸ leading to unique transport properties. In the past few years, various experimental^8–10 reports on heat transport in graphene have shown an ultrahigh thermal conductivity, which is important for applications, such as thermal management in electronics and improvement in thermal conductance of composite materials. Properties of phonons in graphene have recently attracted significant attention from the physics and engineering researchers. Acoustic phonons are the crucial heat carriers in graphene in the vicinity of room temperature. Therefore it is important to study its thermodynamic properties in detail, which are governed by phonons, quanta of the lattice vibrational energy.

As mentioned phonon dispersion relations are essential ingredients of thermal conductivity and for pure graphene these have been extensively studied theoretically along with thermal conductivity in recent times.^11–13 Phonon dispersion relations of graphene have also been measured experimentally using inelastic X-ray scattering¹⁴ as well as high resolution electron energy loss spectroscopy (HREELS)¹⁵ and the theoretical estimates seem to be in reasonable agreement with experimental data. Nika et al.¹³ has shown both theoretically and experimentally that transport properties of phonons are significantly different in a quasi-2-dimensional system such as graphene compared to the basal planes in graphite or three-dimensional bulk crystals. They pointed out that unique nature of two-dimensional phonon transport translates into unusual heat conduction in graphene and related materials. They described different theoretical approaches used for phonon transport in graphene discussing contributions of the in-plane and out of-plane phonon modes, and compared with available experimental thermal conductivity data. H. Zhang et al.¹⁶ have performed molecular dynamics simulations to investigate phonon transport in graphene at 300 K with the Green–Kubo method. They have shown a comparison of the phonon dispersions of graphene among the results from density function theory (DFT), the original reactive empirical bond-order (REBO)¹⁷ and the optimized REBO method.¹⁸ Therefore phonon related properties of pure graphene are reasonably well documented and studied.

It has been recently pointed out that B and N doped graphene has interesting and designable electronic and optical properties^5–7 which can be suitable for device application. The issue which remains unresolved is rapidity of heat dissipation for these doped materials. In order to address this issue it is important to estimate thermal conductivity pattern of these materials. As a first step in this direction we attempt to calculate phonon dispersion and related thermodynamic properties of doped materials.

Some studies have been found in literature, which report the effect of B (N) doping on phonon related properties of graphene. Sherajul et al.¹⁹ have presented an extensive and systematic numerical study of the vibrational properties of B- and N-doped graphene with native vacancies using the forced vibrational method. They have computed the change of the phonon density of states for different concentration of B and N impurities or vacancies without explicitly discussing the phonon dispersion relations. Yanagisawa et al.²⁰ have studied the phonon dispersion in stable and metastable BC₃ honeycomb epitaxial sheets both experimentally and theoretically.

There has been some interesting study reported on strain dependence of some phonon modes of graphene²¹ as well as of silicene.²² Further, Jian et al.²³ have done detailed study (upto 50%) of B and N doped graphene and discussed the stability of such structures on the basis of negative frequencies which occur in these structures. They emphasize that if appropriate strain as well as charge is applied to the structures at high doping concentration, it leads to positive frequencies. But to the best of our knowledge, there is not any systematic study, describing the pattern of changes in phonon dispersion curve because of doping and thermodynamic properties resulting from doped structure like specific heat, entropy and free energy, found in literature so far. So in the present work, we attempt to analyze the effect of boron and nitrogen substitution on stability and phonon dispersion curves of graphene sheet and calculated various thermodynamic properties of pure and doped graphene. While discussing thermodynamic properties, as pointed out by Jian et al., negative frequencies do occur indicating the instability of heavily doped graphene sheet in the transverse direction. Assumption of a strain in the transverse direction can be used to provide stability and remove negative frequencies.

The paper is structured as follows. After summarizing the theoretical formulation and computational details in first two sections, we discuss the results in section three. The results are summarized and concluded in Section 4.

2. Theory and computational details

To calculate the phonon dispersion and thermodynamic properties we have made use of VASP (Vienna ab initio Simulation Package)^24,25 code based on density functional theory (DFT) in interface with phonopy²⁶ code. In density-functional theory, the ground state electronic density and wave functions of a crystal are found by solving self-consistently a set of one-electron equations. These procedure and equations are well described in literature.²⁷ The Kohn–Sham equation is solved self consistently to find the total energy of the system. The VASP program is run to accumulate data of forces at displaced structures which will be used in subsequent sections.

2.1 Dynamical matrix

In order to obtain phonon frequencies we need to obtain the interatomic force constants Φ_αβ,^28–30 i.e.where lκ are positions of any atom w. r. t other atom at l′κ′. α and β are Cartesian indices, r is position vector of atom. In order to obtain the force constants, we further express them in terms of forceswhere F_β(lκ) is the force on the lκ atom in β direction and Δr_α is finite displacement. Since the VASP program gives us the forces, therefore eqn (2) helps us to obtain force constants Φ_αβ and the dynamical matrix.²⁶ Diagonalizing the dynamical matrix gives us phonon frequencies ω_qi, for wave vector q and branch i.

2.2 Thermodynamical properties

In a harmonic crystal, the structure does not depend on temperature. The expressions for heat capacity per unit cell at constant volume and entropy (S)²⁶ are obtained by summing over all the phonon branches.

The expression for free energy term²⁹ used here for calculating free energy as given below

where ϕ(V₀) is the ground state energy at the relaxed volume V₀.

2.3 Computational parameters

The graphene sheet consisting of 32 atoms has been used for simulations and sheets are separated by larger than 9.8 Å along perpendicular direction to avoid interlayer interactions. We adopted the Perdew–Burke–Ernzerhof (PBE)³¹ exchange correlation (XC) functional of generalized gradient approximation (GGA) in our calculations. The plane wave cut-off energy was set to 750 eV for graphene. The Monkhorst–Pack scheme³² is used for sampling Brillouin zone. The structure is fully relaxed with Gamma centered 7 × 7 × 1 k-mesh. The partial occupancies were treated using the tetrahedron methodology with Blöchl corrections.³³ For geometry optimizations, all the internal coordinates were relaxed until the Hellmann–Feynman forces were less than 0.005 Å.

The phonopy software is used in combination with VASP to produce phonon frequencies. To obtain the force constants (FCs) in phonopy, we manually displace the atoms in the center cell by Δx^α (α denotes the three directions), and then calculate the IFCs of the perturbed structure.

3. Results and discussion

3.1 Phonon dispersion

The investigation of lattice vibrational modes is the initial step to study the thermodynamic properties of a material. The simplest approximation is harmonic approximation, where all the oscillators in a solid are harmonic and have fixed frequencies. Using eqn (2) and diagonalizing the dynamical matrix we obtain phonon dispersion curve. The dispersion is basically the relationship between the phonon frequency and the phonon wave vector q. In this section phonon dispersion results for pure graphene have been presented and compared with available theoretical and experimental data. After the successful reproduction of results for pure graphene, the phonon dispersion curves for doped graphene have also been obtained.

3.1.1 Graphene unit cell (2-atoms). Graphene is a hexagonal 2D sheet of carbon atoms arranged in a honeycomb lattice with two atoms in the unit cell as shown in Fig. 1 where a₁ and a₂ are lattice vectors such that

|a₁| = |a₂| = lattice constant = a₀ = 2 × b cos 30° = √3b, (b = bond length)

Fig. 1 Schematic representation of 2D lattice structure of graphene with a₁ and a₂ represent the lattice vectors of unit cell and b shows the bond length at the edge of hexagon.

In a lattice with a basis of two atoms in the primitive cell, there will be four allowed frequencies of wave, two upper branches known as the longitudinal and transverse optical branches abbreviated as LO and TO, and two lower branches called longitudinal and transverse acoustical branches abbreviated as LA and TA respectively. Acoustic phonons have frequencies that become smaller at long wavelengths, and correspond to sound waves in the lattice. We artificially managed a distance of approximately 10 Å between two graphene sheets and displaced the atoms in z-direction also which generates out of plane transverse acoustic (ZA) and optical (ZO) varieties respectively. Thus we get 6 phonon branches in all for a unit cell of graphene. We now discuss in the following subsections phonon dispersion of pure and doped graphene.

3.1.1.1 Pure graphene. In this section, using the theoretical approach described above, we compute the phonon dispersion in single layer graphene which is shown in Fig. 2. The key feature of phonon dispersion in graphene is the existence of ZA mode also termed as flexural mode³⁴ which is the lowest frequency mode and easiest to be excited. The origin of this mode is only due to the surface interactions. The behavior of ZA mode is quadratic in the vicinity of gamma point. There are other important key features of such modes which are derivable from longitudinal modes and exist only for low dimensional systems.³⁵ In fact both ZA and ZO originate from longitudinal motion but ZA is more important because it is the lowest frequency mode.

Fig. 2 Phonon dispersion curve calculated for pure graphene (solid lines) as compared to that in ref. 16 (theory) (dashed lines) and ref. 15 (symbols) (experimental). LA(LO), TA(TO) and ZA(ZO) are longitudinal, transversal and out-of-plane acoustical (optical) branches respectively.

The phonon dispersion curves obtained above are in good agreement with phonon dispersions calculated using VASP by Zhang et al.¹⁶ Zhang et al. have plotted the phonon dispersion curve using VASP and REBO potential for pure graphene. Our results are also similar to³⁶ where the phonon dispersion is calculated by constructing a dynamical matrix using the force constants derived from the second-generation reactive empirical bond order potential by Brenner and co-workers. Finally we also compared our results with experimental data obtained by Mohr et al.¹⁴ They have studied phonon dispersion by Raman 2D-band splitting in graphene and Yanagisawa et al.¹⁵ have used high resolution electron energy loss spectrometer technique for calculating the phonon dispersion in graphene. Result of experimental technique¹⁵ matches with our calculation as shown in Fig. 2.

We also calculated the phonon frequencies for different approximations of exchange correlations. The data is presented along with values provided by Mounet et al.³⁷ (theory) and Yanagisawa et al.¹⁵ (experimental) in Table 1. All the acoustic modes vanish at Γ point.

Table 1 Calculated phonon frequencies for the Γ, M and K points

High symmetry point	Frequency (THz)
	LDA	GGA-PBE	GGA-91	Ref. 37 (GGA-PBE)	Ref. 15 (experimental)
ΓZO	26.90	26.13	26.23	26.41	26.02
ΓLO	45.99	46.45	46.05	46.58	47.36
ΓTO	46.11	46.57	47.02	46.58	47.36
MZA	14.44	14.00	13.97	14.12	13.52
MTA	18.67	18.61	18.52	18.76
MLA	39.33	39.54	39.48	39.81	39.81
MZO	19.47	18.88	19.11	19.03	19.03
MLO	39.73	40.25	40.01	40.17	39.66
MTO	40.78	41.36	41.17	41.67	41.67
KZA	16.44	15.90	16.10	16.03	15.49
KZO	16.73	16.20	16.18	16.03	17.62
KTA	29.49	29.47	29.48	29.88
KLA	35.33	35.64	35.41	36.36	36.36
KLO	35.91	36.23	36.27	36.36	38.61
KTO	38.90	39.59	39.47	38.61	38.61

---

The comparison of phonon frequencies for different exchange correlation shows no significant variation. Therefore this establishes that our procedure is appropriate for computation of phonon frequencies and therefore we extend our computation to phonon dispersion of doped graphene.

3.1.1.2 Doped graphene. We started with a doping of 50% in pure graphene which is obtained by replacing one atom of carbon with boron (nitrogen) in a unit cell of graphene and then repeating the unit cell along xy direction to prepare a 32 atom h-CB/h-CN (hexagonal) sheet. We analyzed the ordered h-CB and h-CN in its neutral state and their phonon dispersion curves, obtained by using similar the technique as for pure graphene, are presented in Fig. 3. Fig. 3(a) shows the phonon dispersion of relaxed structure with equilibrium lattice constant a₀ as 2.698 Å. Fig. 3(b) shows the phonon dispersion of relaxed h-CN structure with equilibrium lattice constant a₀ as 2.397 Å. Phonon dispersion curves for h-CB shows that the lower ZA branch has imaginary frequencies (following conventional notation, we plot imaginary frequency as negative values) indicating that the structure is unstable. The phonon dispersion curve of h-CN shows two branches having imaginary frequency modes. The lower soft branch corresponds to out-of-plane acoustic (ZA) mode, and the upper branch is out-of-plane optic (ZO) mode.

Fig. 3 Phonon dispersion curves for 50% (a) B-doped graphene (b) N-doped graphene.

The imaginary frequency modes indicate that the neutral h-CB and h-CN are dynamically unstable, that is why it is difficult to synthesize these structures. Although we relaxed our structure but the nature of phonons confirms that it is not a stable structure. We therefore conclude that 50% doped structure is unstable under normal conditions.

As discussed earlier Zhou et al.²³ have shown that these structures can be stabilized by applying appropriate strains alongwith doping with opposite charge. They modify the bond length to introduce strains. Therefore in the following section we assume that such doped structures are possible to be stabilized by applying strain alongwith the charges.

3.1.2 Graphene unit cell (8-atoms). In order to vary the concentration of doping, it is necessary to take bigger unit cell. We take unit cell comprising of 8 atoms which allows us to vary the concentration of dopants from 1/8 i.e.12.5% onwards upto 4/8 i.e. 50%. The possible configurations considered by us for the present calculations are shown in Fig. 4. Each structure is initially relaxed and lattice constant is obtained for minimum energy configuration. The electronic properties of such structures have already been investigated by our group.⁵

Fig. 4 Configurations of (a) 12.5% (b) 25% doped, (c) 50% and (d) 37.5% doped (B, N) relaxed structures.

3.1.2.1 Pure graphene. The phonon dispersion of pure graphene for a graphene unit cell of 8 atoms has been shown below in Fig. 5, which gives rise to 3N phonon branches i.e. 24. This dispersion is similar to the one calculated above for a unit cell (2-atoms) of graphene. However there are artificially generated (degenerate) branches which can be correlated with six branches of Fig. 2.

Fig. 5 Phonon dispersion curves of pure graphene for 8-atoms unit cell.

3.1.2.2 Doped graphene. By varying the concentrations as described in Fig. 4 we obtained the graphene structures with different doping concentrations. We optimized these structures for minimum forces but they do not give a correct picture of stability in z-direction, this will be manifested once we discuss the results of phonon frequencies. The lattice constant for an 8-atom unit cell becomes double the original lattice constant. i.e. a′₀ = 2a₀. In case of B-doping, the lattice constant a′₀ keeps on increasing from (2 × 2.46) 4.92 Å for pure graphene to 5.34 Å for 50% doped structure. Inversely the lattice constant a′₀ for N-doped structure decreases from 4.92 Å to 4.79 Å for 50% doped structure. The structure with 12.5% and 25% doping of B and N atoms separately in graphene sheet are stable with no imaginary frequencies and it becomes slightly unstable as we increase the percentage of doped atoms to 37.5% and highly unstable for 50% doped structure as shown in Fig. 3 earlier. The variation of lattice constant with doping is given in Table 2.

Table 2 Variation of lattice constant with doping concentration

Doping conc.		12.5%	25%	37.5%	50%
Lattice constant (a′0)	B-doped	5.04	5.17	5.25	5.34
	N-doped	4.90	4.86	4.82	4.79

---

The study of phonon dispersion of B doped graphene for these doping concentrations is shown in Fig. 6.

Fig. 6 Phonon dispersion curve of B doped graphene (a) 12.5% (b) 25% (c) 37.5% and (d) 50% respectively.

Similarly, the phonon dispersion curves for N doped graphene for various doping concentrations have been obtained and shown in Fig. 7.

Fig. 7 Phonon dispersion curve of (a) 12.5% (b) 25% (c) 37.5% and (d) 50% N-doped graphene respectively.

We observe that for the case of 12.5% doping (Fig. 6(a) (B-doping) and in Fig. 7(a) (N-doping), all the frequencies are positive but they no longer are degenerate branches as in pure graphene case in Fig. 5. Thus we have separate branches with even minimal of doping.

When we move to 25% doping in Fig. 6(b) (B-doping) and in Fig. 7(b) (N-doping), the phonon dispersion shows a softening of acoustic and optical branches. The frequency of lowest ZA branch is decreased in both the cases.

For the higher doping concentrations i.e. 37.5% and 50% as shown in Fig. 6(c and d) (B-doping) and in Fig. 7(c and d) (N-doping), we get imaginary frequencies which are more in number in 50% doped case. There is a noticeable feature that the maximum phonon frequency range largely reduces in case of B doped structure from 48 THz to 35 THz. Although there is decrease in maximum frequency in N doped structures also but there is no prominent change. Therefore as far as for these concentrations are concerned, our structure is unstable. Further, the modes which became imaginary are ZA and ZO branches only.

3.2 Phonon DOS

The DOS curves are plotted in Fig. 8 for doping concentrations from 12.5% to 50% as compared to pure graphene. In obtaining DOS and various thermodynamic properties, we assumed that by applying proper value of strain and charging as reported by the work done by Zhou et al.,²³ we get the structures with nearly zero frequencies of ZA and ZO branches in case of 37.5% and 50% doped structures which have imaginary frequencies.

Fig. 8 Phonon density of states for (a) B-doped and (b) N-doped structures with varying doping concentration along with pure graphene.

The comparison of phonon density of states clearly indicates the difference from pure graphene case. As the B and N doping is increased, the peaks corresponding to pure graphene are shifted towards lower frequency and the intensity of peaks also increased. The DOS also indicates formation of some new peaks in addition to pure one as reported earlier in ref. 19.

3.3 Lattice thermodynamical properties

Thereafter we computed thermodynamical properties of pure graphene and compared with doped graphene. The specific heat of pure graphene as calculated by us has been compared with experimental data³⁹ as well as with earlier calculations^37,38 and presented in Fig. 9. It is observed that there is good agreement with the experimental data. The specific heat, entropy and free energy are calculated using eqn (3)–(5) and have been presented in Fig. 10 and 11. For the case of high doped graphene where the imaginary frequencies are encountered, we take them to be zero to calculate thermodynamical properties. This assumption is justifiable in view of an earlier proposal²³ of disappearance of negative frequencies under appropriate strain and charging. We have assumed that these heavily doped graphene are under appropriate strain and charging which turns the negative frequencies of these structures to nearly zero frequencies.

Fig. 9 Specific heat of pure graphene at various temperatures compared with experimental³⁹ and other theoretical calculations.^37,38

Fig. 10 Variation of specific heat, entropy and free energy with temperature for B doped graphene with increase in doping concentration.

Fig. 11 Variation of specific heat, entropy and free energy with temperature for N doped graphene with increase in doping concentration.

It is easily observed from the graph as shown in Fig. 10(a) and 11(a) that specific heat for pure graphene as well as doped graphene reaches a constant value 21 J K⁻¹ mol⁻¹ which is close to its Dulong–Petit value at high temperature (1000 K),^38,39 however near room temperature (300 K) it is around 7 J K⁻¹ mol⁻¹ (0.6 J g⁻¹ K⁻¹) for case of pure graphene. We notice that entropy and free energy show systematic pattern in their dependence in case of B doped graphene. The behavior of all the thermodynamic properties for the case of N doped graphene show some deviations at certain temperatures as shown in Fig. 11(a–c). Due to soft modes in the transverse direction, there is a tendency for specific heat to increase at low temperature and decrease at high temperatures for highly doped configurations. Unfortunately, there is no experimental or other theoretical data on the lattice dynamical properties of graphene for the sake of comparison.

The specific heat, entropy and free energy values of 12.5% B and N doped graphene are close to that for pure graphene. The specific heat of doped graphene increases with doping at low temperatures for all doping concentrations. The increase in specific heat can be accounted to the higher phonon density of states at low frequencies as shown in Fig. 10. The behavior of specific heat for doped graphene decreases to a constant value above cut off (Debye) temperature. But the specific heat for 50% B and N doped graphene structures is strikingly becoming constant at a much lower value at high temperatures. 50% doping structure shows peculiar behavior in N doped structure. This may be due to highly suppressed transverse modes.

The entropy is increasing with increase in doping concentration for both B and N doped structures. As the B doping concentration is increased to 50%, there is large increase in entropy as compared with pure graphene case. But there is a significant decrease in entropy of 50% N doped graphene at higher temperatures approaching the pure graphene case, indicating a favorable condition for formation of the h-CN sheet at higher temperatures, whereas there is no such decrease in value of entropy in case of B-doped graphene. This is probably due to similar size of C and N, which favors the formation of stabilized h-CN sheet at high temperatures which is indicated by the increase in entropy.

The variation of free energy with temperature for various doping concentrations (Fig. 10(c) and 11(c)) shows a constant decrease with increase in B(N) doping concentrations. Unlike, in case of N doping, the behavior is similar for all doping concentrations except 50%. The 50% doped curve for free energy is lower at low temperature values and cuts the 37.5% curve at around 600 K and becomes higher at high temperatures. Thus free energy curve of 50% N doped structure does not decrease as fast as it should happen. Overall we observe a systematic behavior (increase or decrease) for all the thermodynamic quantities with increasing the doping concentration starting from (12.5%). However, this trend changes unusually at 50% doping concentration, which is more prominent in case of N-doped structures.

We also performed preliminary calculations for Grüneisen parameter defined as

The Grüneisen parameters are a measure of anharmonicity of the potential and are useful parameters for determining thermal expansion and implicit anharmonic properties.³⁰

As shown in Fig. 12 where we present an estimate of Grüneisen parameters of various modes and thus calculated thermal expansion at various temperatures compared with available theoretical and experimental data. Although ZA mode is likely to result in negative value of Grüneisen parameter,³⁵ but such high negative value of Grüneisen parameter for ZA mode is unusual and requires further investigation. Although our results reproduce the results for thermal expansion qualitatively, but there is significant variation in experimental and other theoretical data. There is negative thermal expansion reported in the low temperature region but there are significant variations. The discrepancies between our study and other theoretical^37,41 and experimental⁴⁰ studies also accordingly needs to be fully understood. Once this work is investigated, it will be extended to the doped graphene.

Fig. 12 Grüneisen parameters of various modes in pure graphene at various q values (a) and resulting thermal expansion at various temperatures (b) including other results.

4. Summary and conclusions

We have presented a detailed first-principles study of pure and B and N doped graphene at the GGA-PBE level under harmonic approximation to derive the finite temperature behavior of several thermodynamic quantities. All our results for pure graphene sheet are in very good agreement with previous theoretical and experimental data. The 2-D graphene structure which was optimized under static conditions does not necessarily result in stable 2-D configuration, especially on large doping concentrations. This is evidenced by the results of transverse phonon frequencies which begin to appear as imaginary above critical concentrations of B as well as N. In order to retain these structures as stable which return positive frequencies in the transverse directions as well, we follow Zhou et al.²³ and assume a strained and charged 2-D lattice compensates for this instability. Therefore the role of imaginary frequencies can be overlooked as has been done here. Various thermodynamic properties of B and N doped graphene at different concentrations have been compared with the pure graphene. In general, the heat capacity and free energy of the structures decreases and entropy increases with increasing doping concentration as compared to pure graphene while the entropy (heat capacity, free energy) of the 50% doped structure is dramatically low (high) which is possibly due to the higher symmetry and hence stability of this structure. The heat capacity of doped graphene is higher at low temperatures and decreases to a constant value at high temperature limit as compared to pure graphene. The higher specific heat of doped graphene can be explained on the basis of higher value of density of states at low frequencies. The increasing value of entropy fairly explains the increasing disorder and decreasing stability with doping. The decrease in entropy and increase in heat capacity of the 50% N doped structure is due to symmetry and similar sizes of carbon and nitrogen. Our results suffer from not using appropriate strained lattices under high doping which stabilizes the 2-D structure in the transverse direction and we intend to examine this.

Acknowledgements

We express our gratitude to VASP and phonopy team for providing the code, the HPC facilities at IUAC (New Delhi) and the departmental computing facilities at Department of Physics, Panjab University, Chandigarh.

References

1. K. S. Novoselov, A. K. Geim, S. V. Morozov and D. Jiang, Science, 2004, 306, 666.
2. A. K. Geim and K. S. Novoselov, Nat. Mater., 2007, 6, 183.
3. A. K. Geim, Science, 2009, 324, 1530–1534.
4. P. Avouris, Z. Chen and V. Perebeinos, Nat. Nanotechnol., 2007, 2, 605.
5. P. Rani and V. K. Jindal, RSC Adv., 2013, 3, 802.
6. P. Rani and V. K. Jindal, Appl. Nanosci., 2014, 4, 989.
7. P. Rani, G. S. Dubey and V. K. Jindal, Phys. E, 2014, 62, 28.
8. S. Ghosh, I. Calizo, D. Teweldebrhan, E. P. Pokatilov, D. L. Nika, A. A. Balandin, W. Bao, F. Miao and C. N. Lau, Appl. Phys. Lett., 2008, 92, 151911.
9. W. Cai, A. L. Moore, Y. Zhu, X. Li, S. Chen, L. Shi and R. S. Ruoff, Nano Lett., 2010, 10, 1645.
10. J. H. Seol, I. Jo, A. L. Moore, L. Lindsay, Z. H. Aitken, M. T. Pettes, X. Li, Z. Yao, R. Huang, D. Broido, N. Mingo, R. S. Ruoff and L. Shi, Science, 2010, 328, 213.
11. D. L. Nika, E. P. Pokatilov, A. S. Askerov and A. A. Balandin, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 155413.
12. L. Lindsay, D. A. Broido and N. Mingo, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 115427.
13. D. L. Nika, L. Denis and A. A. Balandin, J. Phys.: Condens. Matter, 2012, 24, 233203.
14. M. Mohr, J. Maultzsch, E. Dobardžić, S. Reich, I. Milošević, M. Damnjanović, A. Bosak, M. Krisch and C. Thomsen, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 76, 035439.
15. H. Yanagisawa, T. Tanaka, Y. Ishida, M. Matsue, E. Rokuta, S. Otani and C. Oshima, Surf. Interface Anal., 2005, 37, 133.
16. H. Zhang, G. Lee and K. Cho, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 115460.
17. D. W. Brenner, O. A. Shenderova, J. A. Harrison, S. J. Stuart, B. Ni and S. B. Sinnott, J. Phys.: Condens. Matter, 2002, 14, 783.
18. L. Lindsay and D. A. Broido, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 205441.
19. M. S. Islam, K. Ushida, S. Tanaka, T. Makino and A. Hashimoto, Comput. Mater. Sci., 2014, 94, 35–43.
20. H. Yanagisawa, T. Tanaka, Y. Ishida, E. Rokuta, S. Otani and C. Oshima, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 045412.
21. Y. C. Cheng, Z. Y. Zhu, G. S. Huang and U. Schwingenschl, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 115449.
22. T. P. Kaloni, Y. C. Cheng and U. Schwingenschl, J. Appl. Phys., 2013, 113, 104305.
23. J. Zhou, Q. Sun, Q. Wang and P. Jena, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 92, 064505.
24. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169.
25. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758.
26. A. Togo, F. Oba and I. Tanaka, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 78, 134106.
27. W. Kohn and L. J. Sham, Phys. Rev. A, 1965, 140, 1133.
28. A. A. Maradudin, P. A. Flinn and R. A. Coldwell-Horsfall, Ann. Phys., 1961, 15, 103.
29. V. K. Jindal and J. Kalus, J. Phys. C: Solid State Phys., 1983, 16, 3061–3080.
30. R. Bhandari and V. K. Jindal, J. Phys.: Condens. Matter, 1991, 3, 899–907.
31. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
32. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Condens. Matter Mater. Phys., 1976, 13, 5188.
33. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953.
34. J.-W. Jiang, B.-S. Wang, J.-S. Wang and H. S. Park, J. Phys.: Condens. Matter, 2015, 27, 083001.
35. V. K. Jindal, 2015, cond-mat.mes-hall: arXiv:1510.02921.
36. D. Kumar, V. Verma, H. S. Bhatti and K. Dharamvir, Pramana–J. Phys., 2013, 81, 1021–1035.
37. N. Mounet and N. Marzari, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 71, 205214.
38. E. Pop, V. Varshney and A. K. Roy, MRS Bull., 2012, 37, 1271.
39. D. Gray, A. McCaughan and B. Mookerji, Physics for Solid-State Application, 2009, 6, 730.
40. D. Yoon, Y.-W. Son and H. Cheong, Nano Lett., 2011, 11, 3227.
41. C. Sevik, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89, 035422.

---