# Atomic Layering at the Liquid Silicon Surface: a First-Principles Simulation

## Abstract

We simulate the liquid silicon surface with first-principles molecular dynamics in a slab geometry. We find that the atom-density profile presents a pronounced layering, similar to those observed in low-temperature liquid metals like Ga and Hg. The depth-dependent pair correlation function shows that the effect originates from directional bonding of Si atoms at the surface, and propagates into the bulk. The layering has no major effects in the electronic and dynamical properties of the system, that are very similar to those of bulk liquid Si. To our knowledge, this is the first study of a liquid surface by first-principles molecular dynamics.

Liquid metal surfaces have attracted much attention during the last years . Their particular properties, very different from those of non-metals, have been investigated by Ångström-resolution experiments, and simulated by different approaches. One of their most interesting features is the atomic layering, a density oscillation that originates at the sharp liquid-vapor interface, and extends several atomic diameters into the bulk. Although some experiments supported an increased surface density at $`l`$-Hg , and some kind of atomic layering was predicted theoretically , its existence has been demonstrated unambiguously only recently by X-Ray reflectivity in Hg , Ga  and Ga-In alloys . Most of the experiments have been done close to the low melting temperatures of these metals, but Regan et al  studied the Ga surface up to 170C. They found that capilarity waves strongly decrease the reflectivity peak heights, but not the peak widths, suggesting that the decay length of local layering is temperature-independent, and that surface layering is still present quite above the melting point.

Different explanations have been proposed to account for this effect . Rice et al  have argued that the abrupt decay of the delocalized electron density forms a flat potential barrier against which the ions lay orderly, like hard spheres against a hard wall. Tosatti et al  have used the glue model of metallic cohesion to argue that surface atoms, trying to effectively recover their optimal coordination, alternatively increase and decrease their density. Surface layering effects, like surface-enhanced smectic ordering, have also been observed in liquid crystals . In this case, its origin is the tendency of the highly nonspherical molecules to present a particular orientation towards the surface.

Liquid silicon ($`l`$-Si) is a rather peculiar system. Silicon transforms, at 1684 K, from a covalent semiconductor solid, with diamond structure and coordination 4, to a liquid metal. Experiments  and MD simulations  show that its coordination ($`67`$) is lower than that of typical liquids ($`12`$), due to the persistence of directional bonding in the liquid phase. In spite of the enormous literature on its solid surfaces, very little is known on the structure of the liquid silicon surface. Measurements are very difficult because of its high reactivity and melting temperature. Model calculations, and computer simulations with semiempirical potentials, are also difficult because of the mentioned coexistence of covalent and metallic bonding, and its unknown interplay at the surface.

In this letter we present a study of the $`l`$-Si surface by first-principles molecular dynamics (MD) simulation . This approach deals equally well with covalent and metallic bonding, and it is therefore very well suited for this problem. Electrons are treated by solving the Kohn-Sham  equations selfconsistently for each ionic configuration, using the local density approximation for exchange and correlation. The quantum mechanically obtained forces are then used to generate the classical trajectories of the ion cores. The calculations were performed with the SIESTA program  using a linear combination of numerical atomic orbitals as the basis set, and norm-conserving pseudopotentials . A uniform mesh with a planewave cutoff of 40 Ry is used to represent the electron density, the local part of the pseudopotential, and the Hartree and exchange-correlation potentials. Only the $`\mathrm{\Gamma }`$ $`k`$-point was used in the simulations, since previous work  found cell-size effects to be small. For the present calculation we used a minimal basis set of four orbitals (1 $`s`$ and 3 $`p`$) for each Si atom, with a cutoff radius of 2.65 Å. We have extensively checked the basis with static calculations of different crystalline Si phases and solid surfaces, and MD simulations of the bulk liquid . The energy differences between solid phases are described within 0.1 eV of other ab-initio calculations. The diamond structure has the lowest energy, with a lattice parameter of 5.46 Å (0.5 $`\%`$ larger than the experimental value). Adatom- and dimer-based (111) and (100) surface reconstructions found in other ab-initio calculations  are well reproduced, with geometries and relative energies changing less than $`0.1`$ Å and $``$0.15 eV when moving from a $`\mathrm{\Gamma }`$-point calculation with a minimal basis set, to a converged $`k`$-sampling with double-$`\zeta `$ and polarization orbitals. The structural, electronic and dynamical properties of $`l`$-Si are in good agreement with other ab-initio calculations at the same density and temperature . The calculated diffusion constant ($`1.5\times 10^4`$ cm<sup>2</sup>/s) is somewhat smaller than that obtained with a double-$`\zeta `$ or polarized basis ($`1.72.0\times 10^4`$ cm<sup>2</sup>/s) which is in agreement with other ab-initio simulations. We interpret that the minimal basis overestimates the energies of the saddle point configurations occurring during diffusion, but we consider that this is not critical for the present application. Also, we leave for a future work the inclusion of spin fluctuations, which affect significantly the diffusion constant but not the structural properties of the liquid .

