Computational Study of a Multistep Height Model

Some time ago, two of the authors introduced a statistical mechanical model of a random surface embedded in a three dimensional space, suspended above a planar substrate with which it interacts Abraham and Newman (1988, 1991); Abraham et al. (1995). A detailed specification permitted the exact mapping of configurations of the three-dimensional surface onto those of the planar Ising model in its Peierls-contour representation. When treated by equilibrium statistical mechanics, this model displayed a phase transition; at low enough temperatures, typical configurations indicate that the surface is bound to the substrate. On raising the temperature sufficiently, the interface unbinds from the substrate and a layer of the bulk phase intercalates between the random surface and the substrate. Although it might be thought, as the authors originally did, that this model affords an example of wetting in 2-d Abraham (1986), it turns out that the exact critical indices obtained are quite unlike those normally associated with wetting. Rather, this model is much more accurately considered as an example of either Stransky-Krastanov (SK) Stranski and Krastanow (1939), or Volmer-Weber (VW) Volmer and Weber (1926) behavior.

We now go back sixty years to the seminal work of Burton, Cabrera and Frank Burton et al. (1951), who were the first to propose that a surface phase transition in a 3-d uniaxial system, say the spin-1/2 Ising model, should display singular behavior like the 2-d Ising model. They considered symmetry-breaking boundary conditions for which boundary spins are fixed to point up in the upper half space and down in the lower half space. They reasoned that below the 3-d critical temperature the upper half space will be in the up magnetized state and the lower half space in the down magnetized state. In the mean field approximation, the 2-d lattice of spins at the interface between the oppositely magnetized bulk phases is subjected to a mean field that cancels out. Hence this layer should be strictly 2-d in zero magnetic field and behave accordingly. This scenario is known to be incorrect; the BCF transition as originally conceived is actually a roughening transition Weeks et al. (1973) with quite different characteristics. The work in Abraham and Newman (1988); Abraham et al. (1995) affords examples in which the Burton, Cabrera and Frank type of transition is obtained by exact analysis, but in a different scenario and without making a mean field approximation.

Our first step will be to describe the model as originally formulated Abraham and Newman (1988, 1991); Abraham et al. (1995) and to review the results obtained for it. We will then indicate why we think the SK or VW scenarios are more appropriate than the wetting one. After that, we will point out a number of limitations that were necessary to obtain exact results and then show how Monte Carlo simulations can be used to get information when these limitations are relaxed and the exact result route is no longer available to us.

In the spirit of Kossel and Stransky (KS) Kossel (1927); Stranski (1928), configurations of molecules adsorbed on the substrate plane are constructed by regarding each molecule as a unit cube the lower side of which meshes with the underlying simple square lattice of the substrate, denoted $\Lambda\subset Z^{2}$. Molecular rafts are assembled as close packed arrays of the KS cubes. It is clear that, because it has no interior holes, a raft can be described just as well as a simple closed loop on $\Lambda$. Loops can meet at vertices of $\Lambda$ but not on edges, as this would imply over-counting. The upper faces of the KS cubes have height 1. To extend this model further, we permit placing rafts on top of rafts, without overhangs. In this way, we encounter loops within loops on $\Lambda$. As an additional restriction, the significance of which will soon become apparent, we do not allow any edge of a raft to lie directly above one of any lower raft, thus excluding multiple height jumps. This model thus manifests stepped towers erected on the plane; it was termed the “multi-ziggurat” (MZ) Abraham et al. (1995) model in recognition of its architectural precursor in ancient Sumeria. To complete the definition of the MZ model for equilibrium statistical mechanics, we must specify the configurational energy. The surface of a collection of molecular rafts is composed of two parts, one in contact with the substrate and the other in contact with the surroundings. Representing this collection as the equivalent loop form, labelled $\Gamma$, the energy $E(\Gamma)$ is given by

where $\tau$ is the surface tension, $\epsilon$ is the adhesion energy of a molecule to the substrate, and $L(\Gamma)$ is sum of the lengths $\ell(\gamma)$ of the simple closed loops of $\Gamma$,

The term $A_{1}$ is the area of contact of the plaquettes with the substrate, which is exactly the same as the “roof” area because of the ziggurat construction. Thus the second term in Eq. (1) is the energetic contribution of plaquettes in $\Gamma$ with normal perpendicular to the substrate plane. The remainder of the surface tension contribution is given by the first term. For stability against detachment, we evidently require that $\tau>-\epsilon$. We now briefly mention the results that can be obtained for this model. Notice that when $\tau=\epsilon$ we have recaptured precisely the planar Ising model, so the transition temperature and the free energy singularity are known. This is also the case if $\epsilon>\tau$; then the first layer is completely covered and subsequent rafts are laid on top of this layer as with $\tau=\epsilon$. In the remaining region of stability, characterised by $-\tau<\epsilon<\tau$, much less detailed information is available Abraham and Newman (1991).

