Abstract
Despite silicon being of great technological importance, an understanding of its behavior across the phase diagram is still lacking, especially near liquid-solid coexistence. The difficulty in describing silicon near coexistence from first principles lies in discriminating between the metallic and covalent bonds present in the material. Using the strongly constrained and appropriately normed (SCAN) density functional, which can describe a wide variety of bonds with quantitative accuracy, we report a thorough investigation of liquid silicon in the vicinity of liquid-solid coexistence using ab initio molecular dynamics simulations. We observe a structural transition in the supercooled regime that is rooted in a change in the electronic structure of the material. This transition is found to occur at a higher temperature than previous predictions. We also discuss implications of the observed change in interatomic interactions for empirical models of transitions between two distinct liquids.
| Original language | English (US) |
|---|---|
| Article number | 140103 |
| Journal | Physical Review B |
| Volume | 97 |
| Issue number | 14 |
| DOIs | |
| State | Published - Apr 26 2018 |
| Externally published | Yes |
All Science Journal Classification (ASJC) codes
- Electronic, Optical and Magnetic Materials
- Condensed Matter Physics
Access to Document
Other files and links
Fingerprint
Dive into the research topics of 'Refined description of liquid and supercooled silicon from ab initio simulations'. Together they form a unique fingerprint.Cite this
- APA
- Author
- BIBTEX
- Harvard
- Standard
- RIS
- Vancouver
}
In: Physical Review B, Vol. 97, No. 14, 140103, 26.04.2018.
Research output: Contribution to journal › Article › peer-review
TY - JOUR
T1 - Refined description of liquid and supercooled silicon from ab initio simulations
AU - Remsing, Richard C.
AU - Klein, Michael L.
AU - Sun, Jianwei
N1 - Funding Information: This work was supported as part of the Center for the Computational Design of Functional Layered Materials, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award No. DE-SC0012575. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy. Additional computational resources were provided in part by the National Science Foundation through major research instrumentation Grant No. CNS-09-58854 and major research instrumentation Grant No. CNS-1625061. R.C.R. acknowledges insightful discussions with S. Sastry, J. D. Weeks, and J. P. Perdew. We also thank A. Kohlmeyer for assistance with cp 2 k for the calculations in the SM [37] . Funding Information: Remsing Richard C. * Klein Michael L. Institute for Computational Molecular Science and Department of Chemistry, Temple University , Philadelphia, Pennsylvania 19122, USA Sun Jianwei Department of Physics and Engineering Physics, Tulane University , New Orleans, Louisiana 70118, USA * rremsing@temple.edu 26 April 2018 April 2018 97 14 140103 27 August 2017 31 January 2018 ©2018 American Physical Society 2018 American Physical Society Despite silicon being of great technological importance, an understanding of its behavior across the phase diagram is still lacking, especially near liquid-solid coexistence. The difficulty in describing silicon near coexistence from first principles lies in discriminating between the metallic and covalent bonds present in the material. Using the strongly constrained and appropriately normed (SCAN) density functional, which can describe a wide variety of bonds with quantitative accuracy, we report a thorough investigation of liquid silicon in the vicinity of liquid-solid coexistence using ab initio molecular dynamics simulations. We observe a structural transition in the supercooled regime that is rooted in a change in the electronic structure of the material. This transition is found to occur at a higher temperature than previous predictions. We also discuss implications of the observed change in interatomic interactions for empirical models of transitions between two distinct liquids. U.S. Department of Energy 10.13039/100000015 Basic Energy Sciences 10.13039/100006151 DE-SC0012575 Office of Science 10.13039/100006132 National Science Foundation 10.13039/100000001 CNS-09-58854 CNS-1625061 Silicon forms the basis for the semiconductor industry and is a material of great technological significance. Thus, understanding its physical properties across the phase diagram is of paramount importance to more efficiently control the synthesis of Si-based materials. However, silicon, like many other tetravalently bonded semiconductors, has a rich phase diagram and becomes a complex metallic liquid when melted. Moreover, the atoms in the liquid do not interact solely by metallic bonds, but a non-negligible fraction of covalent bonding persists, leading to a variety anomalous properties [1–14] . The delicate balance of covalent and metallic bonding in liquid silicon ( l - Si ) and related materials results in unique coordination structures that are suggested to underlie a phase transition in the supercooled region of the phase diagram. This postulated transition is between two distinct metastable liquid phases: a disordered, high-density liquid phase and a more ordered, low-density liquid phase [4,6–9,14–16] . However, describing this balance of interactions from first principles has proven quite challenging, and density functional theory (DFT)-based simulations of l - Si have not been able to quantitatively and sometimes even qualitatively predict essential structural and dynamic features [1,2,10–13,17–19] . In this Rapid Communication, we move toward a quantitative description of l - Si using the strongly constrained and appropriately normed (SCAN) meta-generalized gradient approximation (meta-GGA) [20] . SCAN has recently been shown to provide a quantitative representation of diversely bonded systems, including metallic, covalent, and even intermediate-range many-body van der Waals interactions [21–23] . Using SCAN, we provide a refined description of the behavior of l - Si at low temperatures, including an improved estimate of the melting point and the onset of features reminiscent of a liquid-liquid phase transition (LLPT). These refined estimates will likely aid future investigations of the properties of silicon near liquid-solid coexistence. We simulated l - Si at a range of temperatures using Born-Oppenheimer molecular dynamics (MD) simulations within the Vienna ab initio simulation package ( vasp ) [24] , following previous work [21,22] . Electronic structure calculations were performed using DFT within the framework of the projector augmented-wave method [25] , employing the SCAN meta-GGA [20] . Nuclear degrees of freedom are treated classically. The large mass of Si leads to nuclear quantum effects being negligible above temperatures of roughly 800 K [26] , unlike lighter nuclei, such as hydrogen, where such effects can be significant even at similar high temperatures [27] . Systems consisted of N = 216 Si atoms and were equilibrated for at least 15 ps before production runs of at least 22 ps were performed. All simulations were conducted in the isothermal-isobaric ensemble at zero pressure with a time step of 3 fs using a Parrinello-Rahman barostat [28] and Langevin thermostat [29,30] . Employing the N P T ensemble allows us to determine the equilibrium density and the structure and dynamics at that average density; we do not fix the system density at the experimental one. Previous work suggests that l - Si may undergo a LLPT in the deeply supercooled region of the phase diagram [4–9,14,15] , although the existence of a LLPT in some liquids remains a topic of debate [31–34] . In this scenario, one can cross the LLPT line by cooling at constant pressure, and a transition between a metallic liquid with high density and a semimetallic liquid with low density occurs [4,5] . Modeling the low-temperature behavior of l - Si with high accuracy is of great interest to describing these complex features and may aid the interpretation of similar behaviors in other liquids [35] . We first examine the temperature dependence of the bulk density, shown in Fig. 1 . Upon cooling, the predicted densities are in quantitative agreement with the experimental results [36] (solid line), before we observe a transition to low density between 1350 and 1300 K. After subsequent cooling to 1200 K, we heat the system until the reverse transition to high density is observed between 1400 and 1500 K, demonstrating hysteresis in the density on the simulated timescales; time series illustrating these transitions in density can be found in the Supplemental Material (SM) [37] . The transition between high and low density occurs at a significantly higher temperature than previously predicted by the local density approximation (LDA) [6] . However, the temperature difference between this transition and the melting point ( T m ) is comparable in the two models: T m LDA ≈ 1200 – 1300 K [13,19] while T m SCAN ≈ 1652 ± 46 K [37] . 10.1103/PhysRevB.97.140103.f1 1 FIG. 1. The bulk density ρ ( T ) displays hysteresis in the transition between high and low densities as the system is cooled from high T (cooling) or heated from low T (heating). The solid line corresponds to the parabolic fit obtained from experimental data [36] , while the dashed region is its extrapolation to low T . We emphasize that T m estimated by SCAN is in agreement with the experimental value of 1685 K and is a significant improvement over previous estimates within the LDA [13] and the Perdew-Burke-Ernzerhof (PBE) GGA ( 1540 ± 50 K) [38] . Upon moving from LDA to PBE to SCAN, the more covalent phase is stabilized, relative to the less covalent one, resulting in an increase in T m . This is consistent with a similar raising of the critical pressure of the semiconductor-metal transition at zero temperature from LDA to PBE to SCAN, but to a lesser extent, because the liquid retains some covalent character. [21] . The transition in the density is accompanied by significant structural changes. As the metallic liquid is cooled, we expect an increase in covalent character and structural signatures of this to appear. Indeed, the peaks in g ( r ) become more well defined as the liquid is cooled, undergoing a structural change between high and low densities, with an increase in the degree of order at low density [Fig. 2(a) ]. The coordination number N C ( T ) also undergoes a transition that closely follows the density, as shown in Fig. 2(b) . Here, N C ( T ) is defined as the number of particles within a sphere of radius r min around a central particle, and r min is defined as the location of the first minimum in g ( r ) . Four-coordinated tetrahedral coordination structures dominate at low T and ρ , while more disordered high-coordination structures dominate at high T and ρ . Note that N C predicted by SCAN at high densities is in agreement with experiments [39] , further supporting the accuracy of this functional, while the experimental N C values near the transition region are between those of the low- and high-density liquids, but still in reasonable agreement with the simulation predictions. 10.1103/PhysRevB.97.140103.f2 2 FIG. 2. (a) The transition in the density is accompanied by significant structural changes, highlighted by pair distribution functions g ( r ) at high and low temperatures. (b) The temperature dependence of the coordination number N C ( T ) closely follows that of the density. Experimental data (orange points) was obtained by x-ray diffraction [39] . Arrows indicate the direction of changing T . (c) Significant changes in liquid structure are further highlighted by the structure factor S ( k ) , in which the shoulder in the main peak at high temperatures becomes the dominant peak at low T . Also shown is S ( k ) / 5 for crystalline Si obtained from an MD simulation with the SCAN functional. (d) The fraction of nearest neighbors that are covalently bonded as a function of temperature η ( T ) illustrates that the structural changes are accompanied by a change in bonding. We also examine the temperature dependence of the structure factor, (1) S ( k ) = 1 N ∑ i = 1 N ∑ j = 1 N e − i k · r j − r i , where k is the wave vector and r i is the position vector of atom i . The structure factor of the high T liquid phase displays a major peak at k ≈ 2.75 Å − 1 and a shoulder at k ≈ 3.35 Å − 1 [see Fig. 2(c) ]. The balance of the two subpeaks that make up the first peak in S ( k ) is sensitive to temperature changes. At high T = 2200 K, a single broad peak appears at small k , that splits upon reducing the temperature to 1800 K. Reducing the temperature further increases the intensity of the subpeak at larger k relative to the peak at small k . At low densities, S ( k ) displays a significant separation of the first peak. To interpret the behavior of S ( k ) , we consider l - Si to be a mixture of two components, atoms bound by (i) metallic and (ii) covalent interactions. The latter leads to tetrahedral coordination environments, while metallic bonding can lead to disordered structures. The relative composition of these components varies with temperature, such that metallic bonds dominate at high T and covalent bonding dominates at low T . Previous work found that Si atoms with an interatomic distance less than r c ≈ 2.49 Å are covalently bonded [12] . We use this distance cutoff to quantify the number of covalently bonded nearest neighbors N B [37] . The fraction of covalently bonded nearest neighbors, η = N B / N C , as a function of temperature [Fig. 2(d) ] indicates that the high-density state is predominantly noncovalent and transitions to a low-density state dominated by covalent bonding at low T . Note that the behavior of η ( T ) across the transition region is a result of sharp transitions in both N C ( T ) and N B ( T ) [37] . The increase in η at low T results in tetrahedral coordination structures [1,11,12] , while the high-density state is significantly more disordered [37] . The splitting of the first peak of S ( k ) into two well-defined peaks at low T is a result of this increase in covalently bonded neighbors, and the resulting positions of these peaks are similar to those in the crystalline [Fig. 2(d) ] and amorphous [40] solid phases of silicon. The increase in covalent bonding when l - Si transitions from high to low density is suggested to manifest an accompanying metal-to-semimetal transition in the electronic structure [4] . This behavior is observed here in the temperature dependence of the electronic density of states N ( E ) , shown in Fig. 3(a) , which behaves in accord with previous findings [4] . The liquid undergoes a metal-to-semimetal transition when the density transitions to lower values, evidenced by the dramatic increase in the depth of the minimum in N ( E ) at the Fermi energy E F . 10.1103/PhysRevB.97.140103.f3 3 FIG. 3. (a) The density of states N ( E ) at 1800 and 1200 K indicate that l - Si undergoes a metal-to-semimetal transition as the temperature is reduced. (b) The order parameter N ̃ F ( T ) , the scaled density of states at the Fermi level, shows behavior qualitatively similar to the density, indicating that a transition in electronic structure accompanies the structural transition as the system is cooled. (c), (d) The transition in electronic structure can be visually observed through isosurfaces of the electron density (red) in configurations taken from simulations. The snapshots show that interatomic covalent bonds are mostly absent at (c) high T = 1800 K , while Si atoms (yellow) are nearly always connected by electron density at (d) low T = 1200 K . We quantify the transition in electronic structure through the order parameter N ̃ F ( T ) = N ( E F ; T ) / N ( E F ; T = 2200 K ) , which is equal to unity at the highest temperature studied here and zero when the system is semiconducting. The T dependence of N ̃ F ( T ) [Fig. 3(b) ] indicates that the electronic structure of the system changes in concert with its density, and that a metal-to-semimetal transition occurs when the system transitions from high to low density. This transition in interatomic interactions can be observed qualitatively by visualizing isosurfaces of the electron density in representative simulation configurations, shown in Figs. 3(c) and 3(d) . At high temperature, metallic bonding dominates and Si atoms are rarely connected by well-defined regions of electron density. In contrast, the tetrahedral coordination structures observed at low T are the result of covalent bonding, evidenced by well-defined regions of electron density between atoms [Fig. 3(d) ]. Note that we do not find evidence of a pseudogap in N ( E ) at the Fermi level in the high-density phase, in contrast to a previous report that used the empirical Stillinger-Weber (SW) potential [41] to predict the location of the liquid-liquid transition and the volumes of the low- and high-density phases as input to GGA-level calculations [5] . This may be attributed to both states studied in that work being below the LLPT and in the low-density state when described within the GGA, but not with the SW potential, as supported by the results presented here. From this analysis, it is clear that the nature of the interparticle interactions in l - Si changes significantly at a temperature between 1300 and 1500 K, manifesting a transition between high- and low-density states. Phenomenological descriptions of LLPTs and their associated features can be obtained through an empirical description of interparticle interactions that involve two length scales [42,43] . Such two-length-scale potentials can yield systems with a well-defined liquid-liquid critical point and associated phase transitions between high- and low-density liquids [42–44] . Previous work has even provided a quantitative mapping between a liquid and its coarse-grained, two-length-scale counterpart [45] . The use of temperature-independent pair potentials neglects the change in interparticle interactions with temperature observed in l - Si [4,5] , as discussed further below. We first illustrate that a coarse-grained, pairwise description of l - Si yields a two-length-scale potential by using the inverse Boltzmann inversion technique [46] to obtain effective pair potentials between Si atoms at two temperatures, T = 1800 and 1200 K , referred to as high-density and low-density states with effective potentials u H ( r ) and u L ( r ) , respectively, shown in Fig. 4(a) . These interaction potentials were obtained by matching the g ( r ) of the effective system to that from ab initio MD (AIMD) at the state point of interest in an iterative procedure. Both potentials display two length scales, defined by the minima near r = 2.5 and 4 Å; the former corresponds to the length scale r c below which covalent bonding is observed. The minima and maxima are more prominent in u L ( r ) , which has a significant energy barrier between the two minima. The difference in potentials may not be surprising, because they are effective representations of qualitatively different interatomic bonding environments; u H effectively describes metallic bonds, while u L describes the tetrahedral structures resulting from s p 3 hybridized covalent bonding. The distinction between the interactions underlying the effective potentials renders them both unable to reproduce the structural changes that occur in l - Si with temperature [Figs. 4(b) and 4(c) ]; the g ( r ) are not qualitatively reproduced at state points for which the potentials were not parametrized. This suggests that it may be informative to incorporate this behavior in phenomenological descriptions of LLPTs governed by changes in electronic structure, through, for example, more complex empirical potentials [41,47–49] or state-point-dependent coarse-graining procedures [50–53] . 10.1103/PhysRevB.97.140103.f4 4 FIG. 4. (a) Effective pair potentials u H ( r ) and u L ( r ) that reproduce the g ( r ) of l - Si at high density ( T = 1800 K) and low density ( T = 1200 K ) , respectively. (b), (c) The effective potentials each yield nearly identical g ( r ) for different temperatures and are unable to describe the structural changes observed upon cooling/heating the liquid. Colors correspond to those in (a). Solid lines indicate the state at which the potential was parametrized, while data points indicate the state point where a transition should be observed, but is not found, for these effective potentials. We have not addressed the dynamics surrounding the observed transition. It has been suggested that features ascribed to LLPTs in many simulations are due to a misinterpretation of behaviors that occur upon coarsening of a liquid phase on its way to becoming a solid [33,34] . The dynamics of this process are slow, and any claims regarding the existence or nonexistence of a LLPT requires advanced sampling techniques and a thorough characterization of the timescales underlying liquid-solid and possible LLPTs in the material of interest. While such an investigation is beyond the scope of this work, estimates of timescales from our simulations and a discussion of their importance can be found in the SM [37] . We note that convergence of averages in the supercooled regime is not guaranteed, because our simulations are only several relaxation times in length. Despite this fact, we have demonstrated that the SCAN functional improves the description of the low-temperature behavior of l - Si from ab initio simulations, evidenced by the improved estimate of the melting temperature T m = 1652 ± 46 K , in agreement with the experimental value of 1685 K. Additionally, structural changes that have been attributed to a low-temperature phase transition are found to occur between 1300 and 1500 K, in contrast to the 1060 K predicted previously with LDA. Our improved estimates suggest that experiments searching for proposed liquid-liquid transitions in silicon should focus on temperature ranges 150–350 K below the melting point of Si, which require significantly less supercooling than previous LDA-based estimates would suggest. Such temperatures can be accessed in experiments and are above estimates of the homogenous nucleation temperature of the crystalline solid phase [54] . Funding Information: This work was supported as part of the Center for the Computational Design of Functional Layered Materials, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award No. DE-SC0012575. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy. Additional computational resources were provided in part by the National Science Foundation through major research instrumentation Grant No. CNS-09-58854 and major research instrumentation Grant No. CNS-1625061. Publisher Copyright: © 2018 American Physical Society.
PY - 2018/4/26
Y1 - 2018/4/26
N2 - Despite silicon being of great technological importance, an understanding of its behavior across the phase diagram is still lacking, especially near liquid-solid coexistence. The difficulty in describing silicon near coexistence from first principles lies in discriminating between the metallic and covalent bonds present in the material. Using the strongly constrained and appropriately normed (SCAN) density functional, which can describe a wide variety of bonds with quantitative accuracy, we report a thorough investigation of liquid silicon in the vicinity of liquid-solid coexistence using ab initio molecular dynamics simulations. We observe a structural transition in the supercooled regime that is rooted in a change in the electronic structure of the material. This transition is found to occur at a higher temperature than previous predictions. We also discuss implications of the observed change in interatomic interactions for empirical models of transitions between two distinct liquids.
AB - Despite silicon being of great technological importance, an understanding of its behavior across the phase diagram is still lacking, especially near liquid-solid coexistence. The difficulty in describing silicon near coexistence from first principles lies in discriminating between the metallic and covalent bonds present in the material. Using the strongly constrained and appropriately normed (SCAN) density functional, which can describe a wide variety of bonds with quantitative accuracy, we report a thorough investigation of liquid silicon in the vicinity of liquid-solid coexistence using ab initio molecular dynamics simulations. We observe a structural transition in the supercooled regime that is rooted in a change in the electronic structure of the material. This transition is found to occur at a higher temperature than previous predictions. We also discuss implications of the observed change in interatomic interactions for empirical models of transitions between two distinct liquids.
UR - http://www.scopus.com/inward/record.url?scp=85046808473&partnerID=8YFLogxK
UR - http://www.scopus.com/inward/citedby.url?scp=85046808473&partnerID=8YFLogxK
U2 - 10.1103/PhysRevB.97.140103
DO - 10.1103/PhysRevB.97.140103
M3 - Article
AN - SCOPUS:85046808473
SN - 2469-9950
VL - 97
JO - Physical Review B
JF - Physical Review B
IS - 14
M1 - 140103
ER -