We first perform a long simulation of bulk $`l`$-Si at $`T`$=1800 K , using a cubic 64-atom cell with periodic boundary conditions. The fixed cell size (10.58 Å) was adjusted to obtain zero mean pressure, and corresponds to a density 3$`\%`$ smaller than the experimental density near the melting point. We then construct our initial 96-atom slab by repeating one bulk unit cell in the $`x`$ and $`y`$ directions, and one and a half cells in the $`z`$ direction, plus 10 Å of vacuum. No particles leave the slab during the 30 ps simulation. After a relaxation of 10 ps, the system reaches equilibrium and the averaged magnitudes are essentially the same for the next or the last 10 ps, and for both sides of the slab. These long relaxation and observation times are required because the calculated density autocorrelation time at the surface ($`1`$ ps) is considerably longer than the typical bulk-liquid correlation times ($``$ 0.1 ps) . The average surface energy (836 $`\pm `$ 40 dyn/cm) is in good agreement with the experimental surface tension (850 dyn/cm at 1800 K) , suggesting a small entropic contribution.

In Fig. 1 (solid line) we present the ionic density profile $`\rho (z)`$. It shows a pronounced atomic layering, with similar features as those reported for the Ga surface . Like in that case, $`\rho (z)`$ can be fitted accurately by a sharp error function at the surface, and a sinusoidal wave with an exponential decay towards the bulk: by superimposing two such functions (not shown), for both surfaces, we obtain similar values of the parameters and an oscillation period of 2.5 Å. To check that the observed layering is not a sign of incipient crystallization, we have computed several magnitudes in the central region of the slab ($`|z|<3.7`$ Å). The radial and angular distribution functions, electronic and vibrational densities of states, and the diffusion constant, are all very similar to those of bulk $`l`$-Si, and have no resemblance of those in the solid phases. As an example, we compare in Fig. 2 the bond-angle distribution function for the bulk and slab simulations, using a bond cutoff distance of $`r_m`$=3.10 Å. For a better understanding of the origin of the layering, we compute the normalized density-density correlation function :

$$c_\rho (z_0,z)=\frac{\delta \rho (z_0,t)\delta \rho (z,t)}{\delta \rho ^2(z_0,t)^{\frac{1}{2}}\delta \rho ^2(z,t)^{\frac{1}{2}}},$$
(1)   
where $``$ denotes time average and $`\delta \rho `$ is the difference between the instantaneous density at time $`t`$, $`\rho (𝐫,t)=_{i=1}^N\delta (𝐫𝐫_i(t))`$, and the average density $`\rho (𝐫)=\rho (𝐫,t)`$, where $`\delta (𝐫)`$ is Dirac’s function. Fig. 1 also shows $`c_\rho (z_0,z)`$ for $`z_0`$ at the positions of the outermost peaks of each side. Its decaying oscillation is clearer than that of the density profile, all whose relevant features match very well with the superposition of the two $`c_\rho `$’s. The apparent lack of decay of $`\rho (z)`$ towards the interior is peculiar to our particular slab thickness, because the superposition is positive at the center of the slab, and negative at the two surfaces. Most important is, however, that the two surface-induced oscillations are clearly independent of each other ($`c_\rho `$’s out of phase), and incommensurate to the slab thickness. The density layering is thus an intrinsic surface feature and not a result of finite size effects.

In order to obtain information about the bond orientations at the surface we calculate the two-particle density:

$$\rho _2(𝐫_0;𝐫)=\frac{N/(N1)}{\rho (𝐫_0)}\underset{i=1}{\overset{N}{}}\underset{ji}{}\delta (𝐫_0𝐫_i)\delta (𝐫𝐫_j).$$
(2)   
To represent $`\rho _2(𝐫_0;𝐫)`$, we first average over the directions parallel to the surface:

$`\rho _2(z_0;z,x)`$ $`={\displaystyle \frac{1}{2\pi xA}}{\displaystyle d^3𝐫_0^{}d^3𝐫^{}\rho _2(𝐫_0^{};𝐫^{})\delta (z_0z_0^{})}`$ (4)   
$`\delta (zz^{})\delta (x\sqrt{(x^{}x_0^{})^2+(y^{}y_0^{})^2}),`$   
where $`A`$ is the area of the simulation cell. In Fig. 3 we show $`\rho _2(z_0;z,x)`$ for $`z_0`$ located at the three peaks of $`\rho (z)`$. Fig. 3(a) shows a clear tendency of surface atoms to form bonds parallel and normal to the surface. The height of the correlation peaks goes well beyond those of the density, which can be seen also in the figure (notice that $`\rho _2(𝐫_0;𝐫)\rho (𝐫)`$ for $`|𝐫𝐫_0|\mathrm{}`$). This shows that the bond-induced correlations are responsible for the layering of the density, and not the other way around. Fig. 3(b) shows a similar, but attenuated tendency, that disappears in the third layer (Fig. 3(c)), which already has a very symmetric, bulk-like pair correlation function.

Further insight can be obtained from the $`z`$-dependent coordination $`n(z)`$, defined as the average number of neighbors within a distance $`r_m`$. At the bulk, we obtain $`n`$=6.4, in agreement with the experimental value . The distribution of local coordinations (DLC) is also very close to those of other ab-initio calculations , showing a maximum at coordination 6. We can also use the bulk simulation to construct an ideally terminated surface, cutting abruptly the system at, say, $`z`$=0. We then find $`n(z)=4.3`$ at $`z=1.0`$ Å, which is the distance between the outermost peak and the inflection point in the slab density profile. At the outermost peaks of the actual slab, we obtain $`n(z)`$=5.3, and a DLC peaked at 5. These values show that surface structural rearrangements increase the coordination of the ideally terminated surface, reaching a value of only one neighbor less than in the bulk. If we associate coordination 6, in the bulk liquid, with an octahedral arrangement, a simple picture can be drawn, in which the surface atoms try to preserve their bulk environment while minimizing the number of broken bonds. As a consequence, the octahedra get oriented in the surface so that only one broken ‘bond’ points towards vacuum, with another bond towards the interior and four bonds laying on the surface plane. This picture is consistent with figures 3(a) and 3(d). In the latter, we have restricted the sum over $`i`$ in eq. (2) to particles having coordination 5, what results in even more pronounced peaks in the $`x`$ and $`z`$ directions. Also, we note that the maximum of $`\rho _2(z_0;z,0)`$ occurs at $`|zz_0|=2.5`$ Å, what explains the same period observed in the density profile. The same distance is found for the in-plane surface bonds (i.e. for $`\rho _2(z_0;z_0,x)`$), and for the bulk bonds. Thus, contrary to other metals, we do not find a shortening of the surface bonds, and the silicon surface layering seems to be related only to the bond orientations. However, it must be emphasized that the bond angle distribution (Fig. 2) is very wide, indicating a large variety of fluctuating atomic environments , so that our ‘oriented octahedra’ should be considered only as a very rough and qualitative picture.

An interesting question is whether the surface structural rearrangements produce a noticeable signature in the electronic structure. In Fig. 4, we compare the local density of states (LDOS) at the outermost peaks of $`\rho (z)`$ and at the center of the slab. Although $`k`$-sampling is important for a converged LDOS , we use here only the $`\mathrm{\Gamma }`$-point eigenvalues, to facilitate the comparison with previous work  and because we focus on its spatial variation. It can be seen that, apart from a slight narrowing, due to the diminished surface coordination, there are no major differences, suggesting that the surface and bulk atomic environments are rather similar, and again pointing towards bond orientations as responsible for surface layering.

In conclusion, we have performed the first study of a liquid surface by first-principles MD simulation. In spite of the high melting temperature of Si, we find a marked layering of the density near the surface, similar to those observed in other metals, like Ga and Hg, with low melting temperatures. However, the surface layering of Si seems to have an origin at least partially different from that in other metals, with remanent directional covalent bonding playing an essential role. In spite of the rather slow decay of the layering towards the bulk, the average structural, dynamical, and electronic properties converge very rapidly to their bulk liquid values. Although more converged simulations would be highly desirable in the future, we consider that this work provides a new qualitative understanding of the complex structure of liquid surfaces.

We acknowledge useful discussions with E. Chacón and M. Weissmann. This work was supported by Argentina’s CONICET and by Spain’s DGES grant PB-0202.