If we are to regard the MZ model as a version of the VW and SK scenarios, then we should point out three shortcomings of the model as it stands. The first, which we will not discuss further here, is that we do not allow for the energy of mismatch, elastic in origin, between the first raft and the substrate. The second deficiency, the examination of which is the main subject of this paper, is the Ehrlich-Schwoebel (ES) Ehrlich and Hudda (1966); Schwoebel (1969) phenomenon. Simply stated, multiple height steps should be allowed, but they are energetically disfavored. In the next section, we will describe in some detail a model of Abraham and Newman (2009) that is simply related to the MZ model but permits multiple height steps and incorporates the ES idea.

The third problem is the formation of corrals. A corral or hole is a region of lesser height surrounded by a region of greater height. The model studied in this paper forbids corrals. Detecting a corral requires non-local information and it is unlikely that a physically reasonable equilibrium model would contain long range interactions that forbid corrals. On the other hand, dynamical considerations may suppress corrals. The formation of a corral from a plateau requires the removal of a molecule that is entirely surrounded by neighbors. Such a molecule will be tightly bound and its evaporation suppressed. A corral might also be formed by the sequential deposition of a wall that eventually encloses a region. However, surface diffusion would tend to favor more compact structures and suppress the formation of walls. Thus, although our model is an equilibrium statistical mechanical model, we believe that the no-corrals rule makes it an appropriate model for non-equilibrium surface growth processes.

From the point of view of equilibrium statistical mechanics, the above height models can be regarded as generalizations of the height representation of the 2-d Ising model, which corresponds to Eq. (1) with $\tau=\epsilon$. Other generalizations are possible. In this work we consider the generalization of Abraham and Newman (2009) that incorporates the ES idea by allowing height steps greater than one. Height steps greater than one are given an extra energy proportional to $\epsilon_{1}$. We study this multistep height model numerically and find that, along its critical curve, the critical exponents vary continuously with $\epsilon_{1}$. Though our model is motivated by surface physics, it appears to be closely related to another class of models in statistical physics, the O($n$) loop models. In the O($n$) loop models the statistical weight depends on the total loop length $L(\Gamma)$ and also contains a factor $n$ for each simple loop. Thus, the energy $E(\Gamma)$ becomes

where $C(\Gamma)$ is the number of simple loops in $\Gamma$ and $n$ is referred to as the loop fugacity. On lattices with vertices of degree 3, the loops are disjoint and on the honeycomb lattice one has the well-known O$(n)$ loop model introduced in Domany et al. (1981). A great deal is known about this O$(n)$ loop model when $n\leq 2$. For a given $n$, there exist three distinct phases: a disordered phase with small loops, a densely-packed phase, and a fully-packed phase. The densely-packed and fully-packed phases are both critical–i.e., the probability for two points to be on the same loop decays algebraically with distance. Between the disordered and the densely-packed phase is a critical curve as a function of $n$. The singular behavior of the O$(n)$ loop model along this critical curve can be described by a set of critical exponents that are functions of $n$. These exponents can be obtained from a mapping to Coulomb-gas theory Nienhuis (1984). Surprisingly, we find strong numerical evidence that the multistep height model studied here can also be described by the Coulomb gas theory and that there is a one-to-one mapping for universal quantities between $n$ and $\epsilon_{1}$. In the continuum scaling limit when the lattice spacing shrinks to zero, critical $O(n)$ loop models and hence presumably also critical multistep height models should be conformally invariant and described by Conformal Loop Ensembles — see, e.g., Schramm et al. (2009) — with a parameter value $\kappa$ mapped onto $n$ and $\epsilon_{1}$.

The plan for this paper is as follows. In Sec. II, we define the multistep height model. In Sec. III, we list quantities of interest for the model. In Sec. IV, we describe the numerical methods that are used in our simulations. In Sec. V, we present our results for the critical behavior of the multistep height model. The paper closes with a discussion in Sec. VI.

