Bond-coordination lattice model for phase transformations in carbon

Diamond & Related Materials 20 (2011) 1310–1314

Contents lists available at SciVerse ScienceDirect

Diamond & Related Materials journal homepage: www.elsevier.com/locate/diamond

Bond-coordination lattice model for phase transformations in carbon M.S. Byshkin ⁎ Akhiezer Institute for Theoretical Physics within National Science Center, “Kharkov Institute of Physics and Technology”, 1 Akademicheskaya str, Kharkov 61108, Ukraine

a r t i c l e

i n f o

Article history: Received 16 September 2010 Received in revised form 5 July 2011 Accepted 25 August 2011 Available online 3 September 2011

a b s t r a c t A lattice model for phase transformations in pure carbon is proposed. Monte Carlo simulation is performed. The simulation results show that the proposed model produces phase transformations similar to phase transformations graphite–liquid, diamond–liquid, graphite–diamond and liquid–liquid. © 2011 Elsevier B.V. All rights reserved.

Keywords: Phase transformation Carbon Lattice model Simulation

1. Introduction Carbon attracts considerable interest of researchers due to the wide structural variety of its pure solid forms having unique properties. In addition to the crystal form (graphite, diamond), carbon can be found in various amorphous and nanostructural forms [1,2]. At a temperature below 1000 K the structure relaxation in amorphous carbon is very slow and many forms of amorphous carbon are stable enough to be used in different applications. The diamond-like amorphous carbon attracts particular interest of researchers. This form of amorphous carbon possesses properties similar to those of diamond and has a high concentration of atoms in the sp 3 configuration. The pressure-induced phase transformations were observed in amorphous silicon and germanium. The behavior of these amorphous samples mimics the densitydriven phase transformation that is expected to occur in amorphous liquids above the glass transition [3]. The pressure (stress)-induced phase transformation between diamond-like and graphite-like carbon was proposed by McKenzie [4]. The phase transformation in liquid carbon [5] and amorphous carbon [6–11] is still debated. The MD study of solid and liquid carbon at high pressure [12–14] has shown that there occur transitions between several high-pressure phases. The phase diagram was proposed. The changes of liquid's character from about 4-fold to about 6-fold coordinated was observed in this MD study within the range from 400 to 1000 GPa. A negative slope of the diamond–liquid coexistence line was observed for a high pressure value as in the phase diagram of silicon. However, the properties of real carbon at these pressure values are almost unknown because to study them experimentally one needs extreme conditions.

⁎ Tel.: + 380 57 3356203; fax: + 380 57 3352683. E-mail address: [email protected]. 0925-9635/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.diamond.2011.08.010

In pure carbon forms the carbon atoms are usually found in three hybridizations [1,2]. In the sp 3 configuration, as in diamond, four valence electrons of carbon atoms are each assigned to the tetrahedrally directed sp 3 orbital, which makes a strong σ bond with an adjacent atom. In the three-fold coordinated sp 2 configuration as in graphite, three of four-valence electrons enter the trigonally directed sp 2 orbital, which form the σ bonds in a plane. The fourth electron of sp 2 atom lies in the pπorbital, being normal to the σ bonding plane. In the sp 1configuration two- of four-valence electrons enter the sp 1 orbitals, and two other enter thepπ orbitals. The π bonds are much weaker than the σ bonds. Under normal conditions the energy of the σ bond between atoms in the sp 3 configuration in diamond is near − 3.6 eV, and the bond length is 0.155 nm [15]. The energy of the σ bond between atoms in the sp 2 configuration in graphite is of about −4.7 eV and the bond length is 0.142 nm. Thus, the σ bond in graphite is stronger than that in diamond, but the bonding between graphite layers is weak.

