Folded network and structural transition in molten tin

The fundamental relationships between the structure and properties of liquids are far from being well understood. For instance, the structural origins of many liquid anomalies still remain unclear, but liquid-liquid transitions (LLT) are believed to hold a key. However, experimental demonstrations of LLTs have been rather challenging. Here, we report experimental and theoretical evidence of a second-order-like LLT in molten tin, one which favors a percolating covalent bond network at high temperatures. The observed structural transition originates from the fluctuating metallic/covalent behavior of atomic bonding, and consequently a new paradigm of liquid structure emerges. The liquid structure, described in the form of a folded network, bridges two well-established structural models for disordered systems, i.e., the random packing of hard-spheres and a continuous random network, offering a large structural midground for liquids and glasses. Our findings provide an unparalleled physical picture of the atomic arrangement for a plethora of liquids, shedding light on the thermodynamic and dynamic anomalies of liquids but also entailing far-reaching implications for studying liquid polyamorphism and dynamical transitions in liquids.


Supplementary Note 1. Extraction of high-q structure factors of liquid Sn
Identifying subtle changes in the 3D liquid structures by means of 1D radial distribution function g(r), or its reciprocal-space counterpart, structure factor S(q), requires knowledge of g(r) or S(q) to a very high accuracy. In this work, the structure factors of l-Sn were derived from high-energy XRD and analyzed with an in-house computational code 1 . In what follows, we introduce the method employed in this work to extract the structure factors.
The diffraction intensity contains scattering signals from different sources, including background, scattering from the sample, and sample absorption effect: where a is the normalization factor to put the signal into atomic units, and I incoh is the incoherent (Compton) scattering from the atoms in the sample which can be computed using the analytical formula provided by Balyuzi 2 . For the empirical correction (which accounts for other spurious effects such as multiple scattering, fluorescence radiation, temperature diffuse scattering, etc.), we adopted a modified function due to Thijsse 3 : where a, b, c are parameters to be optimized through the "correctness" of the structure factor S(q), which is related to samp ( ) I q by: ( ) ( ) ( ) 2 coh I q S q f q = (4) where f is the atomic scattering factor for Sn, obtained from tabulated values.
The following three criteria were used to gauge the accuracy of S(q) (1) The high-q rule: At high-q values (e.g., q > 12 Å -1 ), S(q) should oscillate around 1. Therefore, we set 2 2 1 (2) The sum rule: According to Krogh-Moe 4 and Norman 5 , the structure factor and the number density ρ 0 have the following relation: A second cost function can thus be formulated to serve as an additional constraint: Here, if the number density ρ 0 is not known a priori, it can be optimized together with other parameters. In the present study, it was stipulated from our ab initio MD simulation.
(3) The Kaplow rule Kaplow et al. 6 observed that for small r, the reduced radial distribution function G(r) has Therefore, a third cost function was constructed to refine the optimization: where σ = 2 Å and Lumping the three criteria together, we defined a total cost function: (10) where w i are fitting weights. The parameters a, b, c, s, and ρ 0 were optimized by minimizing χ 2 using the simulated annealing method 7 in this work.
Based on this approach, we were able to obtain accurate high-resolution pair correlation functions of l-Sn, enabling us to probe into subtle structural changes, e.g., the peak-splitting of g ( where the atomic form factor ( ) f q was obtained from tabulated values, and is a normalization factor to bring the S(0) value at 530 K to match the experimental isothermal compressibility data available in the literature 8 .
From thermodynamic fluctuation theory 9,10 , the isothermal compressibility as a function of temperature is linked to the structure factor S(0) (i.e., density fluctuation at the long-wavelength limit) by the relation: where k B is the Boltzmann constant, T is the absolute temperature, and n is the atomic number density 11 . S(0) was obtained by extrapolating the structure factor in the range of 0.18 < q < 0.4 Å -1 to q = 0 using a quartic polynomial, as shown in Supplementary Fig. 8b.

Supplementary Note 3. Isothermal compressibility of liquid Sn derived from ab initio MD
Based on thermodynamic fluctuation theory 9,12 , the mean-square fluctuation in number density ρ ) at a constant temperature T is given by: where κ T is isothermal compressibility, k B is Boltzmann constant, T is temperature, and V is volume. The above equation can be expressed in terms of volume fluctuation: To compute of l-Sn, our goal was to obtain the mean-square fluctuation of volume, ( ) To this end, a large ensemble consisting of 288 atoms (e.g., V = 31.4 Å 3 at 600 K) was subjected to NPT ab initio MD simulations at desired temperatures. At each temperature, the volume was controlled by the method of Parrinello and Rahman 13,14 combined with a Langevin thermostat as implemented in VASP. It was found that the mean-square volume fluctuation follows a normal distribution and the result converges within approximately 200,000 timesteps (the data are not shown), from which the isothermal compressibility was computed and plotted in Fig. 2d. The error bars in Fig. 2d correspond to the standard deviations of volume fluctuations. 6 To cross-check the accuracy of the isothermal compressibility computed with the above method, we further computed κ T from the following relation: In doing so, we established how the ensemble pressure P changes with uniform dilation of l-Sn. At each temperature, by varying the volume from -5% to 5%, a series of AIMD simulations were carried out in NVT ensembles to obtain the statistical mean of the pressure at each volume. The isothermal compressibility was thus calculated by fitting the pressure and volume relationship employing a second-order Brich-Murnaghan equation-of-state 15 . The as-obtained isothermal compressibility data at selected temperatures were juxtaposed with the data obtained from the fluctuation theory in Fig. 2d, showing satisfactory agreement between the two methods.
One should note that the isothermal compressibility values obtained with AIMD are consistently larger than the experimentally reported ones. This discrepancy concerns primarily with the accuracy of the DFT treatment in dealing with physical properties in general. For instance, with the current DFT treatment (e.g., projector augmented-wave treatment and generalized-gradient approximation) 13,14 , the theoretically calculated bulk modulus of a-Sn at 300 K (~33 GPa, this work) is roughly 62 % of the experimental value (53 GPa, ref. 16 ), which justifies the differences observed between the experimental and theoretical isothermal compressibility.