The height models proposed in Abraham and Newman (2009) and studied here generalize the height representation of the Ising model. Consider the two-dimensional Ising model on a square lattice in the spin representation with uniform plus-spin fixed boundary conditions. There is a one-to-one mapping from spins on the direct lattice to heights on the dual lattice. The height of any dual lattice site is the least number of Peierls contours that must be crossed on a path from the boundary. This definition leads to several important consequences. First, the magnetization of the Ising model is simply the number of even height sites minus the number of odd height sites. Second, there is a constraint that no holes or corrals are allowed in the height representation. That is, from every dual lattice site, there exists a path to the boundary that is non-increasing in height. Finally, the height representation of the Ising model will only have height steps of $+1$ or $-1$.

To generalize the model, we break the constraint of single height steps by allowing larger height steps at an extra energy cost, while still disallowing holes or corrals. The probability for a loop configuration is obtained from the energy, $E(\Gamma)$,

Here, $\Gamma$ is a collection of simple closed loops defined on the dual lattice. In the Ising model, these form the Peierls contours, while in the multistep model they can have a more complicated structure. In particular, this model allows Peierls contours to overlap, corresponding to larger height steps. $L(\Gamma)$, previously defined in Eq. (2), is the sum of the lengths of all of the loops (sum of all the edge weights), including now the lengths of all of the overlapping loops. In the loop representation $\tau$ is the energy associated with adding a single edge to $\Gamma$ while $N(1,1)$ is defined as $L(\Gamma)$ minus the length of all edges with weight one (step size one). $\epsilon_{1}$ is the energy cost associated with these larger height steps. $L(\Gamma)$ adds factors of $\tau$ linearly in the size of the height step, while $N(1,1)$ adds factors of $\epsilon_{1}$ linearly in the size of the step but only for height steps larger than one. Figure 1 shows a typical configuration of the multistep height model near the critical temperature for $\tau=2$, $\epsilon=0$, and system size $L=30$. Note that we have not included the term $(\tau-\epsilon)A_{1}(\Gamma)$ in Eq. (1) in the multistep height. It would be interesting to consider its effect in future studies.

Equation (4) completely specifies the loop representation of the model. The height representation is uniquely obtained from the loop representation. The height associated with a lattice site is obtained using the rule that the height change across a dual edge is equal to the weight of the bond. The direction of the height change is determined by the no-corrals rule. Figure 2 shows a possible edge configuration and the associated heights. In this configuration $L(\Gamma)=62$, $N(1,1)=6$ and $E(\Gamma)=62\tau+6\epsilon_{1}$.

In this paper we simulate multistep height models for several values of $\epsilon_{1}$ in the range $-1\leq\epsilon_{1}\leq\infty$ and $\tau=2$. Equation (4) with $\tau=2$ and $\epsilon_{1}=\infty$ corresponds to the height representation of the Ising model with interaction energy $J=1$.

Order parameters for height models de Mendonça (1999); Mazzeo et al. (1992); Alon et al. (1998) can be constructed in analogy to the magnetization for spin models. In the height representation of the Ising model, plus spins correspond to even heights while minus spins correspond to odd heights. One can define magnetization-like quantities, $M_{n}$ for positive integers $n$ as is done in Alon et al. (1998),

where $h(j)$ is the height at site $j$ and $N$ is the number of sites. This family of order parameters is motivated by the idea that in the rough (high temperature) phase, any height is nearly equally likely to occur so $M_{n}\rightarrow 0$. In the smooth (low temperature) phase, a single height will dominate and $M_{n}\rightarrow 1$ as $T\rightarrow 0$. The magnetization $m$ for the Ising model in the spin representation is equal to $M_{1}$. In the critical region of the multistep model, we find that only $M_{1}$ is useful for the small systems studied here. Heights significantly greater than two do not occur often, therefore $M_{n}$ will have strong finite-size corrections for $n>2$. Henceforth we use $m$ to represent either the magnetization for the spin representation of the Ising model or $M_{1}$ for height representations.

We measured the following quantities obtained from the height field $h(j)$:

Binder cumulant

The Binder cumulant $U$ is defined as

$$U=1-\frac{\langle m^{4}\rangle}{3\langle m^{2}\rangle^{2}}.$$

(6)

Crossings of the Binder cumulants as a function of temperature for various system sizes are used to locate the critical point.

Order parameter

From the finite-size scaling of the order parameter $\langle m\rangle$ at the critical temperature we obtain the exponent ratio $\beta/\nu$ from $\langle m\rangle\sim L^{-\beta/\nu}$ where $L$ is the system length.

Susceptibility

The susceptibility is defined as

$$\chi=\beta N\langle m^{2}-\langle m\rangle^{2}\rangle.$$

