Chemical bonding in chalcogenides: the concept of multi-centre hyperbonding

A precise understanding of the chemical bonding in amorphous, and crystalline, chalcogenides is still lacking due to the complexity arising from the delocalization of bonding, and non-bonding, electrons. Although an increasing degree of electron delocalization for elements down a column of the periodic table is widely recognized, its influence on chemical-bonding interactions, and on consequent material properties, of chalcogenides has not been comprehensively understood. Here, we provide a chemical-bonding framework for understanding chalcogenides by studying prototypical telluride non-volatile-memory, 'phase-change' materials (PCMs), and related chalcogenide compounds, via density-functional-theory molecular-dynamics (DFT-MD) simulations. Identification of the presence of multi-centre 'hyperbonding' interactions elucidates not only the origin of various material properties and their contrast between amorphous and crystalline phases, but also the very similar chemical-bonding nature between crystalline and amorphous PCMs, in marked contrast to the existing viewpoint. This new perspective will help in designing chalcogenide materials for diverse applications from a fundamental chemical-bonding point of view.


Introduction
Chalcogenide materials exhibit interesting properties for a wide range of applications, ranging from optical [1] or optoelectronic [2] applications, and encompassing applications of topological insulators [3,4], low-dimensional materials for electronics [5,6] or thermoelectricgeneration devices [7]. This broad capability is feasible due to the wide tunability of material properties, depending on chalcogen-atom types, or by their broad compositional glass-forming ability [8].
Of all the chalcogenides, some tellurides, viz. 'phase-change' materials (PCMs), show unique material properties, i.e. ultrafast crystallization rates and large (opto-electronic) property contrasts between amorphous (a-) and crystalline (c-) phases [9][10][11]. The simultaneous occurrence of both these material properties is counterintuitive: the former is supposed to involve small structural differences between a-and c-phases, while the opposite trend would be expected for the latter. Consideration of electron delocalization is also seemingly necessary, in that a single Lewis structure, with bonding and non-bonding electron pairs, is unable fully to describe the electron distributions in both a-and c-PCM phases [12,13]. In order to answer this question, chemical-bonding models, assuming a drastic change in the nature of chemical bonding between the two phases, have been proposed [14,15], and form the current mainstream point of view.
It is useful to classify existing models in terms of the length scale regarding electron delocalization. The resonant-bonding model [12,14], based on a limited number of Lewisconforming resonant structures, may be suitable for describing simple crystals with minimum disorder: naturally, structurally-disordered a-PCMs are beyond the validity of this model. Likewise, crystals with significant amounts of structural disorder are also difficult to incorporate within this model: PCM types include metastable rocksalt-type Ge 2 Sb 2 Te 5 (c-GST) [16] (or more generally compositions along the pseudo-binary tie-line of GeTe-Sb 2 Te 3 [10]), Ag-In-Sb-Te (AIST) [17], I-V-VI 2 -type compounds [18], or c-GeTe with vacancies [19].
Although the relevant spatial scale was ill-defined, a measure of electron delocalization with substantial metallic character has been proposed, and termed 'metavalent bonding' [20]. On the other hand, a molecular-orbital (MO) approach, involving notably three-centre, fourelectron (3c/4e) interactions, or the valence-bond (VB) theory of hyperbonding, focus on linear triatomic bonding geometries, and have been successfully adopted to describe hypervalent molecules [21][22][23]. Formation of linear triatomic bonding geometries is also recognized in a-GST [13] and c-GeTe-Sb 2 Te 3 alloys [24], and a signature of electron delocalization within this bonding configuration is found in a-GST [13]. The advantage of this model over other models is that it can be characteristically linked to microscopic properties of the materials.
Here, we show that the concept of multi-centre hyperbonding interactions completes the picture of chemical bonding in chalcogenides. Chemical-bonding interactions for a-/c-PCMs, and associated material properties, are comprehensively elucidated with this model. It is shown that, without invoking any new type of bonding interaction, property contrasts between a-and cphases of PCMs can be naturally understood in terms of their different extent of hyperbonding. This is in dramatic contrast with the existing theory for PCMs. Figure 1 shows chemical-bonding-indicator data, measured at the bond critical point (BCP), in simulated models of a-/c-GST. The trends for a-GST reveal that the charge density ( " ) [25], electron-localization function (ELF B ) [26], negative of the local energy density (LED) [27], or negative of the integrated crystal-orbital Hamilton population (ICOHP) [28] all increase monotonically with decreasing interatomic distance; this trend may be readily understood, given that the bonding is covalent-like. This perspective is in accord with the observation that covalent-bonding interactions in molecules exhibit negative energy densities (while positive values are found for ionic or intermolecular interactions) [27]. Other results, however, seemingly render this conclusion vague: for instance, charge densities at the BCP (~0.05 / ' ( ) are rather small to be classified as being associated with purely covalent bonding [29]; and the values of ∇ * " are close to zero with a change of sign, an indication of the bonding character being intermediate between covalent (negative ∇ * " ) and ionic (positive ∇ * " ) [25]. This presumably results from the additional contribution of metallicity, i.e. a tendency for delocalization of valence electrons.

Results and Discussion
Interestingly, a similarly broad distribution of interatomic distances is also found for metastable rocksalt-type c-GST models (Fig. 1), as for a-GST, due to structural disorder [11,30]. Nearly all the bonding-character data-points for c-GST exactly overlap with the wider range of datapoints for a-GST (Fig. 1a-e). This surprising similarity indicates that individual interatomic interactions in c-GST are indistinguishable from those in a-GST for the same bond length; alternatively, chemical-bonding interactions in c-GST belong to a subgroup of the broad spectrum of interactions existing in a-GST. This finding is in stark contrast to the current consensus that the nature of bonding in a-and c-GST is dissimilar [11,14], which has been the basis for rationalizing the property contrasts in GST. Nevertheless, a meaningful difference exists between the peak positions of the interatomic-distance distributions: the peak positions for c-GST lie at longer interatomic distances compared to those for a-GST (Fig. 1f). Such elongated interatomic distances are reminiscent of the relatively elongated axial bonds recently found for defective-octahedral configurations in a-GST [13]. The comparison in Fig. 1f indeed reveals that the bond-length distribution for c-GST resembles the distribution for axial bonds [13], which is, along with the bonding-character data shown in Fig. 1a-e, indicative of a similar bonding nature between bonds in c-GST and the axial bonds in a-GST.
The term 'hyperbond' was coined [32] to emphasize the distinctive characteristics of the 3c/4e bonds relative to those of ordinary 2c/2e covalent bonds (see Methods for a detailed definition The proposed interaction model (Fig. 2b) may be validated as follows. First, bond lengths (or ELF b values) of 2c/2e bonds (B-C in Fig. 2(i)) increase (or decrease) with the formation of 3c/4e interactions (Fig. 2c), which conforms to the hyperbonding concept, as the antibonding level of the B-C bond becomes occupied via this interaction (Fig. 2b). Second, the dependence of the hyperbonding tendency on the band gap of materials (assumed to scale with ∆ /0"1 ) shows the expected trend ( Fig. 2d), i.e. a stronger hyperbonding tendency for materials with smaller band gaps. Lastly, the involvement of Te LPs in forming hyperbonds is directly supported by the observation (Fig. S5) that 96% (or 93%) of hyperbonds centred at Ge (or Sb) atoms have at least one, or two, three-fold coordinated Te atoms as ligands, each of which can donate the LP electrons required for the 3c/4e interactions (Figs. 2a, 2b). These percentages are much higher than the percentage of three-fold Te atoms (~50% [13]) found in these models.
In c-GST, a significant difference with a-GST is that the crystalline symmetry requires Ge and Sb atoms to reside in (distorted) octahedral-ligand environments, allowing for the formation of a maximum of three perpendicular hyperbond pairs, i.e. six bonds. This means that the crystalline structure of metastable c-GST inherently provides conditions favourable for strong hyperbonding, with near-linear alignments of p orbitals. The sharp increase in the percentage of hyperbonds for c-GST, compared to a-GST (Fig. 3a), is hence due to this crystal-structure amplification effect. This finding, together with the high polar covalency of hyperbonds, is manifested in a higher degree of ionicity in c-GST than in a-GST due to the charge transfer from Ge and Sb atoms to Te atoms ( Fig. S1) with the very pronounced increase of hyperbonds in c-GST. This reveals an interesting character of hyperbonds that, although their electronegativity difference remains the same, simply switching between a bond and a hyperbond changes the ionicity of the same bond (e.g. the B-C bond in Fig. 2a).
It is now clear that the bonding difference between a-and c-GST arises from the relative abundances of the bonding types; 2c/2e covalent bonding predominates in a-GST, with a minor We can link the hyperbonding concept to many material properties. The first to be considered is a salient feature of PCMs, namely the large optical-property contrast between a-and c-phases, e.g. of GST [14]. The relationship between hyperbonds and Born effective charges (BECs) is the key finding that associates these seemingly unrelated properties (Figs. 3a,b). The BEC, averaged over Ge, Sb, and Te atoms in a-GST (c-GST) models was found to be 2.72±0.98 (6.40±2.30), 3.72±1.76 (8.74±3.60), and -2.57±1.30 (-6.16±3.41), respectively. These values are comparable to those given in a previous report [33]. The averaged BEC values for the atomic constituents in a-GST are slightly higher than their nominal ionic charges, but the atomic constituents of c-GST exhibit much higher BEC values. The sharp increase of BEC values from a-to c-GST follows the similar trend of increasing ratios of proportion of hyperbonds to that of ordinary covalent bonds (Fig. 3a), indicative of an intimate connection between BEC values and the presence of hyperbonds. To obtain a deeper insight into their relationship, we decomposed the whole distribution of BECs into contributions from atoms belonging to groups of atoms involved in forming a specific number of hyperbonds. As the AOs of each atomic constituent are involved in forming more hyperbond pairs, the atomic constituent tends to show a higher BEC value (Fig. 3b). In other words, hyperbonding facilitates a larger BEC than does ordinary covalent bonding; the high sensitivity of electronic polarization to atomic displacements [34] therefore constitutes another characteristic of hyperbonding. Considering the reported proportional relation between BECs and dielectric constants [35], it is the drastically increased number of hyperbonds in c-PCMs that escalates BEC and dielectric-constant values. Thus, crystal symmetry, permitting hyperbonding, is an essential requirement for large optical contrasts between a-and c-phases: the resulting linear atomic alignment hence plays a crucial role. The importance of aligned p orbitals in giving high dielectric constants was also discussed by Huang et al. [36], supporting the present conclusion.
The formation, or decomposition, process of hyperbonds, depicted in Fig. 2a, provides an additional insight into the dynamical properties of PCMs. The first example of interest is the collective bond-breaking behaviour in PCMs found from laser-assisted-evaporation experiments [37], presented as supporting evidence for the existence of metavalent bonding [20]. The essential observation identified by us is that the reported multiple-event probability [37] appears to scale with the hyperbond content. For hyperbonding-dominant materials, external perturbations causing bond-breaking processes ((iii) → (i) in Fig. 2a) would generate -bonded dimers or structural motifs [13], embedded in a sea of relatively weaker hyperbonds at the surface: the probability of atoms being evaporated in the form of strongly bonded molecules (referred to as 'multiple events' in [37]) is likely to be high, compared to the case without the presence of such a bond-strength hierarchy. According to this scenario, collective bond-breaking can therefore be attributed to the unique response of strongly resonating hyperbonds.
Another finding is that, unlike in crystalline Ge-Sb-S (c-GSS) or Ge-Sb-Se (c-GSSe), interatomic distances in c-GST show a continuous distribution without any significant gap (or dip) that, otherwise, separates contributions from strong and weak interatomic interactions ( Fig.   1f and Fig. S2). This is because hyperbonding bridges the gap between strong 2c/2e covalent bonding and weak LP interactions. We postulate that the presence of hyperbonding can provide low-energy-barrier routes for bond formation and breakage: for instance, the transition from a  [13,15]. Finally, a hint of a correlation with phonon anharmonicity is given by the recent observation by Lee et al. [38] that a strong phonon anharmonicity has its fundamental root in the cubic (rocksalt) crystalline structure, enabling a linear alignment of p AOs. The same condition for maximizing hyperbonding suggests that strong anharmonicity may therefore accompany pronounced hyperbonding. It is also interesting to note that strong phonon anharmonicity is often linked to the presence of shallow double-well potentials due to weak Peierls distortions [30,39], which corresponds to the presence of (incomplete) hyperbonds. Now let us compare with other chalcogenides containing different types of chalcogens, as already partly discussed for Fig. 2d. On going from GSS to GSSe, and to GST, rather general trends in bonding characteristics are apparent: i) the difference between GSS and GSSe is relatively small, yet a drastic change occurs from GSSe to GST; and ii) the contrast between a-and c-phases becomes amplified. For instance, a rather small increase of the hyperbond ratio ( Fig. 3a) occurs from GSS to GSSe, followed by a much larger increment from GSSe to GST, together with the amplifying contrasts between a-and c-models. Although the contrast becomes weaker, basically the same trends are observed for the averaged BECs (Fig. 3a). In principle, the origin of this distinction between chalcogenides may be described in terms of the factors affecting hyperbonding (Fig. 2). However, we also found that the extent of sp hybridization of the constituent AOs conveniently discriminates these differences. For the amorphous phases, the trend of sp hybridization, as measured by the distance from the central atom to its LP MLWF centre (WFC LP ), is clearly in inverse proportion to the degree of hyperbonding (Fig. 3a). The same conclusion is reached for crystalline models (see Supplementary Information for details). This sp hybridization, stimulating stereochemical bonding activity of LPs [40,41], leads to structural distortion, followed by the breakage of linear atomic alignments in c-GSS, and less so in c-GSSe, i.e. the loss of hyperbonding. The local atomic-coordination change involved in the distortion in the crystalline models is found to involve transitions from the initial octahedral coordination to lower coordinated polyhedral units, notably the trigonal-pyramidal geometry [13], that can be described by the sequential decomposition processes of (iii) → (ii) → (i) in Fig. 2a, along the three hyperbonding directions. This is not the case for c-GST (see Fig. 3). The stability of multi-centre hyperbonding interactions, along with the sustainability of the cubic crystalline form, is, therefore, in increasing order of GSS, GSSe, and GST. The close correlation between the hybridization of AOs, multi-centre hyperbonding, and crystal structure rationalizes why most PCMs are tellurides rather than sulphides or selenides. With the same reasoning, it can also be understood why sulphides, or selenides, of binary crystalline IV-VI compounds adopt the blackphosphorus crystalline structure, which is more distorted from its hypothetical cubic rocksalt structure than are tellurides with a rhombohedral structure (e.g. GeTe or SnTe) [40].
A distinction between chalcogenides is also found in their extent of electron delocalization (Fig.   3d). Electron delocalization occurs to a significant extent in a-GST, best illustrated from the substantial increase of the spread (Ω) of bonding MLWFs with the heavier Gp. VI elements (Fig. 3d), an indication of increasing metallicity [40]. Such a metallic contribution explains the unusually small values of chemical-bonding indicators for a-and c-GST in Fig. 1. The additional contribution of hyperbonds to delocalization (see Fig. S6 by the coupled shortening, and lengthening, respectively, of long, and short, bonds (Fig. 4). As a result, stronger linkages across Ge-Te long-bond layers occurs via the process of (ii) → (iii) in Fig. 2a. Shortening of a small number of Ge-Te short bonds is also found. Although the overall changes in terms of atomic positions is marginal [46], the detailed charge redistribution is significant, as indicated by the large disorder in the ELF distribution induced by the introduced vacancies. This disorder is induced as a way of minimizing the system energy following 3c/4e interactions. Therefore, multi-centre hyperbonding interactions provide a useful means of stabilizing defects in c-GeTe, and presumably also in GST.
We also investigated chemical-bonding indicators for I-V-VI 2 -type compounds and obtained basically the same conclusions as illustrated above.

Conclusions
The concept of multi-centre hyperbonding extends conventional chemical-bonding theory to give a correct understanding of the electronic structure and properties of chalcogenides. The combination of hyperbonding and ordinary covalent bonding can elucidate the origin of chalcogen-dependent structural differences and material properties associated with heavy Group VI elements. The established connection between chemical bonding, crystal structure, and material properties provides a completely new perspective in understanding, and hence designing, chalcogenide materials for various applications, including phase-change-memory or thermoelectric-generation materials.

Ab initio DFT and molecular-dynamics (AIMD) simulations
We employed constant-volume ab initio simulations based on density-functional theory (DFT) with the generalized-gradient approximation (GGA), as implemented in the VASP code [47].
The projector augmented-wave (PAW) method [48] with the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional [49], or SCAN functional [43], were used. The HSE06 hybrid functional [50] was used to calculate the band gap of models. The outer s and p electrons of constituent atoms were treated as valence electrons. The plane-wave energy cutoff was 300 eV or 400 eV. The Brillouin zone was sampled at the Γ point for amorphous models, while a Monkhorst-Pack 2x2x2 grid was used for crystalline models. For AIMD simulations, the temperature was controlled by a Nosé thermostat algorithm, and the MD time step was 3 fs.
The amorphous and crystalline models were simulated in cubic supercells with periodic boundary conditions. Structural relaxation was performed using a conjugate-gradient method, until the force on any atom was below 0.01 eV/Å. The Born effective charges (BECs) were determined using density-functional perturbation theory, as implemented in the VASP code [47].

Amorphous and crystalline models
The amorphous models (315 atoms) of GSS, GSSe and GST were generated via the conventional 'melt-quench' method. Initial random configurations of atoms were mixed at 2000 K for a few ps, and then kept at liquid temperature for tens of ps. For each type of chalcogenide, three amorphous models were generated by quenching different liquid configurations, sampled at different MD trajectories with at least a 10 ps interval to ensure different initial quenching configurations. The quenching rate was 15 K/ps. The three different amorphous models were then used to collect statistics on static or dynamic material properties for each chalcogenide system. The density used for a-GST models was the experimental density of 5.88 g/cm 3 . The densities of a-GSS and a-GSSe were determined via two steps. First, an amorphous model for each type of chalcogenide was generated with an approximate density, estimated with reference to the density of a-GST, and then the cell volume was relaxed to find the minimum-energy density. This density was then used to regenerate amorphous models by repeating the first step. The pressure of models generated in this way was close to zero.
Similarly, for statistics, three rocksalt-type crystalline models (196 atoms) for each crystalline chalcogenide were generated. The anion sites were filled with relevant chalcogen atoms. The cation sites were randomly filled with Ge, Sb, and 20% of vacancies. For a proper comparison, the distributions of cations and vacancies were the same for all models. The initial crystalline structures with approximate densities correspond to the perfect rocksalt-like structures without any distortion. However, structural relaxation towards a minimum-energy configuration leads to a distortion of the local coordination, with the extent of the distortion strongly depending on the chalcogen type, as noted in the main text. These relaxed models were then analyzed to get an insight into the impact of sp hybridization of AOs on the hyperbonding tendency and on the crystal structures.

Chemical-bonding indicators and definition of bonds
Five different chemical-bonding indicators were analysed to characterize interatomic interactions. The electron charge density, its Laplacian, local energy density, and ELF were analysed at bond-critical points (BCPs) for the relevant atom pairs within the framework of Quantum Theory of Atoms in Molecules (QTAIM) [25]. The BCP for each pair of atoms was determined by the critic2 code [51]. The orbital-based indicator of ICOHP was determined by the LOBSTER code [52]. Maximally-localized Wannier functions and their spread were computed by the wannier90 code [53]. We analysed ELF attractors [54], representing LPs, which were then used to characterize the magnitude of sp hybridization of LPs for the crystalline models. We employed several different exchange-correlation functionals, namely PBE [49] and SCAN [43] for GST and GeTe models, and additionally HSE06-D3 [44] and rev-DFT-D2 [45] for GeTe models. Except for PBE, all the other functionals include the effect of van der Waals interactions. The same conclusions were reached for all functionals considered.
The term 'hyperbond' was used specifically to represent each bond of a bond pair constituting a linear triatomic bonding geometry with near-identical bond distances, while the term '3c/4e bond' was used in a more relaxed manner, as related to the bond angles and distances. In this study, we define a hyperbond when ELF b for both linear bonds exceeds a value of 0.5 [13], but other chemical-bonding indicators can be equivalently employed, due to the definite relationships among those indicators (Fig. 1). Similarly, an ordinary covalent bond was defined when the corresponding ELF b value exceeds 0.5. In this study, we focused on Ge-Te and Sb-Te 'correct' bonds, rather than 'wrong' (homopolar and Ge-Sb) bonds that are not observed in ideal crystals. As their concentrations, related to linear triatomic bonding geometries, are relatively very low, they were not considered further here.
The coordination of Ge and Te atoms in rhombohedral c-GeTe consists of three strong and three weak Ge-Te bonds. Strictly speaking, these long (weak) bonds in c-GeTe do not satisfy the bonding criteria adopted in this study (Fig. 4), but for simplicity, they will be called 'long bonds' in the main text. Models, generated with diverse exchange-correlation functionals, reproduced the experimental short and long bond lengths to within an error of ~±2%. The rhombohedral symmetry persists with pressure (Fig. 4), but the length ratio of short to long bonds continuously diminishes. In particular, the reduction of long-bond lengths is more pronounced than for short bonds.

Comparison between a-and c-GST models: Bader charge
In the main text, the interatomic interactions were characterized in terms of various chemicalbonding indicators [1][2][3][4] in order to emphasize the similarity of chemical-bonding interactions between amorphous and crystalline GST models for the same bond distances. Also, from the intrinsic properties of hyperbonds [5], the contrast in various microscopic material properties between amorphous and crystalline phases was elucidated: one of these properties is the difference in Bader charges. The formation of hyperbonds from the interaction between a LP and an antibonding state of a nearby bond (Fig. 2a)

1) Coordination change with relaxation
As noted in the main text, the degree of sp hybridization is a useful quantity for describing the structure of crystalline chalcogenides. For this, we compared structural-relaxation behaviour between crystalline models with the same initial rocksalt-like configurations. The local coordination change involved with the relaxation is mostly a transition from an initially octahedral coordination to the trigonal-pyramidal, seesaw, or square-pyramidal geometries, with the formation of a LP on each cation atom [5]. These transitions can be described by the sequential decomposition processes of (iii) → (ii) → (i) in Fig. 2a,

3) ELF attractors for LPs
The sp hybridization of AOs in c-GST was investigated by analyzing the distance from the central atom to the ELF attractors representing LPs (Fig. S3) [6].

The mechanism of 3c/4e hyperbonding
The mechanism of hyperbonding is described in Fig. 2. Figure S5 supports the proposed hyperbonding mechanism by showing that the ligand Te atoms are mostly the Te(3,1) type [5], that is, three-fold coordinated Te atoms with a single LP. The additional LP, compared to twofold coordinated Te atoms, is involved in 3c/4e hyperbonding interactions, as illustrated in the main text.

The spread of bonding MLWFs
The spreads of bond MLWFs (Ωs) for different chalcogenide models are shown in Fig. 3d. The larger standard deviations of Ωs for a-GST than a-GSS or a-GSSe models are attributed in the main text to the stronger delocalization of electrons involved in 3c/4e hyperbonding than that of electrons involved in ordinary 2c/2e covalent bonding, which is supported by the distribution of Ωs, as shown in Fig. S6.