Molecular Simulation of Covalent Bond Dynamics in Liquid Silicon
Richard C. Remsing Michael L. Klein
physics.chem-ph cond-mat.mtrl-sci cond-mat.stat-mech























[]rick.remsing@rutgers.edu
Department of Chemistry and Chemical Biology, Rutgers University, Piscataway, NJ 08854

[]mike.klein@temple.edu
Institute for Computational Molecular Science and Department of Chemistry, Temple University, Philadelphia, PA 19122

















Many atomic liquids can form transient covalent bonds reminiscent of those in the corresponding solid states.

These directional interactions dictate many important properties of the liquid state,
necessitating a quantitative, atomic-scale understanding of bonding in these complex systems.

A prototypical example is liquid silicon, wherein transient covalent bonds give rise to local
tetrahedral order and consequent non-trivial effects on liquid state thermodynamics and dynamics.

To further understand covalent bonding in liquid silicon, and similar liquids,
we present an ab initio simulation-based approach for quantifying
the structure and dynamics of covalent bonds in condensed phases.

Through the examination of structural correlations among silicon nuclei and maximally localized
Wannier function centers, we develop a geometric criterion for covalent bonds in liquid Si. 

We use this to monitor the dynamics of transient covalent bonding in the liquid state and estimate
a covalent bond lifetime. 

We compare covalent bond dynamics to other processes in liquid Si
and similar liquids and suggest experiments to measure the covalent bond lifetime. 









Molecular Simulation of Covalent Bond Dynamics in Liquid Silicon
    Michael L. Klein
    December 30, 2023
================================================================











§ INTRODUCTION




Many metallic and semi-metallic atomic liquids contain significant numbers of dynamic covalent bonds
reminiscent of the static bonds formed in the corresponding solid state.

The covalent bonds in these liquids are dynamic, readily breaking and forming
on molecular timescales, and the characterization of these processes is complicated 
due to the interplay between electronic and nuclear structures.

Such metallic liquids include molten silicon <cit.>, boron <cit.>,
gallium <cit.>, and hydrogen at high pressure and temperature <cit.>,
as well as many alloys, including those
of importance in phase change random access memory materials <cit.>.

These liquids play important roles in fuel cells, catalysis, and electrochemistry, where the dynamic
covalent bonds are expected to play an importance role in chemical reactivity <cit.>.

Many of these liquids are also found in planetary cores <cit.>
and understanding their structure and dynamics is of importance to planetary and geophysical sciences.




In all of these fluids, a complete understanding of their properties requires knowledge of the fundamental
interactions and timescales governing their chemical and physical transformations. 

The relative abundance of dynamic covalent bonds is expected to play a role in determining
the thermodynamic properties of the liquid state, as well as the kinetics of phase transformations.

For example, liquid silicon displays many of the hallmark anomalies found in water, because
the covalent bonds in liquid silicon lead to tetrahedral structures, analogous to hydrogen-bonding
in water.




Due to the importance of silicon to the semiconductor industry, and technology as a whole, l-Si
is well-characterized.

However, past research has mainly focused on the structure and thermodynamics of liquid silicon (l-Si).

In the solid state at ambient conditions, silicon is a covalently-bonded semiconductor in the diamond lattice,
and upon melting it undergoes a semiconductor-to-metal transition <cit.>.

Early simulations indeed predicted that l-Si is metallic, in agreement with experiments,
but they also uncovered a non-negligible
fraction of covalent bonds that persist in the liquid state <cit.>. 

The existence of these remnants of the solid phase were
later confirmed through a combination of computer simulations and Compton scattering experiments <cit.>.




The formation of covalent bonds in the disordered liquid state indicates the presence of a competition
between metallic and covalent interactions in silicon.

Indeed, a quantitative description of the solid phases of Si necessitates an accurate model
for the balance of metallic and covalent interatomic interactions <cit.>.

This competition between metallic and covalent interactions is also predicted to underlie
a metallic-to-semimetallic liquid-liquid phase transition in silicon <cit.>. 

