Microscopic modeling of confined crystal growth and dissolution

We extend the (1+1)-dimensional fluid solid-on-solid (SOS) model to include a confining flat surface opposite to the SOS surface subject to a constant load. This load is balanced by a repulsive surface-surface interaction given by an ansatz which agrees with known analytical solutions in the limit of two separated flat surfaces. Mechanical equilibrium is imposed at all times by repositioning the confining surface. By the use of kinetic Monte Carlo (KMC) we calculate how the equilibrium concentration (deposition rate) depends on the applied load, and find it to reproduce analytical thermodynamics independent of the parameters of the interaction ansatz. We also study the dependency between the surface roughness and the saturation level as we vary the surface tension, and expand on previous analyses of the asymmetry between growth and dissolution by parametrizing the linear growth rate constant for growth and dissolution separately. We find the presence of a confining surface to affect the speed of growth and dissolution equally.

Figure 1

The confinement model illustrated. A periodic array of $L$ heights $h_i\left(t\right)$ makes up the crystal surface. All heights are relative to a fixed arbitrary point (here a zero point below the surface) along the vertical axis. A flat surface, subject to a constant applied force $F_0$, is located above the crystal surface at a height $h_l\left(t\right)$ at which the repulsive surface-surface force $F_{\psi }$ cancels out $F_0$ (mechanical equilibrium for all times $t$). We here skip some function arguments to increase the readability. The free energy $G_{\psi }$ depends on the surface configuration and $h_l$, and is generally unknown. In this work we describe this term by an ansatz (introduced in Sec. 2b) which depends only on the local separations $d_i$. The two open arrows denoted by 1 and 2 describe the two types of allowed transitions from a current surface configuration, that is, dissolution and deposition, respectively. Filled arrows represents forces.

Figure 2

Simulations done in equilibrium for a system with size $L=256$ for different values of $P_{\lambda }\propto F_0/L$, where $F_0/L$ is the applied pressure. For these simulations, we used ${\sigma }_0=1$ and ${\lambda }_D=2$. The solid lines represent the closed form solutions to the unconfined system. The statistical errors are increasing with $1/\alpha$, but overall small enough to not impact the trends shown here. The left axis label represents the calculated quantity, and for $P_{\lambda }\to 0$ (no confinement), the right axis label represents the respective thermodynamic quantity. From the top and middle panels we see that for low pressures, the average vertical crystal surface size $\langle s_{\uparrow \downarrow }\rangle$ and its fluctuation appear linear in the logarithmic scale of $1/\alpha$ with slopes close to 1, indicating that they are proportional to $1/\alpha$. This ultimately results in a constant heat capacity $C_V/k$ in the limit of no confinement, seen in the bottom panel. It is clear from these results that the effects of increasing the pressure $F_0/L$ is to flatten the surface and dampen the surface fluctuations.

Figure 3

Results of parametrization of the thermodynamic relation from Eqs. (27) and (28). In (a) and (b) we used ${\sigma }_0=1.0$ and ${\lambda }_D=1.61$. The dashed lines represent the result of the linear fitting of the data, whose slope are used to estimate the ratio between the quantities along the vertical and horizontal axes. In (c), the result of averaging the ratio of the axes in (b) for several ${\sigma }_0\in [0.1,5)$ is shown together with the analytical solution given by Eq. (28). We observe a good fit for ${\lambda }_D>0.5$, indicating that the model reproduce basic thermodynamics very well. A total of $\sim 3\times {10}^5$ equilibration simulations were done to generate the bottom plot.

Figure 4

Steady state convergence of the $\sigma \left(\mathbf{h}\right)$ from Eq. (19), the square root of the temporal surface length variance $\sigma \left(s\right)$, and the average distance between the crystal and the confining surface $\langle d_i\rangle$. The time unit of $1/\nu$ comes from the rate expression in Eq. (1). The smooth nature of the graphs are a result of combining 1000 simulations of the same system. For these simulations we used $L=256,\mathrm{\Omega }=-3/4,\alpha =E_b/kT=3.0,P_{\lambda }=0.01,{\lambda }_D=4.0$, and ${\sigma }_0=1.0$. We clearly see that the values converge to a plateau value, which is what we refer to as the steady state value. The steady state value of the calculated quantities are placed in the plots.

Figure 5

Steady state surface roughness $\sigma \left(\mathbf{h}\right)$ from Eq. (19) as a function of the saturation level $\mathrm{\Omega }$ for different values of $\alpha =E_b/kT$. The simulations were performed using $L=256,P_{\lambda }=0.01,{\lambda }_D=4.0$, and ${\sigma }_0=1.0$. We see that as we increase $\alpha$, the overall magnitude of the surface roughness is decreasing. For growth $\left(\mathrm{\Omega }>0\right)$, the roughness tends to increase with supersaturation. For dissolution $\left(\mathrm{\Omega }<0\right)$, we can either have a decrease or an increase in the roughness by increasing $\mathrm{\Omega }$, depending on the value of $\alpha$.

Figure 6

Results showing the linear growth rate $\dot{H}$ vs the saturation level $\mathrm{\Omega }=c/c_{\mathrm{eq}}-1$. Analyses of (a) indicate the existence of two separate linear relationships for growth $\left(\mathrm{\Omega }>0\right)$ and dissolution $\left(\mathrm{\Omega }<0\right)$. The slope represents the linear growth rate constant $k_g$. Without confinement $\left(P_{\lambda }=0\right)$, the trends presented here are identical to those found for $\alpha <0.75$ in Ref. [56]. In (b), the results of the dependency of $k_g$ on $\alpha$ and $P_{\lambda }$ are shown for $\mathrm{\Omega }>0$. The plot suggests the existence of a common point $P^{\prime }$ along the $P_{\lambda }$ axis where $ln{k}_g$ is independent of $\alpha$, and a linear dependency of $\alpha$ outside this. This point has been estimated to be 0.14 for growth and 0.53 for dissolution. For values of $\alpha <1$, the corresponding lines start to deviate from this common point in a nontrivial manner. The slopes of the lines from (b) vs $\alpha$ are shown in (c). This shows the slope to be very close to $\alpha$, indicating that $ln{k}_g\sim \alpha (P_{\lambda }-P^{\prime })$. This is general for all values of ${\sigma }_0$ and for ${\lambda }_D>2$. These simulations were performed using $L=256,{\sigma }_0=1.05$ and ${\lambda }_D=7/3$. In (a), we used $\alpha =4/3$.

Figure 7

The ten first iterations of ${\gamma }_k$ (squares) with iteration steps $k$ for a $L=256$ and $\alpha =1$ basic model without confinement. For this run we used ${10}^4$ cycles per iteration and $N=100$ total iterations. In practice, a good enough estimate is obtained after around ten iterations, however, in order to get a fair estimate of the spread, additional iterations are performed. The initial value is ${\gamma }_0=1$, and the exact solution is ${\gamma }_{\mathrm{eq}}=0$. The final result gave ${\gamma }_{\mathrm{eq}}=0.00\left(6\right)$.