2. Model Neglecting the energy of weak π bonds, the interaction energy of carbon atoms can be written as a sum of all the σ bonds. In the real carbon all the σ bonds have different energies. In crystalline carbon the temperature fluctuations of bond lengths and bond energies occur. In the disordered structure the bond energies and bond lengths differ more essentially. Nevertheless, as well as in crystals the local order in disordered structures is determined by the interaction of valence atomic electrons, and the local order near every atom is determined mainly by its coordination number [7]. Hence, the energy of σbond of two atoms strongly depends on the coordination numbers of these atoms.

M.S. Byshkin / Diamond & Related Materials 20 (2011) 1310–1314

a

Based on the above considerations the following model is proposed. The interaction energy of carbon atoms U(r1, r2, …, rN) is written as U ðr1 ; r2 ; …; rN Þ ¼ ∑ Eci cj i;j

1311

Energy disribution

1

ð1Þ

0.5

where ri is the atom position, Elm is the bond energy between l-fold and m-fold coordinated atoms, ci is the coordination number of atom i, the summation is over all pairs of atoms, bound by the σ bond, i.e. over all the neighboring atoms. 3. Bond-order potential

0 -6

The bond-order potential is a class of empirical potentials used in the molecular dynamics and Monte Carlo simulations. The examples include the Tersoff potential [16], the Brenner potential [17], the second-moment tight-binding potentials [18] and others. These potentials are written in the form

i;j

-2

Between sp2 and sp2 atom Between sp3 and sp3 atom Between sp2 and sp3 atom

b U ðr1 ; r2 ; ::; rN Þ ¼ ∑ Eij ;

-4

Bond energy

ð2Þ

Energy disribution

1

where U is the potential energy of atomic interaction, Eij is the energy of bond between atoms i and j,

 Eij ¼ fR rij þ bij ·fA rij

0.5

ð3Þ

rij is the distance between atoms i and j, the summation is over all pairs of atoms which are spaced within the first coordination sphere, i.e. neighboring atoms. bij is the bond-order term. Thus, the potential is written as a simple pair potential but the strength of this bond is modified by the environment of the atom i via the bij. Several bond-order potentials were proposed for modeling covalently bonded materials and carbon in particular. The first successful bond-order potential was the empirical potential of Tersoff [16]. Been proposed 20 years ago it is still used by many authors. In the Tersoff potential the bond-order term is given by !


 bij ¼ φ ∑ g θijk ;

0 -6

-4

-2

Bond energy Between sp2 and sp2 atom Between sp3 and sp3 atom Between sp2 and sp3 atom Fig. 1. Partial bond energy distributions in the liquid carbon structure under pressure of 10 GPa at a temperature of 4500 K (a) with temperature fluctuations and (b) structure is artificially deprived of temperature fluctuations.

k≠i;j

shown in Fig. 2. It has been found that E34 is well described by expression θijk is the angle between the atomic pairs ij and jk, the summation is over the neighbors of the atom j. This potential and the functions f R, f A, g and φ, were chosen to fit the theoretical and experimental data obtained for realistic and hypothetic silicon [19] and carbon [16] configurations. MD simulation was performed using Tersoff potential. The liquid carbon structures under pressure ranging from 10 GPa to 70 GPa at a temperature of 4500 K were obtained. The details of simulation may be found in [8,9]. Under these pressure values practically all the atoms in the structure were either 3-fold or 4-fold coordinated. The bond energies Eij were studied as a function of the coordination numbers of atoms i and j. Fig. 1 shows a normalized bond energy distribution. Three types of bonds were considered: 1) bond between two 4-fold coordinated atoms, 2) bond between two 3-fold coordinated atoms, and 3) bond between 3-fold and 4-fold coordinated atoms. From the bond energy distribution one can see that the interaction energy of two atoms Eij is in strong dependence on the coordination numbers of the atoms ci and cj. Thus, the energy of atomic interaction (Eq. (2)) may be approximately written in the form of Eq. (1). The mean bond energy between two 4-fold coordinated atoms E44, mean bond energy between two 3-fold coordinated atoms E33 and mean bond energy between 3-fold and 4-fold coordinated atoms E34 were determined for the liquid carbon structure at a temperature 4500 K. The pressure dependences of E33, E34 and E44 are