Supplementary Note 4. Liquid excess entropy obtained with a multiple-histogram method (MH)
The multiple histogram method (MH) is a technique originally proposed by Ferrenberg and Swendsen [17][18][19] to obtain free energies and entropies from a discrete sampling of energy.
In a canonical ensemble, the distribution of potential energy E at a given temperature is given by This equation can also be rewritten as: is the configurational part of the Helmholtz free energy ( ) A T , which has the following expression: where ( ) Λ T is the de Broglie thermal wavelength 20 .
We list two equations based on which our histogram-based code was developed to derive the configurational entropy: The Ehrenfest definition of second-order phase transition requires that the second derivative of the free energy be discontinuous at the transition temperature. Therefore, in our analysis, we sought to derive a twice-differentiable free energy curve with respect to energy or temperature (strictly speaking, in the vicinity of the transition region). To achieve this, we expressed the configurational entropy with cubic spline interpolation 7 : where c k are spline parameters to be determined and E k are spline knots evenly distributed between two predefined energy bounds. In our spline representation of ( ) c S E , 15 knots were used (i.e., n = 16). Given the analytical form of the configurational entropy, the integral form of the free energy can be numerically assessed according to Supplementary Eq. 24. In this approach, the entropy can only be optimized within an additive constant (a shift of S c will be counterbalanced by the change of the integrated free-energy term following Supplementary Eq. 23), therefore, c 0 is arbitrarily fixed a priori.
To optimize the spline parameters c k , we defined a cost function by comparing the theoretical In practice, in order to improve the sampling statistics (the tails of the histograms are poorly sampled), we adopted a histogram reweighting technique, and the final cost function is redefined as: (27) where M refers to the number of independent temperatures at which NVT simulations were conducted, and N is the number of histogram bins. Supplementary Fig. 9 illustrates the histograms sampled at selected temperatures T j . We adopted the simulated annealing technique 7 to optimize the cost function.
To achieve thermodynamically convergent results, 18 temperatures were chosen to ensure adequate overlaps in the potential energy distributions. Satisfactory statistics on energy distribution were achieved by long-time ab initio MD simulation, ca. 300,000 time-steps at each temperature. MD simulations were performed in canonical ensembles (288 atoms) using the Nose-Hoover thermostat for temperature control 21 . The total number of state points for the thermodynamics analysis was tallied to 5×10 6 .
Having obtained the configurational entropy ( ) c S E , its temperature dependence could be derived Supplementary Fig. 10). Note that the canonical configurational entropy