In this case, a high density metallic l-Si that is dominated by metallic bonding can transition to a low density
semimetallic liquid, in which the interatomic interactions are predominantly covalent bonds, albeit transient ones
that readily break and reform in response to thermal fluctuations <cit.>. 





Despite these significant investigations into the structure and thermodynamics of silicon,
quantification of the lifetimes of the transient covalent bonds in l-Si is lacking.

To address this issue, we use ab initio molecular dynamics (AIMD) simulations to characterize covalent bonding kinetics
in l-Si, which serves as a prototypical liquid-state system with dynamic covalent bonds.

After discussing simulation details in the next section, we quantify the structure of
transient covalent bonds in l-Si and present a geometric covalent bond definition. 

We then use this definition to quantify the dynamics of covalent bonding in l-Si
and conclude with a discussion of future directions.













§ SIMULATION DETAILS



We simulated liquid Si at T=1800 K using the CP2K software package
following our previous work <cit.>.

The energies and forces in the MD simulations were evaluated using the QUICKSTEP module <cit.>,
which employs basis sets of Gaussian-type orbitals and plane waves for the electron density.

We used pseudopotentials, to represent the core electrons, and basis sets parameterized by
Godecker-Teter-Hutter (GTH):
GTH-PADE and GTH double-ζ, single polarization (GTH-DZVP), respectively <cit.>.

We explicitly treated the valence electrons of the 216 Si atoms using the strongly-constrained
and appropriately normed (SCAN) meta-generalized gradient approximation (meta-GGA)
density functional <cit.>, as implemented in  version 4.0.1 <cit.>,
with a plane wave cutoff of 650 Ry.

Initial configurations were taken from extensively equilibrated
simulations performed in earlier work <cit.>.

The systems were then further equilibrated at a constant temperature of T=1800 K,
maintained using the canonical velocity rescaling thermostat <cit.>.

Dynamic properties were computed from simulations in the microcanonical ensemble
using a timestep of 0.5 fs.









§ COVALENT BOND STRUCTURE IN LIQUID SILICON




We characterize the covalent bond structure in l-Si through the calculation
of maximally localized Wannier functions (MLWFs) and their centers (MLWFCs) <cit.>. 

MLWFs, in essence, can act as analogs of molecular orbitals for periodic systems, such as crystalline
and amorphous solids and liquids, and provide a useful, local picture of chemical bonding <cit.>.

We use the MLWFCs to represent the position of electrons, as done previously <cit.>,
in order to quantify electron-nuclei correlations
and develop a geometric criterion defining the existence of a covalent bond in an atomic configuration.

We note that although the MLWFs themselves are not unique, the MLWFCs are invariant
with respect to the choice of gauge within a lattice vector, which is a time-independent constant
for simulations in constant volume ensembles, like the microcanonical ensemble used here <cit.>. 

Thus, the MLWFCs can be used to define single bonds in each configuration.




We assume that covalent (single) bonds can be accurately defined by considering two- and three-body
correlations among Si nuclei and MLWFCs (W). 

The relevant two-body correlations are encoded in pair distribution functions,
g(r), involving Si nuclei and MLWFCs.

Three-body correlations are captured by the probability distributions P(θ),
where θ is the Si-Si-MLWFC angle formed by an Si atom, its nearest neighbor Si atom,
and a MLWFC.

Both sets of distribution functions are shown in Figure <ref>.





The pair distribution function for Si-Si and Si-MLWFC correlations, shown in Fig. <ref>a,
are consistent with the formation of covalent bonds between Si atoms.

The first peak in the Si-W g(r) is located halfway to the first Si-Si peak, and subsequent peaks
in the two correlation functions are out of phase.

Note that the small peak in the Si-W g(r) at distances less than roughly 1 Å is
consistent with the existence of non-bonded, lone pair electrons, visually depicted in Figures <ref>a-c.


	