E34 ¼

E33 þ E44 þ Ep; 2

with Ep=0.05 eV. The values of the partial bond energies E44, E33 and E34 in the amorphous carbon structure were not determined as the relaxation time is very long in glass. The main result of this section is that the energy of bond between two atoms is in strong dependence on their coordination numbers. It should be noted, that the results presented in Figs. 1 and 2 concern the model carbon with the Tersoff empirical bond-order potential. Nevertheless, we believe that they mimic the reality, at least qualitatively. 4. Monte Carlo simulation A system of atoms with interaction energies (Eq. (1)) may be studied on a lattice. The variable aij is defined: aij = 1, if the bond between two neighboring lattice sites i and j exists and aij = 0, if the bond between these sites is eliminated. The energy (Eq. (1)) of such a system may be written as U ¼ ∑ Eci cj ·aij ; i;j

ð4Þ

1312

M.S. Byshkin / Diamond & Related Materials 20 (2011) 1310–1314

Mean bond energy, eV

a

and the variables {aij}, obtained for one parameter set, were used as initial variables for a new parameter set. The model (4) was used for MC simulation on the diamond lattice having the size of 60 3 sites. The following model parameters were used: T = 4500 K and bond energy values E33(P), E44(P), E34(P) in liquid carbon presented in Fig. 2. The obtained pressure dependences of energy U and of concentration of 4-fold coordinated sites c4 are presented in Fig. 3. All the partial bond energies Elm except E33, E34, E44 were assigned zero, that is sufficiently higher than the values of E33, E44, E34. This caused all the coordination numbers of sites in the equilibrated structures to be either 3 or 4. The stepwise increase of 4-fold coordinated site concentrations under pressure of 30 GPa is observed in Fig. 3a. The molecular dynamic simulation with the Tersoff potential of amorphous carbon [8] has shown the peculiarities of the pressure dependences of amorphous carbon properties at the same pressure value, about 30 GPa. The stability of diamond crystallites being embedded in supercooled carbon liquid was studied [9]. The pressure value of 30 GPa was found to be the critical pressure for the crystallites stability at 4500 K. Thus, the peculiarities described in these papers are related with the stepwise increase of the concentration c4 of 4-fold sites in the model (4). The peculiarities of amorphous carbon properties are observed in [8,9] when 3E33 = 4E44. The step-wise change of c4 in Fig. 3 is observed exactly when 3E33 = 4E44 (see Fig. 2) and this property of model (1) is temperature independent. The temperature dependence of properties of the system described by Hamiltonian (Eq. (4)) was studied. Two sets of the bond energies E33, E44, E34 were taken for simulation: the bond energies in liquid carbon at pressures 10 and 60 GPa, at 4500 K (Fig. 2). The temperature dependences of the concentration of 4-fold coordinated

-4.2 -4.21 -4.22 -4.23 -4.24 -4.25 -4.26 -4.27

0

10

30

20

40

50

60

70

60

70

Pressure, GPa E33 E44*4/3

Mean bond energy, eV

b

-3.62

-3.64

-3.66

-3.68

0

10

20

30

40

50

a

Pressure, GPa E34 (E33+E34)/2+0.05

0.35 4 3 E44,

sp3 fraction

Fig. 2. Partial bond energies as a function of pressure. Circles — E33, boxes — crosses — E34. Dotted curves are the fits to energy values.

T=4500 K

where ci is the coordination number of site i, ci ¼ ∑ aij ;

0.3

0.25

j

8 1; ΔU ≤ 0 > > < 1   ; ΔU N 0 ; p¼ ΔU > > : 1 þ exp kB T

0.2 10

20

30

40

50

60

50

60

Pressure, GPa

b

T=4500 K

-6.32 U, eV