Supplementary Note 5. Atomic relaxation dynamics in liquid Sn
The self-intermediate scattering function F s (q,t) characterizes the relaxation time of l-Sn. The F s (q,t) is defined as the Fourier transform of the van Hove function 22 . Instead of Fourier transformation, these functions can also be directly computed from the atomic trajectories 22 : where q is the wave vector, and N is the number of particles. The relaxation time τ a is defined as the time it takes for s 1 ( , ) F q t e = (see Fig. 3b in main text).

Supplementary Note 6. Chemical bonding analysis based on electron localization function
The interaction between the center atom and a coordination atom is considered a bond. The primary signature of a covalent bond is the sharing and localization of valence electrons between the bonding atoms. The degree of localization and the relative strength of atomic bonds can be assessed by the electron localization function 23,24 , where D σ is a measure of Pauli repulsion and scales with the probability density of finding another same-spin electron near the reference electron. 0 D σ is the D σ value in a homogeneous electron gas having the same local spin density. ELF = 0.5 represents the same level of Pauli repulsion as in the homogeneous electron gas, and a higher ELF value indicates that the electrons are more localized (e.g., ELF = 1.0 can be interpreted as perfect localization). We used VASP to quantify the ELF values 25,26 on a 100 × 100 × 100 grid in the cubic box. To characterize the bonding between any two atoms, we define the spatially averaged ELF value within an enclosed sphere Ω in the middle of two atoms: In this work, we set a radius of 0.35 Å for the enclosed sphere to obtain χ. Note that bond character value χ universally characterizes the bonding nature of the atoms regardless of the absolute electron charge density of the system in question. A similar method was previously adopted to define the chemical bonding in a chalcogenide glass 27 .

Supplementary Note 7. CRN structure for LDA-Sn versus FN structure for l-Sn
Group IVA elements Si and Ge were demonstrated to have tetrahedrally coordinated low-density amorphous (LDA) structures 28,29 , which can be described by the CRN model 30 . Knowledge of the structure of amorphous Sn, however, was not established. Since Sn belongs to the same group with s 2 p 2 valance electrons, it is thus intriguing to inquire whether the CRN type low-density Sn exists.
Currently, no experimental information is available on the structure of amorphous Sn, which can be synthesized in the laboratory, e.g., by vapor deposition onto very cold substrates 31,32 .
Based on first-principles modeling, we are able to demonstrate the existence of metastable LDA-Sn at low temperatures. We first constructed an ideal tetrahedrally-coordinated CRN structure of Sn employing the WWW bond-switching algorithm 30 . The as-prepared amorphous Sn was fully relaxed with ab initio molecular mechanics at a constant pressure of zero. The relaxed LDA-Sn was found to be mechanically stable and to maintain a nearly perfect tetrahedral coordination; that is, the fully expended random tetrahedral network of Sn. To examine the thermal stability of the LDA-Sn, the as-relaxed CRN structure was heat treated by means of constant-pressure (at zero pressure) ab initio MD at a heating rate of 5×10 12 K/s. Upon increasing temperature, the structure reaches its stability limit at 450 K, followed by a transition to high-density liquid (HDL) Sn accompanied by a volume collapse of ~20%. We now have two disordered structures of Sn, the LDA-Sn and the metallic HDL-Sn. Supplementary Fig. 19 displays the computed structure factor of LDA-Sn at 300K. Also compared are the tetrahedral order parameters for the CRN LDA-Sn and the folded-network structure of l-Sn at 450K (Supplementary Fig. 20).

