Molecular dynamics simulation of melting of finite and infinite size silicene

TL;DR

This study uses molecular dynamics simulations with an optimized Stillinger-Weber potential to determine that free-standing silicene melts at approximately 1500 K, providing visual and quantitative analysis of the melting process.

Contribution

First systematic MD simulation of silicene melting using an optimized potential, offering detailed indicators and a consistent melting temperature.

Findings

Melting temperature of silicene is about 1500 K.

Melting scenario is visually captured and analyzed.

Multiple indicators confirm the melting point.

Abstract

We report the melting temperature of free-standing silicene by carrying out molecular dynamics (MD) simulation experiments using optimzed Stillinger-Weber (SW) potential by Zhang {\it et al.}. The melting scenario of a free-standing silicene is well captured visually in our MD simulations. The data are systematically analyzed using a few qualitatively different indicators, including caloric curve, radial distribution function and a numerical indicator known as `global similarity index'. The optimized SW potential consistently yield a melting temperature of 1500~K for the simulated free-standing, infinite silicene.

Equations

$${\xi }_i=\frac{1}{n}\sum _{s=1}^n(k_{s,i}+1{)}^{-1}$$

$$k_{s,i}=|\sqrt{d_{s,i}}-\sqrt{d_{s,0}}|$$

$$\bar{x}=\frac{x_1+x_2+\cdots +x_n}{n}$$

$$\bar{x}=\frac{1}{\frac{1}{x_1}+\frac{1}{x_2}+\cdots +\frac{1}{x_n}}$$

$$\bar{x}=\sqrt{\frac{x_1^2+x_2^2+\cdots +x_n^2}{n}}$$

$${\bar{\xi }}_i=\frac{1}{N_{\mathrm{COR}}}\sum _{\mathrm{COR}}{\xi }_i^{\mathrm{COR}},$$

Full text

Keywords: Molecular dynamics, Melting, Silicene, Optimized Stillinger-Weber potential

1 Introduction

Silicene, a two-dimensional nanosheet made up of silicon atoms arranged in a honeycomb lattice, has been predicted theoretically by Takeda and Shiraishi [1] in year 1994. Subsequent DFT calculations such as [2] and [7] revive the interest of the prediction by showing that silicene is indeed energetically stable, and a feasible possibility. Silicene, unlike grapehene which prefers sp2 hybridization, is not flat. Rather, due to the preference of the bounding for sp3 hybridization, the silicene sheet has a buckled configuration, where the out-of-plane buckle parameter is predicted to be 0.44 Å according to DFT calculations. Having a close resemblance to graphene, silicene offers many possibilities as a functional material of advanced applications, such as photovoltaic, and optoelectronic devices [3], thin-film solar cell absorbers beyond bulk Si [4] and hydrogen storage [5]. One advantage of silicene over other 2D materials is that it is in principle easier to get integrated into nanodevices which are mainly silicon-based.

Silicene is a relatively new 2D material, and was first synthesized on supported substrates in a series of discovery since 2007 [2]. Following the successful synthesis of silicene on supported substrate, many theoretical studies and simulations on the structural, mechanical, electronic and thermal properties of silicene on supported substrate have been published [6].

The structural properties of a free-standing silicene sheet get modified when grown on a substrate. The silicene experimentally synthesized so far are not free-standing but sitting on a substrate. We have yet to see any report of experimentally synthesized free-standing silicene. Having said that, investigation of free-standing silicene serves the purpose to understand the pristine system in the absence of interactions with surfaces. The understanding of the basic properties of silicene without the interference from surface interaction with substrate shall provide useful insight for higher level manipulation of silicene, such as the ‘van der Waals’ heterostructures envisaged by [8].

As far as we are aware of, very limited work on the melting behavior and thermal stability of free-standing silicene is reported in the literature. Bocchetti et al. simulated the melting behavior of free-standing silicene via Monte Carlo method with original and a modified version of Tersoff potential parameter set (known as ARK) for silicon atom [10]. According to Bocchetti et al., original Tersoff parameters for silicon atom results in a melting of the free-standing silicene at 3600 K, meanwhile the melting temperature obtained using ARK parameter set is $\sim$ 1750 K. Berdiyorov et al. simulated the influence of defect on the thermal stability of free-standing silicene via MD using Reactive force-field (ReaxFF) [9]. He found that pristine silicene is stable up to 1500 K. As a general observation, melting properties and thermal stability of free-standing silicene obtained in MD simulations varies from cases to cases depending on the details of the simulation procedure. Most importantly, the simulation results are strongly forcefield-dependent. Apart from simulating thermal stability, MD simulation has been also been applied to investigate or predict thermal conductivity of both free-standing silicene. A wide range of potentials is employed in these simulations. For example, Zhang et al. developed a set of Stillinger-Weber potential parameters specifically for a single-layer Si sheet to simulate the thermal conductivity [11]. Others use Tersoff potential with original parameter sets. In this paper we wish to contribute to the existing knowledge about the melting temperature of free-standing silicene sheet by carrying out MD simulations using modified Stillinger-Weber (SW) [11] potentials. Our MD simulations are carried out using the LAMMPS package [12]. To our best knowledge, optimized SW forcefield has never been used for simulating the melting scenario of an infinite silicene sheet.

2 Methodology

To obtain an infinitely large free-standing silicene sheet, we first construct a bulk Si structure using a diamond unit cell with a lattice constant of 5.431 Å. This is done by replicating a total of $30\times 30\times 5$ diamond unit cells to form a simulation box that is subjected to periodic boundary condition. The structure is then re-oriented such that the (1 1 1) surface is aligned along the $z$-axis. This is done because the Si diamond lattice in the (1 1 1) direction resembles a silicene sheet with a non-zero buckling parameter. After re-orientation, all atom-containing $xy$-planes perpendicular to the $z$-axis are removed from the bulk until only a single sheet is left. A total of 72000 atoms is removed in the process. The removal of the silicon atoms creates a region of vacuum with 41.5 Å thick along the $z$-direction on both sides of the remaining silicon atom plane, which is now resembling a free-standing silicene sheet with surface area of $130.56\times 150.76$ Å2, consisting of 9600 atoms. The sheet is infinite in extent (in the $x$- and $y$-directions) due to periodic boundary condition imposed in the MD simulation.

The data file of the free-standing silicene sheet is then inserted into the LAMMPS package to determine its melting point by carrying out simulated annealing simulations. SW potential is well known for its wide usage in semiconducting materials. It is first developed in 1985 for various phase of silicon crystal and it is in good agreement in predicting various physical properties regarding silicon crystals. The formalism of SW and its detail parameters is referred to the original Stillinger et al. paper [13]. As it turns out, during our trial attempts to perform silicene melting simulation, conventional SW parameters of [13], as implemented in the LAMMPS codes, do not produce stable free-standing silicene sheet. On the other hand, using modified version of silicon SW potential by Zhang et al. [11] (the so-called optimized SW), we are able to obtain a stable 2D honeycomb structured silicene sheet upon energy-minimization at temperatures below its melting point. Moreover, the sheet turns out to have a buckling value of 0.44 Å out of the $xy$-plane. The average bond length of Si-Si after energy minimization is 2.4 Å. Hence, the optimized SW silicon potential by [11] is used exclusively in the present work. The details of the procedure are summarized as follow:

The free-standing silicene is relaxed via conjugate gradient minimization (energy minimization) technique. The resultant structure is depicted in Fig. 1. 2. 2.

Simulated annealing procedure is proceeded in NVT ensemble using the Nose-Hoover thermostat with a time step $\Delta t=0.5$ fs. The system is then equilibrated at 1 K for a period of $5\times 10^{4}\Delta t$. Temperature is then raised to 300 K at a heating rate of $1\times 10^{8}$ K/s. After reaching 300 K, the system is equilibrated for a period of $2.5\times 10^{4}\Delta t$. 3. 3.

The stage is now set for raising the system’s temperature to a large target temperature $T$ (one at which the silicene will surely melt). To raise the temperature from 300 K to $T$, we have to fix a heating rate. We have performed a convergence test by running trial MD melting simulations to obtain the melting temperature of the free-standing silicene at various heating rates by fixing the target temperature $T=2000$ K. From our convergence test, it is found that by choosing a heating rate of 1.7$\times 10^{12}$ K/s, resultant melting temperature obtained is insensitive to variation in the heating rate pegged at the chosen value. The system is heated from 300 K for a total $2.0\times 10^{6}\Delta t$ (equivalent to $1$ ns) at the mentioned heating rate until it reaches 2000 K. 4. 4.

Evolution of the silicenes’s MD trajectory in the simulation is monitored visually as well as quantitatively. Radial distribution function, and caloric curve of the silicene sheet are numerically sampled and measured. In addition, we also measure a numerical descriptor, known as ‘global similarity index’, $\xi_{i}$, to gauge the melting process. Detailed description of the global similarity index is discussed in the following.

The definition of the global similarity index is based on generic chemical similarity idea for detecting configurational changes along the trajectory during the heating process. The evolution of the global similarity during the simulated annealing process, which carries statistical information of geometrical variations in the melting mechanism of the silicene will be elaborated in Section 3. The functional form of $\xi_{i}$ is proposed to take the form of

[TABLE]

[TABLE]

where $d_{s,i}$ and $d_{s,0}$ represent the sorted distance of atoms relative to the average positions (center of mass) of all the atoms in the system in the $i$-th (denoted as the subscript $i$) and the 0-th “frame” during a MD simulation, while $n$ corresponds to the number of atoms, which is an integer equals to the number of pairs of $k_{s,i}$. The value of $\xi_{i}=1$ corresponds to total identicalness between the configuration of the system in the $i$-frame as compared to that in the zeroth frame (which in practice is chosen to be a reference frame in which the system appears in a relatively perfect and stable configuration). $\xi_{i}\rightarrow 0$ infers vast difference.

Having defined the similarity index in Eqs. (1), (2), one can optimize the sensitivity of the index with respect to degrees of freedom corresponding to few parameters such as translation and rotation of the particles relative to one another. Note that the average positions, a.k.a., mean, that enters the definition of the parameter $d_{s}$ can be non-uniquely defined. Different definitions of mean capture different aspects of configurational information contained in the system. Several definitions of mean are possible, namely, (a) arithmetic mean [see Eq. (3)], (b) harmonic mean [see Eq. (4)] and (c) quadratic mean [see Eq. (5)].

[TABLE]

[TABLE]

[TABLE]

However, during a melting process, variations in the configuration of the atoms a priori are not known. In principle, variation in a certain mode of motion among the atoms could be more sensitively picked up by a particular definition of mean then the other. To cover all possibilities, the parameter $d_{s}$ that enters the definition of $\xi_{i}$ in Eqs. (1), (2) is calculated by averaging overall of the above three types of mean, i.e.,

[TABLE]

where $\mathrm{COR}$ stands for center of reference, $\mathrm{COR}$=$\{$arithmetic mean, harmonic mean, quadratic mean$\}$, $N_{\mathrm{COR}}=3$, $\xi_{i}^{\mathrm{COR}}$ is global similarity index defined based on a $\mathrm{COR}$ mean. Including all types of $\mathrm{COR}$ can in principle increase the sensitivity as well as accuracy of the index in capturing variations in the geometrical configuration of the system. Similarity index of the silicene frames in the MD simulation during its heating process is captured in successive temporal sequence. If the silicene sheet remains intact, its geometrical configuration should be very close to the one at the initial frame, i.e., ${\bar{\xi}_{i}}\approx 1$. If the silicene melts, we would expect ${\bar{\xi}_{i}}$ to deviate from 1 and in the extreme case, becomes 0. In the following we shall drop the bar symbol in the definition of $\xi_{i}$ when the index is referred. As it turns out, chemical similarity is rarely used for very large system [14] In the present case, the system being simulated consists of 9600 Si particles. This is considered a moderately large one. Global similarity index (which is a modified form of chemical similarity index) as defined above has a high sensitivity to detect even minor distortions in the geometrical configuration of silicene, and hence a very convenient tool to pinpoint the location of melting point during the temporal evolution of the MD.

It is commented that caloric curve and the other two quantities, $g(r),\xi_{i}$, are two qualitatively different indicators for gauging the melting process. The former is a thermodynamical quantity whereas $g(r),\xi_{i}$ embed only geometrical information of how atoms distribute in 3D space.

3 Results and discussion

We first present the results of pair-correlation function, $g(r)$, measured in the simulation. Each pair-correlation function shown in Fig. 2 is obtained by sampling the data in the MD simulation at that specific temperature. At room temperature, a sharp peak is measured at $\sim$2.5 Å, indicating an average bond length of silicene at this value. As temperature increases, the peak at 2.5 Å is lowered and begins to widen.

The $g(r)$ curves with $T\leq 1400$ K are clearly distinctive from those with $T\geq 1500$ K. The former has their first peaks located at 2.5 Å, while subsequent peaks are also clearly visible. In addition, $g(r)$ is approximately zero between the first and second peaks. This group of curves essentially share a qualitatively similar geometrical distribution. We hence conclude that silicene sheet with temperature $T\leq 1400$ K remains intact and has not shown signs of complete melting.

However, for curves with $T\geq 1500$ K, the first peak shifted to a larger position at $\sim$2.8 Å, and the secondary peaks disappear. In addition, $g(r)$ beyond the first peak levitated to non-zero values, and remains relatively flat and constant beyond $\sim 3.0$ Å. The disparity between these two groups of curves indicates that a qualitative change in the geometrical distribution of the silicene sheet has occurred when temperature crosses the 1400-1500 K benchmark. At 1500 K melting process has taken place. The features of the $g(r)$ curves for $T\geq 1500$ K are indicative evidence that the silicene sheet has melted into a liquid state at $T\geq$1500 K. The liquid state of the melted silicon system beyond 1500 K can be visualised in the last few snaphots in Fig. 4. These are liquid blobs made up of silicon atoms sticking together with an average bond length of $\sim$ 2.8 Å.

Figures shown in Figs. 3 and 4 visually depict progressive melting process of the silicene sheet in the simulated annealing simulation obtained when the system is heated from 300 K towards 2000 K. At room temperature, the silicene sheet is thermally stable and maintains its hexagonal shape at 300 K. It is not perfectly flat though, with some puckering distortion along the direction perpendicular to the surface if viewed from the side. The equilibrium bond length of Si-Si is 2.5 Å. At room temperature, the Si particles vibrate at their respective equilibrium positions, maintaining a defect-free silicene sheet. As the temperature increases, the particles begin to vibrate more violently and vigorously around $\sim$1100 K. Visually, it can be observed that vacancies begin to appear at selected areas on the silicene sheet as Si particles are ejected from their existing sites due to the thermal energy gained from the reservoir. At this temperature, however, the structure remains overall intact. For temperature at or above 1500 K, disintegration of the silicene sheet is clearly sighted. For temperatures less than the melting temperature, no sighting of silicene melting is observed. In other words, visual sighting provides another direct indication of melting at and above $T=1500$ K.

Caloric curves sampled for the simulated silicene sheet at various temperatures provide an independent measurement of the melting point. It is seen that in Fig. 5 an abrupt drop in the potential energy plot occurs at $\sim$1500 K. From here onward, the sheet begins to tear apart and lose molecular cohesion, signaling the onset of full-fledge melting process. Concurrent to the abrupt drop in the caloric curve, disintegration of the sheet at $\sim$1500 K is also visually sighted in MD simulation video. The atoms of the silicene sheet break apart, forming disconnected 3D amorphous blobs, such as that displayed in the last few snapshots in Fig. 4.

The melting process is now analyzed in terms of global similarity index, which is yet another independent measurement of the melting temperature. The structure during heating process is frame-by-frame compared with the reference structure at 300 K (the ‘0-th frame’). The curve of global similarity index as a function of temperature is shown in Fig. 6. As temperature increases, the index gradually decreases. At 1500 K the global similarity index drastically drops from $\sim 1$ to 0.02, inferring that the configuration of the system at or above 1500 K is drastically different from that at the reference configuration at 300 K. The melting point as indicated by global similarity index coincides with that obtained from other gauging indicators discussed above.

At this point we wish to comment that the emergence of local point vacancies, visually sighted at $\sim 1100$ K, when silicon atoms leaving their existing sites, cannot be effectively captured by all melting indicators.

The final structure obtained when the system is quenched to 300 K is as shown in Fig. 7. It is found that the resultant structure does not revert to its original planar form. Instead, these atoms condensed into some kind of 3-D amorphous “columns”.

We have attempted to run another round of independent MD simulation on the melting of silicene sheet, with all the simulation specifications being exactly the same as that reported but with the NVT ensemble replaced by the NPT ensemble. It is reckoned that the implementation of NPT ensemble in a MD simulation effectively mimics an ‘adjustable periodic boundary’. The volume of the simulation in the NPT ensemble is monitored as a function of temperature. The variation in the volume as the temperature increases is as depicted in Fig. 8. It is found that the volume of the silicene fluctuates along the course where the temperature is annealed from below toward the melting temperature. Once it reaches the melting temperature the volume displays an abrupt drop, resulting in an apparent decrease in the volume of the silicene simulation. This observation is consistent with that found in bulk silicon melting process [15]. It is found that the melting temperature as measured in the NPT and NVT ensemble remains consistently similar. The only slight observational difference of the simulation in both ensembles is that the silicene in the NPT appears to display an undulation of a relatively larger amplitude before the melting point. In both ensembles, the top view of silicene reveals that the sheet remains intact in its intrinsic hexagonal honeycomb shape before melting.

4 Conclusions

We have performed a MD experiment to measure the melting point of a finite but large free-standing silicene sheet using optimized SW potential by Zhang et al.. The melting scenario of the silicene sheet has been visually captured. Analysis of the statistical data sampled from the MD simulation, including caloric curve, radial distribution function and global similarity index, reveals that the silicene sheet melts at 1500 K. This result serves the purpose of being an independent measurement of the melting temperature of free-standing silicene. Our value is consistent with that obtained by Berdiyorov et al who uses ReaxFF [9], but is in contrast to [10] which predicts 1750 K using ARK parameter set. The prediction of free-standing silicent sheet will serve as a reference information to be verified by experiments when they become available.

List of Figures

1. 1 Left: An infinitely large free-standing silicene sheet can be created by cutting out a silicon atom plane from the bulk Si diamond structure after re-orientating the (1 1 1) surface of the bulk into the $z$-axis. Right: The resultant structure has a 2D honeycomb shape with a dimension of 130.56 Å$\times$150.76 Å and buckling parameter of 0.44 Å. The average bond length of Si-Si is 2.4 Å.
2. 2 Pair correlation function of the system at various temperatures. Two groups of curves are clearly visible. The group comprises of 300, 1100, 1200, 1300 and 1400 K curves (continuous, colored lines) is qualitatively distinctive from the group comprises of 1500 - 2000 K curves (filled and open circles and triangles).
3. 3 Simulated annealing of silicene by using optimized Stillinger-Weber potential. Initially the silicene has equilibrium bond length of 2.5 Å. An abrupt change can be viewed at 1500 K in which the sheet is tearing apart and the melting process has begun. After the melting point, the Si particles settle down to form four smaller “islands” which has equilibrium bond length of 2.8 Å.
4. 4 Continuation from Fig. 3.
5. 5 Potential energy plot (caloric curve) of the system (eV) against temperature (K) measured for target temperature $T=2000$ K. The arrow indicates the point at which particles begin to vibrate more violently and vigorously such that vacancies begin to appear at selected areas on the silicene sheet as Si particles are ejected from their existing sites. The sharp and abrupt drop of potential energy at 1500 K indicates occurrence of melting at that temperature.
6. 6 Global similarity index plot against temperature. The silicene structure is compared with is 300 K state. The value of 1 means both comparing structure are identical. As the temperature increases, the comparing structures become more and more dissimilar as the kinetic energy of system increases, the Si particles vibrate (there maybe rotational or translational motion based on statistical influence of choice) more violently. The melting occurs at 1500 K, and the structures is immediately unable to identify as a silicene as compare to its 300 K structure.
7. 7 Amorphous silicon “cylinder” after quenching of the melted silicene. The structure is unable to revert to a sheet form.
8. 8 Variation of simulation volume as a function of temperature in the NPT ensemble during the melting of the silicene sheet.