Triplet correlations, as quantified by P(θ), show that a significant fraction of nearest-neighbor MLWFCs are consistent
with covalent bonding, evidenced by the peak in P(θ) as θ approaches zero (Fig. <ref>b), consistent
with linear Si-MLWFC-Si arrangements.

Large values of θ correspond to Si-Si-W triplets not involved in covalent bonds.

We additionally consider the possibility that lone pair-like MLWFCs may also appear at low values of θ.

Thus, we place the additional constraint on the Si-Si-W angle that the Si-W distance must be between 0.75 Åand 1.75 Å to avoid counting lone pair MLWFCs, resulting in the angle θ_t.

This more tightly-constrained angle removes contributions from lone pairs to the angular
distribution P(θ_t), Fig. <ref>b.

With this additional constraint, we suggest a reasonable geometric definition of a Si-Si covalent bond
corresponds to a Si-Si distance less than 3 Å and θ_t<30.

One such covalent bond is highlighted in Fig. <ref>c.









In the snapshot shown in Fig. <ref>c, as well as those in Figs. <ref>a-c,
one can observe a range of bonding and coordination environments. 

In particular, the number of MLWFCs associated with a single Si atom varies,
and, similarly, the number of Si atoms associated with a MLWFC also varies. 

Characterizing the fluctuations in Si and MLWFC coordination structures
quantifies the probability of forming lone pairs, covalent bond pairs, and metallic or diffuse pairs <cit.>. 

However, positional correlations alone do not suffice to characterize the nature of the MLWFs, as discussed
above in the context of angular correlations. 

Lone pairs and covalent bond MLWFs are high localized in space, such that their spreads are small,
such that bonded MLWFs are generally more localized than lone pair MLWFs <cit.>.

In contrast, the spreads of metallic or diffuse MLWFs are large, corresponding to delocalized pairs.

Therefore, we quantify the coordination structure of Si atoms and MLWFCs through distances,
as well as a range of MLWF spreads. 

We introduce a MLWF spread cutoff, σ_ W, such that MLWFs with spreads above this value are not included
in averages; σ_ W=∞ indicates that all MLWFCs are included in the calculations.





We start by examining the impact of σ_ W on the Si-MLWFC pair distribution function,
g_ Si-W(r), shown in Fig. <ref>a.

Introducing a finite σ_ W removes the contribution of diffuse MLWFs from the pair
distribution function.

A spread cutoff of 5 Å^2 reduces the peak at small distances and slightly increases
the first major peak, as well as the second peak. 

Further decreasing σ_ W to 3 Å^2 removes diffuse and nearly all lone pairs,
evidence by an absence of a peak for r<1 Å. 

The MLWFs remaining in the calculation are predominantly covalently bonded pairs, 
and the intensity of the first major peak in g_ Si-W(r), corresponding to bonded pairs,
increasingly significantly.

The second peak increases slightly as well. 

Therefore, σ_ W can be tuned to remove diffuse and lone pair MLWFCs,
and this tuning has the mainly impacts the structure of the first coordination shell, r<1.75 Å.




The probabilities of observing N_ W MLWFCs in the first coordination shell of an Si atom
for various σ_ W are shown in Fig. <ref>b.

When all MLWFCs are included, N_ W ranges from 1-7, and a maximum is observed at
N_ W=4, as expect from the sp^3 hybridization of the Si atoms.

As the cutoff is decreased, P(N_ W) shifts toward lower values, as delocalized MLWFs
are removed from consideration.

For σ_ W=5 Å^2, the maximum shifts to N_ W=3, with similar probabilities
at N_ W=2 and N_ W=4.

The appearance of finite probability for N_ W=0 is consistent with Si atoms that have all their
electrons in diffuse, metallic states.

Further reduction of σ_ W to 3 Å^2, which limits
the set of MLWFCs almost entirely to covalently bonded pairs, shifts the maximum of P(N_ W)
to zero, consistent with a small fraction (roughly 30 percent) of Si atoms involved in covalent bonds. 