Supplementary Note 8. Relaxation dynamics and lifetime of atomic bonds in liquid Sn
Atomic packing in l-Sn departs from the hard-sphere model and hence the Voronoi method 33 is no longer suitable for evaluating the coordination number (CN). Alternatively, the CN was directly obtained by counting the number of atoms within a cutoff distance. At the melting temperature of 580 K, Sn atoms are found to accommodate 8~9 atoms in the first atomic shell.
The lifetime of the atomic bonds was analyzed through the partial bond autocorrelation function, defined as: (31) where N(i,t 0 ) is the atom number of the i th nearest neighbor at time t 0 ; Θ(N,t 0 ) indicates whether atom N is a coordination atom at time t 0 . If yes, The angular bracket refers to the ensemble average. It can be seen that, for an arbitrary atom, if the i th nearest neighbor remains a coordination atom after time t, we say that the i th atomic bond persists within time t, that is, Supplementary Fig. 15a, the bond relaxation dynamics of l-Sn can be briefly summarized as follows: (1) The bond relaxation time, τ b , defined as ( ) e , is found to vary with bond lengths.
The covalent bonds exhibit a slightly longer relaxation time than the metallic bonds. (2) During short times t < τ b (t < 1.0 ps), the covalent bonds relax noticeably slower than the metallic bonds. (3) For extended times t > τ b , the relaxation dynamics of the long-lived metallic and covalent bonds converge, owing to the fluctuating bonding nature of the atomic bonds discussed below.
The lifetime distribution of the atomic bonds can be derived from which is plotted in Supplementary Fig. 15b as a probability profile. It can be seen that the distribution profiles for the covalent and metallic bonds show slight discrepancies in the short-time window but exhibit almost identical behavior with prolonged relaxation time. This is consistent with the observation that the chemical bonding has a bifurcated nature and fluctuates between covalent and metallic bonding, as confirmed by the cross bond-bond correlation function, defined as: (32) where superscript c denotes covalent bonding and superscript m denotes metallic bonding. Our analysis of ( ) Ψ c B t confirms the mutual transition between covalent and metallic bonding of l-Sn.
Interestingly, the bond relaxation dynamics of l-Sn resembles the bifurcated hydrogen bonding observed in water. A schematic diagram is drawn to show the transition from covalent to metallic bonding in l-Sn as in Supplementary Fig. 16. The heat capacity C P as a function of temperature from experiment and AIMD. The AIMD C P of l-Sn was derived from the slope of the enthalpy-temperature curves. A structural transition of l-Sn is inferred from the specific heat changes, which agree favorably with the experimental data. The density changes of the same simulations are given in Supplementary Fig. 5. (red solid line). In calculating the local bond-angle distribution (3) g ( ) local θ , we only consider the bonds that belong to the same tessellated tetrahedron in the first coordination shell (Note that the coordination polyhedron can be spatially tessellated into tetrahedra using the Voronoi method). As shown for the case of l-Cu where the atoms are closely packed, the g (3) (θ) displays a peak around 60°, typical of a HS-like packing (regular tetrahedron). On the other hand, for the covalent-bonded CRN structure, it shows a peak around the tetrahedral bond angle of 109.5°. By contrast, the g (3) (θ) of l-Sn shows a bimodal distribution at bond angles lower than 109.5° (a clear peak centered around 60° and a broad peak centered around 100°).   Fig. 18). For VASP simulations, we used the PW exchange-correlation functionals together with the PAW treatment. Valance electron configurations used in the AIMD are also provided.