the summation is over all neighboring sites of site i, Elm is the bond energy between l-fold coordinated and m-fold coordinated atoms. Thus the lattice model with Hamiltonian (Eq. (4)) is obtained. The Metropolis Monte Carlo method [20] was used to study the model (4). Each simulation step in this method is described as follows. The random variable aij is chosen and the trial change of its value is made: the bond is being eliminated if it exists and the bond is being set if it is eliminated. The change of structure energy ΔU upon this change is found. The change of aij at this simulation step is accepted with the probability

ð5Þ

where T is the temperature and kB is the Batsman constant. After sufficient number of simulation steps the variable set {aij} and the energy U reach their equilibrium values, which are independent of initial values of {aij} and are defined by the model parameters: bond energies Elm and temperature T. The simulation with different parameter sets was performed. The parameter values were changing smoothly

-6.34

10

20

30

40

Pressure, GPa Fig. 3. The temperature dependence of (a) the concentration of 4-fold coordinated sites c4 and of (b) the energy U of the system at 4500 K.

M.S. Byshkin / Diamond & Related Materials 20 (2011) 1310–1314

a

1

0.5

0

T=1000 K 1

sp3 fraction

sp3 fraction

a

1313

0

1000

2000

3000

4000

0.5

5000

Temperature, K

0 10

E33(P=10 GPa), E44(P=10 GPa) E33(P=60 GPa), E44(P=60 GPa)

20

30

40

50

60

50

60

Pressure, GPa

b

b

T=1000 K

-6.34 -6.35

U, eV

U, eV

-6.32

-6.4

-6.36

-6.38 0

1000

2000

3000

4000

5000

Temperature, K E33(P=10 GPa), E44(P=10 GPa) E33(P=60 GPa), E44(P=60 GPa)

-6.4 10

20

30

40

Pressure, GPa Fig. 4. The temperature dependence of (a) the concentration of 4-fold coordinated sites c4 and of (b) the energy U of the system. Solid curves: E33 = −4.262 eV, E44 = −3.163 eV. Dashed curves: E33 = −4.219 eV, E44 = −3.184 eV.

sites c4 and of the energy U are presented in Fig. 4. The stepwise changes of these properties at a temperature of about 1500 K are observed, the hysteresis is observed too, and thus the explicit indications of phase transformations are found. At low temperature and low pressure all the sites are 3-fold coordinated similarly as in graphite. At low temperature and high pressure all the sites are 4-fold coordinated similarly as in diamond. Thus the analogy with real carbon is observed. As temperature is increased the phase transformations similar to phase transformations graphite–liquid and diamond–liquid occur. The MC simulation at a temperature of 1000 K was made. The bond energies values are those presented in Fig. 2. The results are shown in Fig. 5. Here the indications of phase transformation are also observed. One can see that as the pressure is increasing the phase, all sites of which are 3-fold coordinated, transforms to the phase, all the sites of which are 4-fold coordinated. Thus, the analogy with graphite–diamond phase transformation is observed. Note that the values of bond energies E33, E44, E34 in real carbon at a temperature other than 4500 K have values different from those used for MC simulation and presented in Fig. 2 (Figs. 4 and 5 are presented only for the qualitative study of the model proposed in this paper). The MC simulation of the system with carbon atoms having interatomic energy (Eq. (1)) was also performed in the pure form without using the lattice model. A continuous-space Monte Carlo study of the system was carried out: the trial atom displacement made on each step of simulation was accepted with a probability (Eq. (5)). The energy U was given by Eq. (1) with the summation over all pairs of atoms which are within 1.8 Å from each other. All the results of these simulations are in qualitative agreement with results presented in Figs. 3–5. These simulations have taken a long time. In contrast, the results presented in Figs. 3–5 were calculated in rather short period of time, approximately 10 min. The interaction energy

Fig. 5. The pressure dependence of (a) the concentration of 4-fold coordinated sites c4 and of (b) the energy U of the system at 1000 K.