In this limit, P(N_ W) essentially corresponds to probability of an Si atom
having N_ W covalent bonds.

This spans zero to four bonds per Si atom, and monotonically decreases with N_ W,
indicating that Si atoms fully coordinated by covalent bonds are less probable than partially
and non-bonded Si atoms.

This is consistent with the metallic nature of liquid Si.










We now turn our attention to the coordination structure of the MLWFCs, particularly
those that are in the first coordination of a Si atom.

To quantify this, we compute the joint probability distribution, P(N_ Si, N_ W),
corresponding to the probability that a central Si atom is coordinated by N_ W MLWFCs
and one of those coordinated MLWFCs is in turn coordinated by N_ Si Si atoms
(including the central Si atom). 

A schematic for N_ W=4 and N_ Si=2 for one MLWFC is shown in Fig. <ref>a. 

For ease of visualization, we focus on the conditional probability,

    P(N_ Si | N_ W)=P(N_ Si, N_ W)/P(N_ W),

physically corresponding to the probability that a MLWFC in the coordination shell of a central Si
atom is coordinated by N_ Si Si atoms, given that the central Si has N_ W MLWFCs
in its first coordination shell. 

These distributions are shown in Fig. <ref>b-c for varying σ_ W.





When all MLWFCs are included in P(N_ Si | N_ W), the conditional probability
has a maximum at N_ Si=2 for all values of N_ W>1. 

This maximum corresponds to a covalent bond between the central Si and
a neighboring Si atom.

Reducing N_ W from 6 to 2 primarily impacts the probably of observing singly coordinated
MLWFCs (to the central Si), which increases as N_ W is lowered.

For N_ W=1, P(N_ Si | N_ W)
differs significantly from the rest, with N_ Si=1 being most probable,
indicating that the MLWFC is most likely a diffuse or lone pair. 





Introducing a finite MLWF spread cutoff of σ_ W=5 Å^2
removes diffuse pairs and results in distributions that are similar for all N_ Si,
except N_ Si=1.

As N_ W is decreased, the probability of observing singly-coordinated, lone pair
MLWFCs increases.

Further reduction of σ_ W to 3 Å^2 results in sharp
P(N_ Si | N_ W) distributions that are nearly independent of N_ W.

These distributions span 1≤ N_ Si≤3, with a large maximum at N_ Si=2,
corresponding to covalent bond MLWFCs.

We note that there is a small probability for lone pairs when σ_ W=3 Å^2,
but lowering the spread cutoff further will increasingly remove these MLWFCs
and select only covalently bonded pairs <cit.>.





In summary, the local (bonding) coordination structure of Si atoms significantly fluctuates
and involves diffuse metallic, lone pair, and covalently-bonded states.

We have characterized these states, and the use of a MLWF spread cutoff, σ_ W,
can be used to systematically tune the involvement of these states when necessary.

Alternatively, Si-Si-W angular correlations, in addition to Si-Si and Si-W distances, used
in the above-described geometric definition of covalent bonds can be used to uncover
covalently-bonded MLWFCs in a similar manner, because linear Si-W-Si structures
are consistent with N_ Si=2 covalently-bonded pairs.









§ COVALENT BOND DYNAMICS IN LIQUID SILICON
	





Our simulations suggest that covalent bonds in l-Si rapidly break and reform on sub-picosecond timescales.

Figures <ref>a-c highlight one such covalent bond breakage and reformation event.

There, we show the time evolution of the MLWFCs (blue spheres)
for the atoms involved in the highlighted covalent bond exchange.

Initially (a), the central yellow Si atom is bonded to the left pink Si atom; a linear Si-MLWFC-Si
structure indicates a single covalent bond expected from the sp^3 hybridization of Si.

This bond breaks at a later time (b) due to thermal fluctuations, before the MLWFC of the central Si rotates
and forms a new bond with the rightmost red Si (c). 





	
		


