Yielding under compression and the polyamorphic transition in silicon
Jan Grie{\ss}er, Gianpietro Moras, Lars Pastewka

TL;DR
This study uses molecular simulations to reveal that amorphous silicon undergoes a discontinuous polyamorphic transition under high pressure, characterized by elastic instabilities and plastic events similar to shear yielding.
Contribution
It demonstrates that the polyamorphic transition in amorphous silicon is a nonequilibrium yield-like transition with discrete plastic events and elastic instabilities during compression.
Findings
Discontinuous density and elastic constant changes at 13-16 GPa.
Plastic events accompanied by shear modulus vanishing.
Transition differs from gradual structural change in quench simulations.
Abstract
We investigate the behavior of amorphous silicon under hydrostatic compression using molecular simulations. During compression, amorphous silicon undergoes a discontinuous nonequilibrium transition from a low-density to a high-density structure at a pressure of around -~GPa. Ensemble-averaged density and elastic constants change discontinuously across the transition. Densification of individual glassy samples occurs through a series of discrete plastic events, each of which is accompanied by a vanishing shear modulus. This is the signature of a series of elastic instabilities, similar to shear transformation zones observed during shear yielding of glasses. We compare the structure obtained during compression with a near-equilibrium form of amorphous silicon obtained by quenching a melt at constant pressure. This gives structures identical to nonequilibrium compression at low and…
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsHigh-pressure geophysics and materials · Diamond and Carbon-based Materials Research · Material Dynamics and Properties
Yielding under compression and the polyamorphic transition in silicon
Jan Grießer
Department of Microsystems Engineering, University of Freiburg, Georges-Köhler-Allee 103, 79110 Freiburg, Germany
Gianpietro Moras
Fraunhofer IWM, MicroTribology Center TC, Wöhlerstraße 11, 79108 Freiburg, Germany
Lars Pastewka
Department of Microsystems Engineering, University of Freiburg, Georges-Köhler-Allee 103, 79110 Freiburg, Germany
Cluster of Excellence livMatS, Freiburg Center for Interactive Materials and Bioinspired Technologies, University of Freiburg, 79110 Freiburg, Germany
Abstract
We investigate the behavior of amorphous silicon under hydrostatic compression using molecular simulations. During compression, amorphous silicon undergoes a discontinuous nonequilibrium transition from a low-density to a high-density structure at a pressure of around - GPa. Ensemble-averaged density and elastic constants change discontinuously across the transition. Densification of individual glassy samples occurs through a series of discrete plastic events, each of which is accompanied by a vanishing shear modulus. This is the signature of a series of elastic instabilities, similar to shear transformation zones observed during shear yielding of glasses. We compare the structure obtained during compression with a near-equilibrium form of amorphous silicon obtained by quenching a melt at constant pressure. This gives structures identical to nonequilibrium compression at low and high pressure, but the transition between them occurs gradually rather than discontinuously. Our observations indicate that the polyamorphic transition is of a nonequilibrium nature, and it has the characteristics of a yield transition that occurs under compression instead of shear.
††preprint: APS/123-QED
I Introduction
Materials can exist in structurally distinct forms in their crystalline state, a property which is called polymorphism [1]. Which form is actually present depends on the external conditions temperature and pressure and on the way the material is formed and processed afterwards. One form can transform into another, for example when increasing the pressure above some critical value. In crystals, pressure-induced phase transitions are commonly between equilibrium states of the system [2, 3, 4]. Amorphous solids, however, are intrinsically out-of-equilibrium, yet pressure-induced transformations between distinct amorphous forms are possible [5, 6, 7, 8, 9, 10].
Materials with multiple amorphous modifications are called polyamorphic [5, 6, 7, 8, 10]. One prominent example of a material featuring such a polyamorphic phase transition is silicon [11, 12, 13, 14, 15, 16, 17, 18, 19]. Silicon crystallizes into the diamond cubic phase under ambient conditions and transforms into the -Sn phase upon hydrostatic compression [20]. An indication for a polyamorphic transition in the corresponding amorphous form is that the melting temperature in silicon decreases with increasing pressure, (see Fig. 1a), meaning that the liquid phase is denser than the low-temperature diamond cubic crystal. A maximum in the melting curve, which for silicon is expected to occur at negative pressure [21, 7, 22], is typically explained by a “two-state” model [23]. This model assumes that low-density and high-density phases coexist in the liquid state, with the proportion of high-density domains increasing with pressure, thus leading to a pressure-driven transition from a low-density to a high-density liquid phase [13, 21, 7, 10]. A transition between two amorphous modifications is expected to reflect this liquid-liquid phase transition.
Experimental studies revealed that amorphous silicon (aSi) can indeed exist in two distinct structures. Upon compression, aSi transforms from a low-density structure (LDA) to a high-density structure (HDA) at a critical pressure between GPa and GPa [11, 15, 24, 25]. These experiments were supported by early numerical simulations, which showed a transition from an initially four-fold coordinated LDA structure to an HDA structure with five-fold coordinated atoms [14], or to a very high density amorphous (VHDA) structure with even higher coordination [12]. This structural transformation of aSi has been numerically investigated by multiple authors since, always showing a coexistence of LDA and HDA regions that sequentially transform to HDA (or VHDA) during compression [12, 14, 18, 19].
There are multiple well-documented mechanisms that lead to structural transitions between polymorph. Besides equilibrium phase transitions, crystalline materials have ultimate stability limits where the crystalline lattice collapses as a consequence of an elastic or dynamic instability [26, 27, 28, 29]. The result of such an instability is that all atoms in a perfect and defect-free crystal rearrange at once, forming a new stable crystal. Plasticity during shear loading of amorphous structures is associated with similar, yet localized mechanical instabilities [30, 31, 32, 33, 34, 35]. The rearranging regions are called shear-transformation zones and their occurrence can be predicted through the identification of “soft spots” [36, 37, 34, 35]. In aSi, these soft spots have been identified as local regions of HDA character [38, 39].
This raises the question whether plasticity in amorphous solids and the polyamorphic transition have different microscopic signatures, or conversely, whether this transition is nothing else than some sort of yield point.
In the theory of equilibrium phase transitions, the distinction between first and second order is made on the behavior of susceptibilities across the transition point [40, 41]. A partial answer to the question on the nature of the polyamorphic transition may therefore rest in how the (linear) mechanical response of aSi behaves across the transition. We here use atomic-scale computer simulations of periodic representative amorphous volume elements, to study the macroscopic (elastic) and microscopic (shear transformation) response of aSi across the polyamorphic transition. To drive the transition, we compress (and decompress) these volume elements affinely in small increments, and search for the next local minimum after the deformation step. We call this procedure athermal quasistatic compression (AQC), in analogy to the athermal quasistatic shear procedure used in the amorphous plasticity literature [42, 43]. It mimicks experiments where pressure is induced either through a pressure medium in a loading cell or by contacting with an indenter [11, 15, 24, 25, 16, 20].
Our results indicate, that the collapse of the aSi structure has signatures similar to the yield point during shear loading. The transition can be viewed as a sequential series of spatially localized compression transformation zones or plastic events. To emphasize the nonequilibrium character of the transition, we also probe the properties of a form of aSi that is compressed in the liquid phase and then quenched at constant pressure to become amorphous. The resulting amorphous samples show a gradual, rather than a discontinuous, variation of thermodynamic properties (density, elastic constants) as a function of pressure.
II Methods
We employ molecular dynamics and molecular statics simulations to investigate the behavior of aSi under hydrostatic compression. We first melt a Si-I (diamond cubic) crystal with atoms in the isothermal-isobaric (NPT) ensemble at K and zero hydrostatic pressure. The liquid is cooled down to K at a constant rate of K/s and equilibrated for further ps.
Our final aSi configurations are prepared from this silicon melt in two distinct ways. First, we produce periodic representative volume elements of aSi at zero hydrostatic pressure. We draw new random velocities from a Gaussian distribution at K and equilibrate the structures for ps. This equilibration time is sufficient for the atoms to loose memory of their initial state i.e. for the mean-squared displacement to be in the diffusive regime. After equilibration, we quench the liquid samples from K to K at quench rates of K/s and K/s. To drive the polyamorphic transition, the final glass is then compressed hydrostatically by affinely remapping all atomic positions with a prescribed density increment [42, 43]. We optimize the glass structure after each increment, a procedure we refer to as AQC in the following. Second, we equilibrate liquid silicon at K and a set of constant hydrostatic pressures GPa. The liquid is quenched to K using a quench rate of K/s while maintaining the constant hydrostatic pressure. We generated ensembles of independent configurations for the glasses quenched at vanishing pressure and configurations for each quench at constant pressure.
In all simulations, we used the interatomic potential by Kumagai et al. [44]. Previous publications showed that this interatomic potential correctly reproduces the polyamorphic transition in aSi [18] and delivers results on shear-induced silicon amorphization that are comparable to those obtained with more complex machine-learning potentials [45, 46]. All dynamic simulations employed a timestep fs, a Nose-Hoover chain thermostat and Parinello-Rahman barostat [47, 48, 49]. The relaxation constant of the thermostat and the barostat were chosen as ps and ps, respectively. Energy minimization in static simulations was performed using a conjugate-gradient minimizer with a force tolerance of eV/Å.
III Results
III.1 Phenomenology
We begin by investigating how amorphous silicon prepared at zero pressure behaves under hydrostatic compression. As shown in Fig. 1b, the density initially increases linearly with pressure. The material densifies at constant critical pressure , which manifests as an abrupt increase in the density-pressure diagram. This sudden compression is the signature of the pressure-driven transition from an LDA to an HDA structure in aSi. The critical pressure for this transition is at GPa and GPa for glasses prepared with quench rates of K/s and K/s, respectively. As the system is further compressed, the mean density of glasses prepared with different quench rates converges to the same value. The glasses prepared at a constant non-zero pressure show a smooth dependence on the hydrostatic pressure and therefore a continuous transition between the two phases. These results are consistent with prior work [18, 19].
Additional insight can be obtained by comparing the enthalpy of the two phases [3, 50, 12], , where is the potential energy, is the pressure and is the volume. In Fig. 1c, we show the enthalpy during compression and subsequent decompression. The two curves intersect at GPa, which is the pressure at which an “equilibrium” phase transition could occur. Glasses prepared at constant pressure follow the amorphous structure with the lowest enthalpy. We conclude that below GPa the LDA structure is thermodynamically preferred, while above GPa the HDA structure is more favorable. The critical pressures obtained from compression are more than twice this “equilibrium” pressure .
We now analyze how the atomic structure changes across this transition. Figures 2a and b show the probability of finding an atom with coordination .
For glasses prepared at a quench rate of K/s we show while we only show for the slower quench rate of K/s. During hydrostatic compression, before is reached, we observe only small changes in the mean coordination . In this regime, slowly quenched glasses have a larger fraction of four-fold coordinated atoms than samples prepared with a fast quench rate. Upon reaching the critical pressure, drops and an atomic configuration with mainly five-fold coordinated atoms (large ) emerges. These five-fold coordinated configurations are transient and vanish rapidly as atoms with even higher coordination numbers 6, 7, 8 and 9 emerge continuously for increasing hydrostatic pressure. The snapshots of compressed configurations in Fig. 2c (LDA) and d (HDA) also clearly show this increase in average coordination across the LDA-HDA transition. Note that some authors denote structures with high coordination numbers as VHDA [12, 14, 51, 19], but given the gradual transition towards highly coordinated structures we refer to the high pressure structures as HDA.
For glasses quenched at a constant pressure, the coordination number changes continuously with hydrostatic pressure. Thereby, low-coordinated atoms are continuously replaced by atoms with higher coordination as the pressure increases. At the pressure GPa, the configurations consist mainly of five- and six-fold coordinated atoms.
III.2 Elastic properties
The abrupt change in coordination during compression, as opposed to a continuous pressure dependency, strengthens the assumption of a nonequilibrium phase transition during compression. To further substantiate this hypothesis, we investigate the elastic properties across the polyamorphic transition. It is well known that in crystalline materials a mechanically driven transition occurs when an eigenvalue of the elastic tensor disappears [52, 53, 27]. In the limit of zero temperature and finite stress, the tensor of elastic moduli (also known as the Birch coefficients [54, 55, 28, 29]) is composed of three contributions: the affine or Born moduli , the non-affine moduli and a stress dependent contribution . The full elastic tensor is computed as [56, 57, 58, 29]
[TABLE]
with
[TABLE]
where is the energy function, is the current volume of the simulation cell, is the position of atom in direction , is the Green-Lagrange strain tensor, is the Cauchy stress, is the Kronecker delta and is the Hessian matrix. More details on this decomposition at finite stress and its analytical computation for many-body potentials can be found in Ref. [29].
Since amorphous solids show isotropic material behavior for sufficiently large system sizes and we consider a hydrostatic stress, the elastic tensor reduces to two independent elastic constants [56, 59, 60, 61]. We report the ensemble-averaged bulk modulus and shear modulus , since they are the eigenvalues of the elastic tensor. They determine the limit of elastic stability that is reached when either modulus vanishes [26, 27]. Only the shear modulus can vanish under hydrostatic compression.
Figure 3 shows the bulk modulus and shear modulus for the full range of pressures.
The bulk modulus of the glasses increases under compression for the full pressure range, independent of preparation. For aSi prepared at vanishing pressure, the bulk modulus initially increases linearly up to GPa and starts to become independent of pressure up to . During the transformation (at ), the bulk modulus increases discontinuously and subsequently approaches again a linear pressure dependency for higher compression. For glasses prepared at constant, nonzero pressure we observe a continuous, but sublinear, increase of the bulk modulus with increasing pressure. The pressure dependency of the bulk modulus resembles the pressure dependency of the density in Fig. 1b. The difference between the quenched and AQC glasses is purely due to their difference in density. The bulk modulus data collapses when plotted versus density rather than pressure.
The shear modulus shows more interesting features. During AQC, the shear modulus decreases to a value of about GPa at the critical pressure. Once the is reached, the modulus jumps to a value of roughly GPa, almost identical to the modulus for the stress-free quenched structure. While the initial quench rate leads to a different modulus at , the systems appear to loose memory of the initial state beyond the polyamorphic transition. At even larger compression, the shear modulus starts to soften again. The shear modulus is independent of preparation protocol above the critical pressure of the polyamorphic transition.
For glasses quenched at constant pressure, the shear modulus also has nonmonotonic functional dependency on . The maximum shear modulus is at GPa, which coincides with the equilibrium critical pressure . As the pressure increases above and , the shear modulus decreases and plateaus at GPa.
In Figure 4 we show the individual contributions from Eq. (1) to the bulk modulus and the shear modulus.
For AQC aSi, the pressure dependency of the Born contribution to the bulk modulus is qualitatively similar to what is observed for the full bulk modulus in Fig. 3a. While the Born contribution hardens the elastic response, the non-affine contribution softens it. With increasing compression, the magnitude of the non-affine contribution increases initially and reaches a local maximum at the critical pressure. At the critical pressure, the ratio between the Born and the non-affine contribution has a local minimum of . This minimum explains the flattening of the bulk modulus observed in Fig. 3a. Above the transformation, the non-affine contribution decreases its magnitude with pressure. The aSi samples prepared at a constant pressure show a pressure-dependency similar to the AQC glasses for both the Born and the non-affine contribution. However, their pressure dependency is continuous and no abrupt change of either property is observed. As expected from Eq. (4), the stress contribution increases linearly and is independent of the preparation protocol. The stress contribution additionally stiffens the material, but its absolute value is small compared to the Born contribution.
The Born contribution to the shear modulus increases linearly under hydrostatic compression for the aSi samples prepared at vanishing pressure. The behavior of the non-affine contribution is the reserve of the Born contribution, i.e. the magnitude of the non-affine contribution decreases linearly. At the critical pressure for the transition, the two contributions have a local maximum, or respectively minimum, with a ratio : The non-affine contribution almost compensation the Born term. At the transformation, the pressure dependency changes abruptly. For larger compression, the affine contribution deceases and the non-affine contribution increases. The stress-dependent part of the shear modulus softens the elastic response in addition to the non-affine contribution. This ultimately explains the observed softening of the shear modulus in Fig. 3 before the transformation. Samples prepared at constant pressure follow again the trends of the hydrostatically compressed glasses but without the abrupt change in behavior.
III.3 Elastic instabilities and soft spots
Although the (ensemble-averaged) shear modulus softens under compression, its value at the critical pressure for the phase transformation is still finite. The material does not appear to exhibit an elastic instability. While the ensemble-averaged data does not show an instability, the individual calculations undergo a continuous series of elastic instabilities, i.e. a divergence of the shear modulus in the pressure range of the polyamorphous transition. Figure 3c shows the non-affine displacement field resulting from one of these elastic instabilities. This non-affine displacement field reveals a random displacement of all atoms, but localized on one cluster of atoms. Figure 3d shows the drops in the shear modulus resulting from a sequence of such events. Comparing the shear modulus before and after an instability, we observe that the absolute value barely changes for pressures . At the critical pressure for the transition (which is slightly configuration dependent), the number of elastic instabilities per pressure range increases. Interestingly, each instability near lead to a noticeable stiffening of the shear modulus.
We investigate the origin of this sudden change by looking at the microscopic signature of these individual plastic events, which are reminiscent of shear transformation zones (STZs) but occur under compression rather than shear. Many methods exist to predict and identify soft spots in amorphous solids [35]. In this work we use the method developed by Richard et al. [62], which relies on the pseudo harmonic modes (PHMs). We choose PHMs because they predict the displacement field of a soft spot well before the instability [62]. We compute the lowest PHM for a subset of 50 configurations as a function of hydrostatic pressure, and report their approximate size using the participation ratio. The participation ratio is defined as
[TABLE]
where (which is a -dimensional vector) is the displacement field of the PHM and is the number of atoms in the solid. The participation ratio is a measure of spatial localization and is the approximate number of atoms participating in the plastic event. For the special case where only one atom is involved in a plastic event, it yields . In contrast, if all atoms are involved and they contribute equally, .
Figures 5a and b show the product of the PHM as a function of hydrostatic pressure.
The mean size of the PHMs increases slightly with compression before the phase transformation. Independent of quench rate, the mean size of the PHMs in this regime involves approximately atoms. At the transition pressure, the mean size of the soft spots starts to increase up to atoms. For larger compression, the number of atoms involved in a plastic event increases up to approximately atoms. The displacement field for one of the modes in the low and high pressure regime is shown in Fig. 5c and d. It is clearly visible that the size of the soft spots increases with compression and that the displacement field has an Eshelby-like character [63, 64, 35].
IV Discussion
First, we note that equilibrium arguments appear to be incapable of capturing the polyamorphic transition in aSi. In particular, the polyamorphic transition does not occur near the “equilibrium” pressure . However, does capture aspects of the transition from low- to high-density structures when compressing the melt before quenching – in particular the maximum in shear modulus (Fig. 3b). This indicates a large kinetic barrier resisting the structural transition, which emphasizes its nonequilibrium character.
Our observations on the polyamorphic transition have similarities to shear yielding of glasses. Shear-yielding of aSi has been reported in the literature, with the a phenomenology as follows: aSi initially deforms elastically but flows after the yield point [39, 65, 66, 64, 67]; approaching the yield point, a large percentage of atoms become -fold coordinated [39, 65, 66]; the shear modulus changes abruptly at the yield point [64]. We observed the same phenomena for the polyamorphic transition, with the difference that atoms with coordination larger than 5 appear during compression.
The microscopic picture of the polyamorphic transition found here is also identical to yielding under shear in glasses. Hydrostatic compression of aSi results in series of instabilities, which can be identified by loss of mechanical stability that manifests as a vanishing shear modulus. Each instability triggers a “soft” spot in the material: The displacement field during these instabilities shows signatures identical to STZs [59, 32, 35], the carrier of plastic deformation in disordered materials under shear. These signatures include a localization of the non-affine displacement field [36] as well as a quadrupolar strain-field surrounding the localized event, as described by Eshelby inclusion theory [63]. The fast-quenched configurations initially have a larger fraction of -fold and higher-coordinated atoms than the configurations obtained at lower quench rates. These are “liquid-like”, i.e. less resistant to deformation and acting as soft spots in aSi [38, 39]. This lower resistance to deformation reduces the critical pressure.
The macroscopic shear modulus, computed here as ensemble averages over configurations, does not vanish. However, it monotonously decreases with pressure below of the polyamorphic transition and jumps to a larger value as the transition is crossed. This is accompanied by an instantaneous change in density, marked by an increase in average coordination. The decrease in shear modulus is driven by the non-affine contribution to the elastic modulus (Fig. 4e). In terms of the microscopic picture painted in the previous paragraph, the non-affine origin of the decrease in shear modulus can be interpreted as a signature of localized softening of the material.
To further emphasize the similarities to yielding, we plot in Fig. 6a the pressure as a function of volumetric strain in a traditional stress-strain diagram. The behavior is qualitatively identical to stress-strain diagrams obtained from shear yielding glasses [43, 38, 39, 65, 66, 68], with pressure taking the role of shear stress. For visual comparison, we also computed the shear response of our aSi samples and show in Fig. 6b the ensemble-averaged shear stress as a function of shear strain . The material exhibits a linear “elastic” response at small strain followed by flow at constant pressure beyond a yield point at und both compression and shear. In particular, the glass prepared at a lower quench rates shows a stress overshoot, explaining the reentrant section of the density-pressure diagram (Fig. 1b) as a compression-softening phenomenon of a well-aged glass [69, 39, 68].
V Conclusions
In summary, we performed extensive simulations of the behavior of aSi during hydrostatic compression. Our simulations show a polyamorphic transition from a low-density to a high-density phase during quasistatic compression of aSi quenched at zero pressure, but a gradual transition for near-equilibrium structures obtained when quenching liquid silicon at finite pressure. Microscopic and macroscopic analysis of the transition revealed similarities to the yield transition of glasses under shear: Localized carriers of plastic events (shear transformation zones), change of local atomic order, and change in elastic properties. We conclude that the polyamorphic transition is essentially a yield transition. Whether yield of glasses itself is a phase transition is an open discussion in the scientific literature [70, 71, 72, 68, 73].
Acknowledgment
We thank Richard Leute and Thomas Reichenbach for fruitful discussions. Molecular dynamics simulations were carried out with lammps [74, 75]. We used ase [76] and matscipy for analysis of results and computation of the elastic properties. Atomic configurations were visualized with ovito [77]. The authors acknowledge support from the Deutsche Forschungsgemeinschaft (DFG, grant PA 2023/2). Simulations were carried out on NEMO at the University of Freiburg (DFG grant INST 39/963-1 FUGG).
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Callister and Rethwisch [2007] W. D. Callister and D. G. Rethwisch, Materials science and engineering , Vol. 7 (Wiley New York, 2007).
- 2Stanley [1971] H. E. Stanley, Phase transitions and critical phenomena , Vol. 7 (Clarendon Press, Oxford, 1971).
- 3Kubo et al. [1998] R. Kubo, M. Toda, and N. Saito, Statistical physics I: equilibrium statistical mechanics (Springer Science & Business Media, 1998).
- 4Schwabl [2006 a] F. Schwabl, Statistical mechanics (Springer Science & Business Media, 2006).
- 5Poole et al. [1997] P. H. Poole, T. Grande, C. A. Angell, and P. F. Mc Millan, Polymorphic phase transitions in liquids and glasses, Science 275 , 322 (1997) . · doi ↗
- 6Brazhkin and Lyapin [2003] V. V. Brazhkin and A. G. Lyapin, High-pressure phase transformations in liquids and amorphous solids, J. Phys.: Condens. Matter 15 , 6059 (2003) . · doi ↗
- 7Wilding et al. [2006] M. C. Wilding, M. Wilson, and P. F. Mc Millan, Structural studies and polymorphism in amorphous solids and liquids at high pressure, Chem. Soc. Rev. 35 , 964 (2006) . · doi ↗
- 8Mc Millan et al. [2007] P. F. Mc Millan, M. Wilson, M. C. Wilding, D. Daisenberger, M. Mezouar, and G. N. Greaves, Polyamorphism and liquid–liquid phase transitions: challenges for experiment and theory, J. Phys.: Condens. Matter 19 , 415101 (2007) . · doi ↗
