Thermodynamic melting point

In order to study the phenomenon of premelting, we have to find the thermodynamic melting point $T_m$ of our model. At least three different methods for estimation of the melting point have been given in the literature. The most elegant and thermodynamically valid is based on the calculation of the Gibbs free energies of the bulk solid and liquid phases as a function of temperature at given pressure. The melting point is then found, by definition, by making the Gibbs free energies of solid and liquid equal [19]. Hoverer, this method is computationally quite elaborate and contains a fairly large source of inaccuracy, due to the fact that the Gibbs free energies usually depend quite weakly on temperature, and $T_m$ is then an intersection point of two flat curves. Therefore, it is quite difficult to find the position of the intersection point precisely and a small error in the obtained position of the point leads to large error in estimation of $T_m$.

In the second method, proposed by J.F. Lutchko et al. [77] one determines the temperature dependent velocity of the moving solid-to-liquid interface $v_{sl}(T)$ at $T>T_m$ and thus, by extrapolating this velocity to the zero $v(T_m)=0$ one can obtain the solid-liquid coexistence point, which is interpreted as $T_m$. This method has been shown to give a melting point comparable to what is obtained from the free energy method.

The third method, is computationally the most straightforward one. In this method, an excess amount of kinetic energy is supplied to the atoms of the surface region of the sample. This results in a temperature gradient and in a nucleation of liquid phase at the surface. Heating is then stopped and the sample is allowed to evolve back to equilibrium. The temperature at which a solid-liquid coexistence takes place is interpreted as the thermodynamic melting point $T_m$ [89]. All these methods can estimate $T_m$ with comparable accuracy.

We chose the second method and determined the velocity of the solid-melt interface at different temperatures (in the range of 2300 - 2500 K) in the samples with various low-index free surfaces. The value of $T_m$ varies with the geometry of the simulation cell, and especially with the number of the surface layers. To test such effects we performed simulations for the samples with the different number of layers, $N_l$, using 64 atoms in a layer (Va(111), Va(001)) and 72 atoms (Va(011)).

A temperature above $T_m$ (the upper limit of $T_m$ was found in the bulk simulations) was reached using the following method. We started the simulation at a low temperature ( usually at $T = 1800~K$) far below $T_m$, and increased the temperature at a high rate (by $\Delta T=50~ K$) after each 1000 steps. In this way the system could be superheated up to $T=2500~K$. A very simple method of the velocity rescaling in each MD step is used, e.g. the instantaneous kinetic energy is calculated as $E_{kin}=\sum_i mv^2_i/2$ and velocity of each atom is rescaled by the same factor to satisfy $E_{kin}=3/2k_BT$. This simple method of the velocity rescaling was chosen instead of the Nose-Hoover mechanism to significantly reduce the fluctuations of temperature. When the desired temperature is reached, the velocity of the solid-melt interface is determined by monitoring the decrease of the structure order parameters, $\eta$, as a function of time (See Fig. 5.8).

Figure: Order parameters of the superheated $Va(011)$ system vs. time (MD steps) at temperature $2420K$. The sample contains $2520$ atoms.

The order parameter was calculated only over the top 10 surface layers to avoid the ordering effects introduced by the static substrate. The temperature dependence of the velocity $v_{sl}$ is governed by the thermodynamic free-energy driving force, given by the free energy difference, $\Delta G=G_s-G_{liq}$, between the solid $G_s$ and the liquid $G_{liq}$ phases [90,91].

$$v_{sl}\sim \frac{D_{liq}d_L}{\Lambda^2}\left( 1-exp\left(\frac{\Delta G}{RT} \right)\right)$$ (6.1)

where $\Lambda$ is the mean free path in the liquid phase, $d_L$ is the thickness of the liquid film, $D_{liq}$ is the self-diffusion coefficient of atoms in the liquid phase. One can write with a good approximation [90] $\Delta G \approx (T-T_m)\Delta S$, where $\Delta S$ is entropy difference between the solid and liquid phases, and obtain at temperatures close to the melting point

$$v_{sl}\sim \frac{D_{liq}d_L}{\Lambda^2}\left( 1-exp\left(\fr... ...T_m)\Delta S }{RT} \right)\right) \approx const\times (T-T_m)$$ (6.2)

Using this expression for $v_{sl}(T)$ as a function of temperature, $T$, we calculate the thermodynamical melting point by a linear fit (See Fig 5.9). The temperature at which the propagation velocity vanishes is considered to be the thermodynamical melting point $T_m$ [89].

Figure: Velocity of propagation of the solid-to-liquid interface, $v_{sl}$, as a function of temperature and extrapolation to zero velocity. The $Va(111)$ sample, $2240$ atoms, $35$ layers.

The melting temperatures obtained for the samples with different number of surface layers and the various types of low-index faces are summarized in Table 5.2. The smaller the number of the layers, the larger the apparent $T_m$. This is due to the influence of the underlying static substrate, which stabilizes and orders the samples. Therefore, to obtain more accurate and reliable estimation for the melting point $T_m$ we use the sample with the largest number of free surface layers.

Table: The melting temperatures for the various surfaces of vanadium, calculated for the samples with different number of the surface layers, $N_l$.

$N_l$	$T_m~~Va(111)$	$T_m~~Va(001)$	$T_m~~Va(011)$
$23$	$2241 \pm 10 K$	$2230 \pm 10 K$	$2260 \pm 10 K$
$26$	$2234 \pm 10 K$	$2200 \pm 12 K$	$2252 \pm 10 K$
$29$	$2220 \pm 8~ K$	$2235 \pm 11 K$	$2260 \pm 12 K$
$32$	$2224 \pm 10 K$	$2237 \pm 8~ K$	$2270 \pm 11 K$
$35$	$2222 \pm 9~ K$	$2228 \pm 10 K$	$2240 \pm 6~ K$

The calculated melting temperature of the sample with the least packed face Va(111) is $T_m=2222 \pm 10 K$ is the close one to the experimental value $T_m=2183 K$. For other surfaces $T_m$ is insignificantly larger due to superheating effects on the solid-liquid interface just above the thermodynamical melting point $T_m$. This effects are very significant for the close packed Va(011) surface. This dependence of the thermodynamic melting point on the geometry of the crystal face is also observed in MD simulations when temperature is increased gradually up to the melting point. It was found the Va(111) and Va(001) surfaces melt approximately at the same temperature $T_m \approx 2230 \pm 30 ~K$, while the sample Va(111), with the close packed free surface, melts at more elevated temperature about $T_m \approx 2245 \pm 10 ~K$.

We want to point out again that in the bulk simulations, where we employed periodic boundary conditions in all directions, melting is observed to happen at a temperature $T_b=2500 \pm 5K$. This melting temperature is not the true thermodynamic melting point $T_m$, but a temperature at which mechanical instability arises, leading ultimately to a mechanical collapse of the crystal lattice. The temperature at which mechanical instability occurs can be determined by using the Born criterion at zero-external stress. This shows that the superheating region, $T_b-T_m$, in MD simulations of metals may be large enough (250 K in case of vanadium) and great care should be exercised while determining $T_m$ from the simulations results of bulk and surface melting.