Using the geometric definition of a covalent bond described in the previous section,
we are able to quantify the kinetics of Si-Si bond breakage,
in a manner analogous to conventional approaches to characterizing hydrogen bond dynamics in water <cit.>,
and, more recently, halogen bond dynamics in solid and liquid chlorine <cit.>.

To do so, we define an indicator function, h(t), which is equal to one when a covalent bond is present and zero otherwise.

We quantify covalent bond kinetics through the reactive flux correlation function <cit.>

    k(t)=-d C(t)/dt=-ḣ(0)1-h(t)/h,

where C(t) is the time correlation function (TCF) characterizing covalent bond lifetimes,

    C(t)=h(t) h(0)/h.








	


The reactive flux correlation function, k(t), is shown in Figure <ref>d.

At short times, vibrational and librational motion manifest the non-trivial, transient behavior
on scales less than roughly 0.3 ps. 

Beyond this time period, k(t) decays in a manner consistent with first-order kinetics,
k(t)∼τ^-1exp(t/τ), where τ≈ 0.6 ps is the covalent bond lifetime
estimated from fitting to the long-time behavior of k(t), shown as a dashed line
in Fig. <ref>d.





The covalent bond lifetime estimated above is similar to other significant timescales in l-Si. 

Orientational correlations, the velocity autocorrelation function,
and self-intermediate scattering functions all decay on timescales
similar to the covalent bond lifetime <cit.>. 

The agreement among various structural relaxation times and the covalent bond lifetime suggests that
covalent bond breakage is a limiting step for structural relaxation. 

This may be expected from the strength of covalent bonds, as well as their propensity to create
local tetrahedral order in the liquid state. 

 






§ CONCLUSIONS
	



In this work, we have presented an approach to estimate the lifetimes of transient covalent bonds
in condensed phase systems with an application to liquid silicon. 

The general strategy utilizes a maximally localized Wannier function (MLWF) approach to chemical bonding,
such that the existence of a covalent bond can be defined using two-body Si-Si correlations
and three-body correlations involving two Si nuclei and the center of a MLWF. 

With this geometric, ab initio definition of a covalent bond, we can estimate covalent bond lifetimes
in the liquid state from equilibrium simulations using the reactive flux formalism.

For liquid silicon at 1800 K and ambient pressure, we estimate a covalent bond lifetime of τ≈ 0.6 ps.

We note, however, that the generic concept of monitoring bond dynamics using a geometric, electronic structure-based
bond definition is not limited to equilibrium and could be used to monitor covalent bond dynamics in
melting or freezing processes at the focus of
laster melting experiments <cit.>, for example.




The covalent bond lifetime in liquid silicon may be measured using time-dependent scattering measurements,
such a time-dependent Compton scattering. 

Compton scattering has been used to shed light on the average bonding properties of metallic liquids <cit.>.

Extensions of this technique to the time domain are expected to uncover similar information about the average
dynamic properties of transient covalent bonds, like the lifetime at the focus of this work <cit.>.




Finally, we note that the covalent bond lifetime in l-Si is similar to bond lifetimes in other conventional
liquids with directional attractive interactions.

The hydrogen bond lifetime in liquid water is on the picosecond timescale <cit.>,
as is the halogen bond lifetime in liquid Cl_2 <cit.>.

Despite the vast differences in these liquids, their directional bonds all share a common thread:
the strength of the isolated bond, which is significantly weakened in the condensed phase, is on the order of 10
at the temperature of the respective liquids <cit.>.

Thus, thermal fluctuations in each of these different liquids are large enough
to cause the directional attractions to exist only fleetingly,
highlighting qualitative similarities of directional bond dynamics in the liquid state.



This work was supported as part of the
Center for Complex Materials from First Principles (CCM), an Energy Frontier
Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy
Sciences under Award #DE-SC0012575.
Computational resources were supported in part by the National Science Foundation
through major research instrumentation grant number 1625061
and by the US ARL under contract number W911NF-16-2-0189.






