Melting temperature of graphene

PDF hosted at the Radboud Repository of the Radboud University Nijmegen The following full text is a preprint version which may differ from the publisher's version. For additional information about this publication click this link. http://hdl.handle.net/2066/140106 Please be advised that this information was generated on 2020-05-06 and may be subject to change. Melting temperature of graphene J. H. Los, K.V. Zakharchenko, M. I. Katsnelson, Annalisa Fasolino* arXiv:1503.03811v1 [cond-mat.mes-hall] 12 Mar 2015 Radboud University Nijmegen/Institute for Molecules and Materials, Heyendaalseweg 135, 6525AJ Nijmegen, The Netherlands and * [email protected] (Dated: March 13, 2015) We present an approach to the melting of graphene based on nucleation theory for a first order phase transition from the 2D solid to the 3D liquid via an intermediate quasi-2D liquid. The applicability of nucleation theory, supported by the results of systematic atomistic Monte Carlo simulations, provides an intrinsic definition of the melting temperature of graphene, Tm , and allows us to determine it. We find Tm ≃ 4510 K, about 250 K higher than that of graphite using the same interatomic interaction model. The found melting temperature is shown to be in good agreement with the asymptotic results of melting simulations for finite disks and ribbons of graphene. Our results strongly suggest that graphene is the most refractory of all known materials. I. INTRODUCTION Surprisingly, understanding of melting is still an open problem in condensed matter physics. For instance, after more than hundred years, there is still no reliable theoretical justification of the Lindemann criterion [1]. Nevertheless, the melting temperature Tm is well-defined as the temperature at which the Gibbs free energy curves of the liquid and solid phase intersect. Graphene, however, is a quasi-2D crystalline membrane [2], which melts into a 3D liquid phase [3]. The different dimensionality of these two phases makes it impossible to determine their free energy difference with existing methods, so that it is a priori not even clear how to define Tm . Here we show that a 2D nucleation theory offers a way to overcome this difficulty and give a reliable quantitative value of Tm . While for 3D systems, there is no consensus about how to describe melting and premelting anomalies, for a 2D crystal there is a commonly accepted microscopic scenario for melting [4], the KTHNY theory [5–9]. In the KTHNY theory, melting occurs via unbinding of topological defects, like vortices in superconductors [10]. This scenario seems not to apply to graphene. Apart from the embedding in 3D space, which affects all its structural and thermodynamic properties [2, 11], the nature of defects is completely different from those considered in the KTHNY theory. In graphene, there are neither observations nor predictions of stable single pentagons or heptagon defects (disclinations). Dislocations (pentagon-heptagon pairs) have been observed mainly at grain boundaries [12]. The 4.6 eV formation energy [14] of a Stone-Wales (SW) defect made of two adjacent pentagon-heptagons (57) pairs is much lower than the formation energy of two well separated 57 pairs as occurring in small angle grain boundaries [13]. In the Appendix we show how the formation energy of two 57 pairs rapidly increases with separation. The relative stability of SW defects make them play a crucial role in the premelting behavior of graphene [3] and prevents them to split into two 57 defects and further, and by this to follow the KTHNY scenario. In ref.3 we have shown that, within the same model of interactions that we use here, spontaneous melting of ∗ graphene occurs at a temperature that we call Tm of 4900 K, giving an upper limit for the true melting temperature Tm of graphene. Here we show that the melting of bulk graphene can be described by nucleation theory, allowing an unambigous definition of the melting temperature, Tm . We find Tm ≃ 4510 K, about 250 K higher than that of graphite[16] and sofar the highest of all materials. The paper is organized as follows. In section II we review previous results [3] on the spontaneous melting of bulk graphene. In section III we introduce an approach based on classical and kinetic nucleation theory to the melting in 2D. This approach is applied in section IV to bulk graphene. The melting of finite disks and ribbons is presented in sections V and VI respectively. Summary and conclusions are given in section VII. II. SPONTANEOUS MELTING In a previous work [3], we have studied the spontaneous melting of graphene by means of Monte Carlo (MC) simulations based on the reactive bond order potential LCBOPII [14]. A suitable Lindemann type order parameter for γ = (1/a2 ) < Pgraphene was defined as √ 12 (ri − (1/12) j rj )2 >, where a = 1/ πρ is the atomic radius, with ρ the 2D particle density, ri is the position of the i-th atom and the sum over j runs over the 12 atoms closest to atom i. For graphene, melting starts at γ12 ≃ 0.1, close to the value found for a strictly 2D, triangular lattice [15] . In Fig. 1a we show a snapshot from a MC simulation in which a bulk graphene system consisting of N =1008 atoms with periodic boundary conditions (PBCs) is heated at a constant rate. Fig. 1b and c show the potential energy E and γ12 as a function of the applied, increasing temperature found in 8 independent MC simulations. Melting is signalled by a jump in E and γ12 , but with large variations in the observed melting tem∗ perature, Tm , at which this jump occurs for different, independent simulations. This jump as well as the coexistence of solid and liquid parts in Fig. 1a are typical 2 a -4.5 c b -5.0 1.6 -5.5 1.2 -6.0 0.8 -6.5 0.4 E (eV) -7.0 4600 4800 5000 5200 T (K) 4800 5000 5200 γ12 0 T (K) FIG. 1. a) A snapshot of simulated bulk graphene at the onset of melting. The system contains 1008 atoms with PBCs. Note the coexistence of solid and liquid parts. b) Energy E and c) and order parameter γ12 as a function of temperature T during a linear temperature ramp for 8 independent simulations. rations. The importance of assuming the existence of a quasi-2D liquid phase with its own thermodynamic properties is that it allows us to formulate a nucleation theory for the melting of a 2D solid embedded in 3D space, the case of graphene. If the melting can only take place via this intermediate, quasi-2D liquid phase, then this intermediate phase is decisive for the location of the melting temperature of graphene. For the analysis of melting in terms of nucleation theory we performed two series of simulations: simulations at constant temperature and simulations applying a linear temperature ramp. To assess the accuracy of our approach, we will compare classical nucleation theory (CNT) and kinetic nucleation theory (KNT), as described hereafter. We use LCBOPII [14] for the interatomic interactions, as in Ref. [3]. We point out that the melting temperature of graphite according to LCBOPII has been accurately determined to be 4250 K [16]. For the application of CNT to melting in 2D, we write the work W (r) done to form a quasi-2D liquid nucleus of radius r as: W (r) = πr2 ρ∆µ + 2πrγsl features of first order phase transitions. Obviously, the ∗ lowest Tm = 4900 K found in Ref. 3 constitutes only an upper limit for the true melting temperature Tm of graphene. III. NUCLEATION THEORY APPROACH TO MELTING IN 2D Here we propose a quantitative approach to determine Tm by considering the melting as a process initiated and dominated by the nucleation of liquid nuclei in the solid sheet, a typical scenario for a first order phase transition. The possible applicability of nucleation theory is furthermore suggested by our observation that the melting process occurs in two steps. First the graphene sheet transforms into a sort of quasi-2D liquid phase consisting of entangled and interconnected chains remaining roughly within the quasi-2D plane of the previously solid sheet. Eventually, the chains start to disentangle and extend to 3D space, with a diverging simulation box in NPT simulation at zero pressure applied in this work. In Fig. 1 this ”graphene → quasi-2D liquid → 3D liquid” scenario is reflected in the energy curve: first, at melting, the energy steeply increases by about 1.05 eV after which it remains roughly constant for a while before it starts to increase further. This two steps scenario via an intermediate quasi-2D liquid phase would imply that the melting temperature of graphene is in fact the temperature at which the 2D solid and the quasi-2D liquid phases are in equilibrium. Looking at the snapshot in Fig. 1a, it can also be argued that before complete melting, when liquid nuclei start to form in the solid sheet, the liquid phase remains connected with the solid sheet at the edges of the nucleus, constraining it to stay within quasi-2D configu- (1) where ∆µ = µl − µs is the chemical potential difference per particle between the quasi-2D liquid phase and the solid phase, and γsl is the solid-liquid interface free energy. Note that πr2 ρ is the number of liquid particles inside the nucleus, which should indeed grow as r2 due to the 2D geometry of the solid. For T > Tm , ∆µ < 0 and W (r) has a maximum at the critical nucleus size with radius rc = −γsl /(ρ∆µ) equal to Wc ≡ 2 W (rc ) = −πγsl /(ρ∆µ). Using the thermodynamic relation ∆µ = ∆h − T ∆s ≈ (∆h/Tm ) (Tm − T ), with ∆h and ∆s = ∆h/Tm the melting enthalpy (latent heat) and melting entropy per particle respectively, the nucleation probability is derived to be:   β αTm CN T Pnucl ≡ K exp (−βWc ) = K exp − (2) (T − Tm ) where K is a kinetic prefactor, β = 1/(kB T ) and we 2 defined α ≡ πγsl /(ρ∆h). If the melting is dominated by nucleation, then the average time t∗ required for melting to occur is given by the solution of the equation: Nnucl (t∗ ) = Z t∗ CN T Pnucl (T (t))dt = 1 (3) 0 where Nnucl (t) is the number of (super)critical nuclei at time t. At constant temperature T > Tm , eq. 3 simplifies to t∗ = 1/Pnucl (T ). Instead, if the system is slowly heated at a heating rate η, such that T = Tm + ηt, we have to solve eq. 3 numerically for t∗ . The temperature ∗ at which melting is expected, Tm , is then readily obtained ∗ ∗ from Tm = Tm + ηt . In principle, CNT does not specify neither the prefactor K nor its temperature dependence. Therefore, we 3 where fc is the attachment rate for a particle to join a nucleus of the critical size, Z is the so-called Zeldovich factor and α as defined before. The attachment rate can be expressed as: fc = νref T √ g nc exp (−βEa ) Tref (5) where νref ≃ kB Tref /h (with h the Planck constant) is an attempt frequency at some reference temperature √ Tref , which we choose to be Tref = 4500 K, g nc accounts for the number of perimeter sites around a critical 2 nucleus with nc = πrc2 ρ = πγsl /(ρ∆µ2 ) and g a geometrical factor, while Ea is an activation energy barrier for the attachment process. For a circular nucleus, the num√ ber of perimeter sites is given by 2πrc dρ = (2d/a) nc where d is the width of one atomic layer (or shell), whence g = 2d/a. Following Ref. [17], the Zeldovich factor Z for our 2D case is found to be:  1/2 β |∆µ|3 Z= . (6) 4π α∆h Comparing the parameters in KNT to those in CNT, and noting that K̃ ≡ νref g replaces K in CNT, there are two additional parameters: ∆h and Ea . However, ∆h can directly be obtained from the simulations as the melting energy (at zero pressure): ∆h ≃ 1.05 eV (see Fig. 1b). The activation energy barrier Ea is typically a few eV, implying that the temperature dependence introduced by exp (−βEa ) is rather weak within the temperature range of interest and has only a very minor effect (as we verified). We have taken Ea = 2.6 eV fixed, equal to the activation energy for diffusion in a mixed sp-sp2 carbon liquid, as found with LCBOPII [16]. LCBOPII include short range (covalent), long range (van der Waals) and middle range (MR) interactions [14]. 5600 5500 * Tm (K) take K constant in our simulation analysis based on CNT. Nucleation is then determined by three parameters: α, K and Tm . Expressions for the temperature dependence of the kinetic prefactor belong to the domain of kinetic nucleation theory (KNT), which allows us to perform a more accurate analysis. The starting point of KNT [17] is quite different from that of CNT. KNT assumes a (liquid) cluster size distribution (CSD) governed by a coupled set of master equations with appropriate boundary conditions. In CNT, the CSD is proportional to exp (−βW (r)) so that it has a minimum for clusters of the critical nucleus radius rc , where W is maximal, while it increases exponentially beyond the critical size. The latter prediction is clearly unphysical. In KNT, this unphysical behavior is avoided by requiring the CSD to satisfy the boundary condition XN = 0, where XN is the number of clusters that incorporate all N atoms of the system. Eventually, KNT leads to a stationary state nucleation rate given by:   β αTm KN T (4) Pnucl = fc Z exp − (T − Tm ) 5400 MC 5300 5200 MD 1 2 3 η (K/ps) 4 5 FIG. 2. Comparison of the average, observed melting tem∗ perature, Tm , of bulk graphene (1008 atoms with PBCs) as a function of the heating rate obtained from MD (blue diamonds) and MC (red triangles) simulations after calibration, using LCBOPII without the middle range interactions. Error bars are based on 48 independent simulations. Lines are guides to the eye. The latter contribution is a correction that improves the reactive properties and, by construction, does not affect the equilibrium carbon structures at low T . Since many simulations were required to obtain sufficient statistics, and the use of the full potential including the MR part is computationally demanding for Molecular Dynamics (MD) simulations, we have performed MC simulations. Following Ref. 18, we then assume that time is proportional to the number of MC cycles tMC (one cycle being N trial moves) so that t = CtMC where C is a time calibration factor assumed to be constant within a limited temperature range provided that the acceptance rate is kept constant. To assess the validity and accuracy of this approximation, we also performed Molecular Dynamics simulations, but without the MR part of the LCBOPII for the sake of feasibility, and compared the results to those from MC simulations using the same potential (i.e. LCBOPII without the MR part). The MC and MD simulations were performed applying a linear temperature ramp. In the MD simulations the temperature was raised by using the Berendson thermostat and scaling the velocities at every time step, while in the MC simulations the temperature was raised every 5 MC cycles. We considered a series of heating rates, and for each of them the ∗ average Tm was determined from 48 independent simu∗ lations. For each simulation, Tm was determined by the intersection of the γ12 curve with the horizontal line at ∗ as a γ12 = 0.2 (see Fig. 1a). Plotting the average Tm function of η, a very good agreement between the MD and MC results was obtained by applying a time calibration factor equal to C = 5.56 10−4 ps/cycle, as is shown in Fig. 2. We note that for the determination of Tm using nucleation theory, it is sufficient that t = CtMC holds and we do not need C. 4 a 25 b 5100 10 * Tm (K) 7 0 d * c tMC (10 cycles) 5 4900 20 best fit: Tm=4510 K 5100 15 5 4900 0.4 0.8 -3 1.2 1.6 ηMC (10 K/MC cycle) 2.0 4800 4900 5000 T (K) 0 5100 ∗ FIG. 3. Symbols: a) and c) average melting temperature Tm as a function of the MC heating rate ηM C and b) and d) MC ’time’ t∗M C required for melting at isothermal conditions as a function of temperature. The error bars are based on at least 48 independent simulations, but those for the isothermal simulations at 4750, 4800 and 4850 K are based on 96 simulations. The solid (red) curves represent simultaneous best fits of both sets of simulation data based on predictions by CNT (a and b) and KNT (c and d). IV. MELTING OF BULK GRAPHENE We apply the CNT and KNT described in the previous section to evaluate Tm of a bulk graphene system consisting of N =1008 atoms with PBCs as in Fig. 1. In Fig. 3 we present the results of two sets of systematic MC simulations, including simulations heating the system at a rate ηMC and isothermal simulations. In view of the large variations in the respective, observed melting ∗ temperature Tm and time required for melting, t∗MC , at least 48 independent simulations were performed in all cases, in order to obtain good statistical averages. Best fits to all 13 data points were then determined by using the described CNT and KNT approaches, with fitting parameters K (or K̃ = νref g in the case of KNT), α and Tm . In both cases good agreement was obtained, with estimates for Tm equal to 4450 and 4510 K according to CNT and KNT respectively. This relatively small difference shows that the exponential factor dominates the nucleation rate and that the role of the temperature dependence of the prefactor in KNT is limited. Since KNT is more accurate, we conclude that Tm ≃ 4510 K. From α we find γsl = 0.134 eV/Å . MELTING OF FINITE GRAPHENE DISKS Searching for verfication of the result from KNT, we also performed simulations of the melting of a series of roughly circular graphene disks of increasing size. Typically, their melting temperature, Tm,d , increases as a function of cluster size, as was already demonstrated for small graphene flakes (N < 55) in Ref. 19. According to the Pavlov model [20] specified below, Tm,d deviates from c -5.6 d 1 2 -6 3 4500 4 4400 5 T(K) 10 5000 V. b 15 5000 5200 a 20 best fit: Tm=4450 K E (eV) 5200 4300 -6.4 Tm,d 3800 4000 4200 4400 4600 T(K) 0 0.01 0.02 0.03 -1 0.04 4200 1/Rs(Å ) FIG. 4. a) Snapshot at 4300 K, illustrating edge melting of a graphene disk containing 1274 atoms heated at a rate of 2.5 10−4 K/cycle. b) Snapshot illustrating the recrystallization after cooling the system in a) from 4300 to 3800 K. c) Melting curves for different graphene disk sizes: 1) 714, 2) 1274, 3) 2470, 4) 3820 and 5) 8020 atoms. For the latter cluster, melting curves from 4 independent simulations are drawn but they are almost overlapping. The crossing dashed lines illustrate the method used to define the melting temperature Tm,d of a disk. d) Best fit (red line) according to eq. 7 of the average simulated disk melting temperatures (symbols). The error bars are based on 4 independent simulations. Tm by a finite size correction that vanishes with the inverse radius of the disk, allowing for the determination of Tm by extrapolation. It should be noticed, however, that the analysis of the disk simulations is arguable because Tm,d is not rigorously defined and the Pavlov model does not consider edge melting. For a disk, nucleation is not required and the melting starts at the edge, as illustrated in Fig. 4a. According to CNT there is no barrier for nucleation at a 1D edge, which is rough at any T . Therefore there is little variation in Tm,d from different simulations, as is illustrated by the almost overlapping melting curves from 4 independent simulations for the cluster with 8020 atoms in Fig 4c. Interestingly, the melting process is reversible to some extent, as long as there are still solid parts in the system. In Fig. 4b we show the recrystallization that has occurred after the system in Fig. 4a was slowly cooled from 4300 to 3800 K. An appropriate approach to analyse our simulation data for disks would be a 2D version of Pavlov’s model [20]. Writing the free energy of the solid (s) and quasi-2D liquid (l) disk as Gd,P = NP µP + 2πRP γP v for P = s, l, with NP the number of particles in the solid/liquid disk, RP their radii and γP v the respective edge free energies of their interfaces with the vacuum, and imposing equality of the chemical potentials through dGd,s /dNs |Ns =N = dGd,l /dNl |Nl =N leads to the following correction for the 5 melting temperature of a finite disk:   γsv − qρ γlv Tm,d = Tm 1 − ρs ∆hRs (7) where qρ ≡ (ρs /ρl )1/2 accounts for the 2D density difference between solid and liquid. The results of the disk simulations are summarized in Fig. 4. A clear trend in Tm,d as a function of 1/Rs is observed and a best fit based on eq. 7 yields Tm = 4503 K, a value very close to that derived from KNT, and γsv − qρ γlv = 0.58 eV/Å . With an (average) value of γsv ≃ 1.05 eV/Å for a graphene ribbon from Ref. [21] and taking qρ ≃ 1 we find γlv ≃ 0.47 eV/Å . These edge energies, and also γsl = 0.134 eV/Å from KNT reported above, are physically sound, with γsv > γlv > γsl . VI. MELTING OF 2D GRAPHENE RIBBONS For further verification of our result for Tm , we have also performed constant temperature melting simulations of a graphene ribbon. We used a ribbon of a size Lx × Ly = 51.12 × 103.23 Å2 containing 2016 atoms with PBC only in the x-direction. For this geometry, as for the disks, the melting proceeds from the free edges, as shown in Fig. 5a, and there is no nucleation barrier. In this a b 4800 K 4700 K 4600 K 0.75 0.5 4400 K 0.4 0.25 0.3 0 ^ 4300 K 0.2 0.3 0.4 tMC (107 cycles) 0.1 0.5 exp(βEa) dNl/dt K ^ Nl 4500 0.2 0.1 4200 4400 4600 T (K) 4800 0 FIG. 5. a) Snapshot of a melting ribbon consisting of 2016 atom at T = 4500 K. b) The simulated melting speed (symbols), exp (βEa )dN̂l /dt, as a function of the temperature T for a graphene ribbon. The inset shows the fraction of liquid particles N̂l as a function of tM C for various simulation temperatures as indicated. For the completely liquid phase N̂l ≃ 0.5 as many particle still have a γ12 < 0.2 due to a rather isotropic neighbor surrounding. case of a straight edge, the melting speed, dNl /dt with Nl the number of liquid particles, which we define as the number of particles with γ12 > 0.2, is expected to behave as K̂(T )Psl (exp (−β∆µ) − 1), where Psl ≃ 2Lx d/(πa2 ) is the number of sites at the solid-liquid interface and K̂(T ) = (νref T /Tref ) exp (−βEa ). Since |β∆µ| << 1, the melting speed can be well approximated as: dN̂l νref ∆h exp (−βEa ) 2Ly dρ = (T − Tm ) dt kB Tref Tm N (8) where we have defined N̂l ≡ Nl /N as the fraction of liquid particles. In Fig. 5b we have plotted exp (βEa )dN̂l /dt as a function of T obtained from simulations (symbols). The intersection of a best linear fit according to Eq. 8 for T > 4500 K with the horizontal axis yields Tm = 4505 K, very close to the previously found values. It has to be noticed, however, that the results for T = 4400, 4450 and 4500 do not follow the linear law, as the system actually melts for these temperatures. Although finite size effects might play a role here, maybe the most important reason is that chain-like molecular units evaporate from the liquid edges, as illustrated in Fig. 5a. Since equilibrium in the present case is characterized by a liquid edge of fixed average width, evaporation must lead to a thinner and thinner solid core and finally to complete melting. To deal with this issue properly would require to consider a three phase equilibrium. From the simulation results for T > 4500 K, using Eq. 8 and the previously found time calibration factor t = 5.56 10−4 ps/cycle, one can directly determine the attempt frequency νref . We find νref = 2.3 1014 Hz, which is quite reasonable when comparing it to kB Tref /h = 0.94 1014 Hz. VII. SUMMARY AND CONCLUSIONS By analyzing systematic simulations we have shown that the melting of pristine bulk graphene, a prototype 2D solid embedded in 3D space, follows a two-stage scenario of which the first step from the 2D solid to a quasi2D liquid phase is well described by nucleation theory for first order phase transitions. As a consequence, graphene has an unambiguously defined bulk melting temperature Tm , found to be 4510 K using the interatomic potential LCBOPII. This value is confirmed by simulations for finite disks and ribbons, which by extrapolation yield Tm ’s close to that from nucleation theory, and show reversibility to some extent. Our finding that Tm is higher (about 250 degrees) than for graphite is not likely to be an artifact of LCBOPII. As a qualitative explanation we suggest that the increased stability of graphene is due to a significant positive entropy contribution for the solid phase due to the possibility of relatively large out-of-plane fluctuations (rippling) at low energy costs, as compared to the situation in graphite. We expect that the melting scenario found for graphene applies to any covalent 2D material. This expectation is supported by recent melting simulations for 2D MoS2 , showing coexistence of the 2D solid and liquid phase during melting [22], as we observed for graphene. While measuring Tm for graphene seems to be a real challenge, measurement of Tm becomes easier for 2D material with a (much) lower Tm than that of graphene. Acknowledgements.The research leading to these results has received funding from the European Union Seventh Framework Programme under grant agreement 6 n604391 Graphene Flagship and from the Foundation for Fundamental Research on Matter (FOM), which is part of the Netherlands Organisation for Scientific Research (NWO). We thank Luca Ghiringhelli for useful discussions. VIII. APPENDIX: FORMATION ENERGY OF STONE-WALES AND 57 DEFECTS. A 5775 Stone-Wales (SW) defect can be created in a perfect graphene layer by rotating a single bond by 90 degrees, followed by a geometrical relaxation of the system (see Fig. 6a). It consists of two adjacent 57 defects. By subsequently rotating another appropriate bond by 90 degrees (and relaxing the system), the 5775 SW defects can be split into two 57 defects separated by one [1] A. C. Mitu and A. Z. Patashinskii, Phys. Lett. A 87, 179 (1982); D. R. Nelson, Phys. Rev. B 28, 5515 (1983); M. Kleman, Adv. Phys. 38, 605 (1989); J. F. Sadoc (Editor), Geometry in Condensed Matter Physics (World Scientific, Singapore, 1990); M. I. Katsnelson and A. V. Trefilov, Phys. Met. Metallography 91, 109 (2001); arXiv:cond-mat/9906402; L. Gómez, A. Dobry, Ch. Geuting, H. T. Diep, and L. Burakovsky, Phys. Rev. Lett. 90, 095701 (2003) [2] M. I. Katsnelson, A. Fasolino, Acc. Chem. Res. 46, 97105 (2013). [3] K. V. Zakharchenko, A. Fasolino, J. H. Los and M. I. Katsnelson J. Phys.: Condens. Matter,23, 202202 (2011). [4] U. Gasser, C. Eisenmann, G. Maret, P. Keim, ChemPhysChem 11, 963 (2010). [5] J. M. Kosterlitz, D. J. Thouless, J. Phys. C 6,1181 (1973). [6] B. I. Halperin, D. R. Nelson, Phys. Rev. Lett 41, 121 (1978). [7] D. R. Nelson, B. I. Halperin Rev. B 19, 2457 (1979). [8] A. P. Young Phys. Rev. B 19, 1855 (1979). [9] K. Chen, T. Kaplan and M. Mostoller, Phys. Rev. Lett. 74, 4019 (1995). [10] G. Blatter, M. V. Feigel’man, V. B. Geshkenbein, A. I. Larkin, and V. M. Vinokur, Rev. Mod. Phys. 66, 1125 (1994). pair of hexagons, as shown in Fig. 6b. This process can be repeated to separate the 57 defects more and more, as shown in Fig. 6c and d. The total formation energy as a function of the separation between the (57) defects according to LCBOPII is given in Fig. 7 and shows that it strongly increases with distance. Similar behaviour was previously reported in Ref. [13], but in that case the 57 defects were created at grain boundaries in a polycrystalline sheet. In the latter work it was shown that the formation energies according to LCBOPII are in very good agreement with those calculated ab initio within DFT. The strongly increasing energy implies that the pathway to disorder by the splitting of a SW defects into two 57 defects and their subsequent further diffusion, a scenario in agreement with the KTHNY theory, is very unfavorable. Indeed, we never observed separated 57 defects in our simulations. [11] A. Fasolino, J. H. Los, M.I. Katsnelson, Nat. Mater. 6, 858 (2007). [12] J. Coraux, A. T. N’Diaye, C. Busse, T. Michely, Nano Lett. 8, 565 (2008). [13] J. M. Carlsson, L. M. Ghiringhelli, A. Fasolino, Phys. Rev. B 84, 165423 (2011). [14] J. H. Los, L. M. Ghiringhelli, E. J. Meijer, A. Fasolino, Phys. Rev. B 72, 214102 (2005); 73, 229901(E) (2006). [15] V. M. Bedanov, G. V. Gadyak and Yu E. Lozovik, Phys. Lett. A 109, 289 (1985). [16] F. Colonna, J. H. Los, A. Fasolino, E. J. Meijer, Phys. Rev. B 80, 134103 (2009). [17] D. Kashchiev, Nucleation: Basic Theory and Applications (Butterworth-Heinemann, Oxford, 2000). [18] Huitema, H. E. A., van der Eerden, J. P., J. Chem. Phys. 110, 3267 (1999). [19] S. K. Singh, M. Neek-Amal, and F. M. Peeters Phys. Rev. B 87, 134103 (2013). [20] P. Pavlov, Z. Phys. Chem., Stoechiom. Verwandtschaftsl.65, 1 (1909);65, 545 (1909). [21] J. M. H. Kroes, M. A. Akhukov, J. H. Los, N. Pineau, and A. Fasolino, Phys. Rev. B 83, 165411 (2011). [22] S. K. Singh, M. Neek-Amal, S. Costamagna and F. M. Peeters, arXiv:1412.1939. 7 a b c d e f FIG. 6. Step by step separation of the 5775 SW defect (a) into two 57 defects with one (b), two (c) and nine (d) hexagonon pairs in between. Panels e) and f) are side views of the configuration in a) and d), showing much larger out-of-plane deformations for the case of separated 57 defects, due to strong stress in that case. (total) formation energy (eV) 25 20 15 10 5 0 0 2 4 6 8 separation between 57 defects (in number hexagon pairs) FIG. 7. Total formation energy as a function of the separation between 57 defects according to LCBOPII, calculated for the same sample as used for Fig.6