(Eq. (1)) has no length and angle terms and the structures obtained with this model have no structure properties inherent to real carbon. The results presented in Figs. 3–5 concern the proposed lattice model (4) with partial energies Elm obtained from the Tersoff potential model. The value of the pressure at which transformations are observed on Fig. 3 is the same as that found using the Tersoff potential, but is overestimated as compared to both Tight Binding model and real experiments [8]. It would be interesting to see the results of model (4) simulation with partial bond energies obtained from more accurate model. A possible range of model applications should be studied too. For example, any evidence of phase transformations was not found in the results of MD simulation of liquid carbon at a temperature more than 5000 K [9]. Furthermore, it is considered that the bond-order potentials are not useful for modeling carbon at pressure as high as 400 GPa and hence the proposed model is hardly adequate for modeling carbon at very high pressure values. In the present paper the lattice model is proposed for modeling the phase transformations in pure carbon under not extreme conditions.

5. Conclusions The energy of σ bond between atoms in liquid carbon was studied by means of molecular dynamics simulation with the Tersoff potential. It is shown that bond energy between two atoms is in strong dependence on their configurations. The pressure dependences of partial bond energies are found. A simple model for the interaction energy of carbon atoms is proposed: the potential energy of atom interaction is completely defined by the coordination numbers of neighboring atoms. The Monte Carlo simulation has shown that the model proposed produces phase

1314

M.S. Byshkin / Diamond & Related Materials 20 (2011) 1310–1314

transformations similar to the phase transformations such as graphite– liquid, diamond–liquid, graphite–diamond and liquid–liquid. References [1] [2] [3] [4] [5] [6] [7]

J. Robertson, Mater. Sci. Eng. 37 (2002) 129. R.C. Powles, N.A. Marks, D.W.M. Lau, Phys. Rev. B 79 (2009) 075430. P.F. McMillan, J. Mater. Chem. 14 (2004) 1506 and references therein. D.R. McKenzie, D. Mulller, B.A. Pailthorpe, Phys. Rev. Lett. 67 (1991) 773. J.N. Glosli, F.H. Ree, Phys. Rev. Lett. 82 (1999) 4659. D.G. McCulloch, D.R. McKenzie, C.M. Goringe, Phys. Rev. B 61 (2000) 2349. A.S. Bakai, M.P. Fateev, Yu.A. Turkin, in: G. Benedek, P. Milani, V.G. Ralchenko (Eds.), Nanostructured Carbon for Advanced Applications, Kluwer Academic Publisher, Netherlands, 2001, pp. 185–198. [8] M.S. Byshkin, A.S. Bakai, A.A. Turkin, Condens. Matter Phys. (in press).

[9] M.S. Byshkin, A.S. Bakai, A.A. Turkin, Diam. Relat. Mater. 19 (2010) 1058. [10] D.W.M. Lau, D.G. McCulloch, M.B. Taylor, J.G. Partridge, D.R. McKenzie, N.A. Marks, E.H.T. Teo, B.K. Tay, Phys. Rev. Lett. 100 (2008) 176101. [11] Y.X. Wei, R.J. Wang, W.H. Wang, Phys. Rev. B 72 (2005) 012203. [12] M.P. Grumbach, R.M. Martin, Phys. Rev. B 54 (1996) 15730. [13] A.A. Correa, S.A. Bonev, G. Galli, PNAS 103 (2006) 1204. [14] J. Sun, D.D. Klug, R. Martonak, J. Chem. Phys. 130 (2009) 194512. [15] S.N. Marinkovich, in: L.R. Radovic (Ed.), Chemestry and Physics of Carbon, vol. 29, Marcell Dekker, New-York, 2004, p. 71. [16] J. Tersoff, Phys. Rev. Lett. 61 (1988) 2879. [17] D.W. Brenner, Phys. Rev. B 42 (1990) 9458. [18] D.G. Pettifor, I.I. Oleinik, Prog. Mater. Sci. (2004) 285. [19] J. Tersoff, Phys. Rev. B 37 (1988) 6991. [20] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, E. Telleret, J. Chem. Phys. 21 (1953) 1087.