(7)

From the finite-size scaling of the susceptibility $\chi$ at the critical temperature we obtain the exponent ratio $\gamma/\nu$ from $\chi\sim L^{\gamma/\nu}$. At the critical temperature for fixed boundary conditions the finite-size scaling relation for $\chi$ is better behaved if the mean is set to zero, the infinite system critical value. Our finite-size scaling fits for $\gamma/\nu$ are made with the subtraction term $\langle m\rangle^{2}$ omitted.

Specific heat

The specific heat $c$ is defined as

$$c=\frac{\beta^{2}}{N}(\langle E^{2}\rangle-\langle E\rangle^{2}).$$

(8)

From the finite-size scaling of the specific heat $c$ at the critical temperature we obtain the exponent ratio $\alpha/\nu$ from $c\sim L^{\alpha/\nu}$.

Bare substrate areal fraction

The areal fraction of the bare substrate $\phi$ is defined as

$$\phi=\frac{1}{N}\sum_{j=1}^{N}\delta_{h(j),0}.$$

(9)

For $T\ll T_{c}$ there are almost no ad-atoms and $\phi\approx 1$ while for $T\gg T_{c}$, $\phi\approx 0$. The bare substrate is also known as the “gasket” in Schramm et al. (2009) where its (mean) fractal dimension at criticality is calculated rigorously for conformal loop ensembles and agrees with earlier non-rigorous results of Duplantier Duplantier (1990) for the area of connected regions in critical O$(n)$ loop models. Following the notation of Ref. Liu et al. (2011) the fractal dimension of these domains is $2-\beta^{\prime}/\nu$ so that the finite-size scaling of $\phi$ is given by $\phi\sim L^{-\beta^{\prime}/\nu}$.

Height

We measured the height at the origin $h(0)$ and the average height of the lattice $\bar{h}=(1/N)\sum_{i}h(i)$. These quantities are expected to diverge logarthmically in the system size.

We developed two Monte Carlo algorithms to simulate these height models. For the Ising case, we simulate the spin representation using the Wolff algorithm adapted to fixed boundary conditions. The spins are then mapped to heights in order to measure properties of the height model. For the multistep height model, there is no spin representation and we use a variant of the worm algorithm to sample the collection of closed loops $\Gamma$. The non-locality of the no-corrals rule creates significant computational difficulties in implementing the worm algorithm for the multistep height model.

We have analyzed a height model that generalizes the height representation of the Ising model by allowing height steps greater than unity with an energy cost parameterized by $\epsilon_{1}$. Our data supports the hypothesis that the critical exponents and prefactor $B$ depend continuously on $\epsilon_{1}$. Although we believe that the system has continuously varying critical exponents, it is conceivable that there are only one or two universality classes and that the $\epsilon_{1}$ dependent exponents reflect a finite-size crossover. The hypothesis of continuously varying exponents is strengthened by the relationship between the measured exponents that allows a mapping via the Coulomb gas parameter $g$ to the O$(n)$ loop model for $n$ in the range $1/2\leq n\leq 1$. We conjecture an exact relationship between $\epsilon_{1}$ and $g$ or $n$ that merits further investigation.

In the multistep height model the energy of a height step is a linear function of its magnitude. It would be interesting to study other energy functions. For example, the energy of a height step could be quadratic in its magnitude. We suspect that this would also yield models in the Coulomb gas universality classes.

A central feature of the model studied here is the no-corrals rule. It would be interesting to study a model with the same energetics but without forbidding corrals. Height models that allow corrals are generally referred to as “solid-on-solid” models. If the energy for height steps greater than one is quadratic and corrals are allowed the corresponding discrete Gaussian model is known to be in the $O(2)$ universality class Nienhuis (1987). The body-centered solid-on-solid (BCSOS) model restricts height steps to one and has an energy that is a sum of next-nearest neighbor height differences. The BCSOS model is exactly solvable and is also in the $O(2)$ universality class van Beijeren (1977); Baxter (2008). In the loop representation, allowing corrals corresponds to oriented loops with both orientations allowed. In the case of oriented loops on the honeycomb lattice with $\epsilon_{1}=\infty$, the loops are disjoint and the orientation degrees of freedom can be simply integrated out. This yields a statistical weight of $2$ for each loop, and hence maps to the $O(n)$ loop model with $n=2$. On this basis, it is tempting to speculate that for $0<\epsilon_{1}<\infty$, these oriented loop models will map to $O(n)$ loop models with loop fugacity in the range $1<n<2$.