# Temperature Dependence of Facet Ridges in Crystal Surfaces.

## I Introduction

Equilibrium crystal shapes typically consist of various flat facets connected by rounded surfaces. The way this structure changes with temperature has been studied for many years . Most aspects are well established. For example, rounded parts are associated with crystal orientations where the surface is rough, and flat facets can disappear from the equilibrium shape at the roughening temperature of those specific facets. Roughening transitions belong to the Kosterlitz-Thouless (KT) universality class . Many details follow from the scaling properties of the KT transition, e.g., that the facet diameter vanishes exponentially with reduced temperature, and that the curvature of the surface is universal just above the KT transition at that facet orientation.

Flat facets typically connect smoothly to rounded, rough surface areas. The surface orientation in the rounded part varies continuously and connects without a jump in angle onto the flat facet. Such phase boundaries are described by Pokrovsky-Talapov (PT) transitions . The angle difference vanishes with a $`\frac{3}{2}`$ power. The most salient feature of PT transitions is their anisotropic scaling. The rounded surface on the rough side of the PT transition can be interpreted as a stepped surface from the perspective of the flat facet. The characteristic lengths along, $`\xi _{}`$, and perpendicular to these steps, $`\xi _{}`$, scale as $`\xi _{}\xi _{}^z`$, with $`z=2`$.

Two dimensional equilibrium statistical mechanical systems are mathematically equivalent to one dimensional quantum systems at zero temperature. The transfer matrix formalism of the first is the path integral representation of the second. PT transitions are equivalent to metal-insulator transitions in one dimensional fermion systems. The steps represent the world lines of fermions, and the lattice direction parallel to the steps plays the role of time. In the rough phase (the metal phase) the fermions behave relativistically, i.e., time, $`\xi _{}`$, scales compared to space, $`\xi _{}`$, with dynamic exponent $`z=1`$, $`\xi _{}\xi _{}`$. At the PT metal-insulator transition the fermions behave non-relativistically, with dynamic exponent $`z=2`$ instead of $`z=1`$ .

The rough-to-facet interface can also be first-order, with a discontinuity in the surface orientation angle instead of the PT-type smooth connection. The critical endpoint between such first-order and PT-type segments provides an example of a quantum field theory with dynamic exponent $`z=3`$. To our knowledge an explicit example of such a PT critical endpoint (PTE) has not been reported before. We will do so in this paper.

Our current interest in these issues stems from the discovery that the scaling properties at facet ridge endpoints (FRE) are related to KPZ-type growth . FRE points are the endpoints of ridges where two facets connect directly. The ridge splits into two PT or two first-order transitions with a rounded, rough surface in between. The exactly soluble body-centered solid-on-solid (BCSOS) model on a square lattice (with only next nearest neighbour interactions) contains such a FRE point . Moreover, the transfer matrix of this model becomes stochastic at the FRE point and maps exactly onto the master equation time evolution operator for the so-called brick model in one dimension lower . The latter describes a growing 1D crystal. The 1D growing interface is rough, and its scaling properties belong to the KPZ universality class. The $``$-direction represents the spatial direction and the $``$-direction represents time, just as in mapping of the PT transtion to 1D fermions mentioned above. The dynamic exponent is equal to $`z=\frac{3}{2}`$. The characteristic time needed to reach the stationary growing state diverges as $`\xi _{}N^z`$, with $`N`$ the size of the 1D lattice.

This translates into non-trivial scaling properties of KPZ-type FRE points . When approached from the rounded phase in the direction parallel to the facet ridge, see Fig. 1, the surface orientation jumps. Moreover, at the rough side the characteristic length scale in the direction along the facet ridge, $`\xi _{}`$, scales as $`\xi _{}\xi _{}^z`$ with $`z=\frac{3}{2}`$.

This exact mathematical mapping between the scaling properties of equilibrium phase transitions to those of dynamic non-equilibrium transitions in one lower dimension is intriguing. But it leaves us with an important question. The above BCSOS model satisfies a special symmetry at its FRE point. The transfer matrix is “accidentally” stochastic. It has full rotational symmetry in the spin-$`\frac{1}{2}`$ operators (the step variables). This implies, e.g., a complete degeneracy of the free energy with respect to all surface orientations. Such a symmetry is too special. FRE points are a more general phenomenon, not limited to KPZ-type stochastic transfer matrices. We need to investigate whether the KPZ-type scaling exponent $`z=\frac{3}{2}`$ is generic. Maybe stochasticity is restored at large length scales, asymptotically close to the FRE point, with non-stochasticity an irrelevant operator. More likely, the FRE point changes its character. In this paper we present such an investigation, a numerical study of a model with a complex crystal shape.

We define this model in section 2. It originates from our earlier study of possible spontaneous low-angle faceting in FCC(111)-type facets . This model has the required property that its crystal shape is complex and that accidental symmetries are unlikely. Section 3 contains an overview of the phase diagram and crystal shapes. In section 4 we describe our numerical method.

The crystal structure is actually much richer than we anticipated. Entire facets vanish as functions of the coupling constants. FRE points become triple points. We also find a more exotic feature, like a first-order (FOR) line inside the rough phase emerging from the FRE point (section 5). At this line the surface orientation angle inside the rough phase jumps. In their paper about the exact solution of the FRE point in the square lattice BCSOS model, Bukman and Shore  speculate about the possible generalizations of Fig. 1. The FOR line structure corresponds to their so-called ridge scenario. However, they doubted that it would be realized. In our model, FOR lines are clearly present although they remain very short.

In a nearby section of the phase diagram, the FOR lines disappear. The rough-to-faceted ridge becomes sharp with a jump in the surface orientation angle. These first-order lines terminate further away from the FRE point in PTE points and the edge continues as a second order PT line. In section 6 we study the scaling properties at these PTE points. They display anisotropic scaling with exponent $`z=3`$. The conclusions are presented in section 7.

## II The Model

Our model emerges naturally at a coarse grained level in a study of spontaneous low angle faceting in FCC(111) facets. This will be described elsewhere . For clarity it is more appropriate to introduce the model from a different perspective. Consider the (001) facet of an HCP-type crystal, with $`ABABA`$-type stacking of triangular slices. The solid-on-solid (SOS) model description of such surfaces leads to a BCSOS model on a honeycomb lattice. The surface heights on the $`A`$-type triangular sublattice are even integers, $`h_A=0,\pm 2,\pm 4,\mathrm{}`$, and those on the $`B`$-type triangular sublattice are odd integers, $`h_B=\pm 1,\pm 3,\pm 5,\mathrm{}`$ The simplest Hamiltonian for such a surface is

$$H=\frac{1}{4}L\underset{(i,j)}{}(h_ih_j)^2$$
(1)   
with $`(i,j)`$ the summation over next nearest neighbour columns.

We break the symmetry between the $`A`$\- and $`B`$-type columns, such that two distinct (001) facets can be cut from this crystal, those with $`A`$’s or those with $`B`$’s on top. $`A`$ and $`B`$ could be different type of molecules, but that is not what we have in mind. The interactions between the next nearest neighbour columns, $`L_A`$ and $`L_B`$, would then be different. Instead, assume that $`A`$ and $`B`$ are identical types of atoms, but that for some reason the bulk structure is reconstructed such that the $`A`$’s move downward from the positions exactly in-between the $`B`$ layers. In that case the bonds between the $`A`$ atoms with the $`B`$’s in the layer beneath them are stronger than with the $`B`$’s in the layer above them. We can model this by adding a nearest neighbour coupling in the BCSOS model

$$H=K\underset{<i,j>}{}(h_ih_j)+\frac{1}{4}L\underset{(i,j)}{}(h_ih_j)^2$$
(2)   
where, in the sum over nearest neighbors, $`i`$ ($`j`$) always refers to $`A`$ ($`B`$) sublattice.

An alternative more generic scenario is the one realized in Naphthalene  (for a different crystal symmetry however). It involves only a surface reconstruction, no bulk reconstruction. Consider a crystal where the $`A`$’s and $`B^{}s`$ are identical molecules, but with an orientation degree of freedom pointing in two different directions. Inside the bulk these two orientations can be equivalent by (a slide-mirror-plane type) symmetry, but the presence of the surface typically breaks that equivalence. The surface molecules rotate due to surface relaxation. Typically they do so differently with the $`A`$’s or the $`B`$’s on top. This breaks the $`A`$ $`B`$ symmetry locally close to the surface and creates an energy difference, $`K0`$. The next nearest neighbour interactions remain virtually equal, $`L_AL_B=L`$, if the orientational aspects of the interactions decay rapidly with distance.

We introduce tilt energies $`E_i`$, $`i=1,2,3`$, to study the crystal structure at and nearby the central (001) facet orientation. The functional dependence of the free energy with respect to these $`E_i`$ yields the crystal shape by means of the Wulff construction .

$`=`$ $`{\displaystyle \frac{1}{4}}L{\displaystyle \underset{(i,j)}{}}(h_ih_j)^2+`$ (3)   
$`(K+E_1){\displaystyle \underset{<i,j>_1}{}}(h_ih_j)+`$ (4)   
$`(K+E_2){\displaystyle \underset{<i,j>_2}{}}(h_ih_j)+`$ (5)   
$`(K+E_3){\displaystyle \underset{<i,j>_3}{}}(h_ih_j).`$ (6)   
The summations run over nearest neighbour columns in the three different directions, see Fig. 2. The three tilt fields $`E_i`$ are related as

$`E_1`$ $`=`$ $`E_x`$ (7)   
$`E_2`$ $`=`$ $`{\displaystyle \frac{\sqrt{3}}{2}}E_y{\displaystyle \frac{1}{2}}E_x`$ (8)   
$`E_3`$ $`=`$ $`{\displaystyle \frac{\sqrt{3}}{2}}E_y{\displaystyle \frac{1}{2}}E_x`$ (9)   
$`E_x`$ and $`E_y`$ are conjugate to the tilt angles in the $`x`$\- and $`y`$-directions, i.e., parallel and perpendicular to one of the three main axes of the honeycomb lattice.

We simplify the model further by not allowing step excitations to touch each other. Fig. 3 shows all the remaining local configurations in this 10-vertex model. The arrows denote whether nearest neighbour $`AB`$ columns differ in height by $`+1`$ or $`1`$. The thin (fat) lines denote that the $`K`$ bond is in the low (high) energy state. The Boltzmann weight contributions are shown in Fig. 3 below each vertex configuration, with $`z_c=e^{2K2L}`$, $`z_s=e^{2K4L}`$, $`z_x=e^{E_x/4}`$, and $`z_y=e^{\sqrt{3}E_y/4}`$. The coupling constants are made dimensionless by absorbing the factor $`k_BT`$.

## III The Phase Diagram

The phase diagram of our model at zero tilt fields $`E_x=E_y=0`$, Eq.(2), is shown in Fig. 4. The (001) facet is stable at all temperatures inside the quadrant where $`K`$ and $`L`$ are both positive. The step excitations do not create sufficient entropy to roughen this surface even at $`T=\mathrm{}`$ (at point $`K=L=0`$). The roughening line lies inside the other three quadrants. We determine its location in the conventional manner for transfer matrix finite size scaling calculations . The step free energy $`\eta `$ in semi-infinite strips of width $`N`$ is equal to the difference in free energy for periodic and stepped boundary conditions in the finite direction. The KT roughening line is obtained by extrapolating the lines where $`N\eta =\frac{1}{4}\pi `$ for increasing strip widths.

At low temperatures for $`L<0`$ the surface spontaneously facets into three coexisting orientations. We refer to these as straight-step (SS) facets, shown in Fig. 5(a)to distinguish them from the three zig-zag (ZZ) structures, shown in Fig. 5(b), that appear elsewhere in the phase diagram at non-zero tilt fields $`E_i`$. The transition from the rough phase into this phase with three coexisting SS facets is first-order. The surface orientation angle jumps. The location of this transition line is determined by the methods described in section 4. Its location is not trivial, unlike the square lattice BCSOS model  where it is simply the meeting point of several $`E_i0`$ PT transition lines at $`E_i=0`$.

An additional ground state appears in the quadrant with $`K<0`$ and $`L>0`$, the triangles-state shown in Fig. 5(c). This state represents a highly degenerate surface reconstruction. Unfortunately this ground state is unstable against thermal fluctuations immediately at $`T>0`$ in this particular model. The surface reconstruction and associated preroughening phenomena associated with this type of ground state could have been quite interesting.

We determine the crystal shape at each value of $`K`$ and $`L`$ by introducing the tilt fields $`E_i`$, see Eq. (3). The technical details of these calculations are described in section 4. In this section we summarize the crystal shape development with temperature. We identified at least 5 distinct crystal structures, labeled I through V in Fig. 6, which shows at which values of $`K`$ and $`L`$ they are realized, as an overlay of the zero-tilt phase diagram Fig. 4.

In region I, located at $`L>0`$, the crystal structure consists of a non-tilted flat facet surrounded by a rough, rounded central region, and fully tilted facets of both types, SS and ZZ, at the edges (each repeated three times). This structure is shown in Fig. 7. All facets are separated from the rough region by PT lines. Their exact location can be calculated easily for the SS facets. Compare the free energy density of the faceted phase with that of a state with one line defect, e.g., as shown in Fig. 8. The PT transition occurs when the free energy of the faceted phase becomes equal to the free energy with the defect. The calculation is straightforward, but leads to a somewhat complex formula.

The PT lines in region I never touch or join. There are no FRE points between the SS and ZZ facets. The following analysis explains why. Suppose facet ridges between the ZZ and SS facets exist. Both faceted structures are frozen, without any fluctuations (in our model) and therefore the free energy in both is equal to the energy, $`E_{ZZ}=2K+2L\frac{1}{2}E_x`$ and $`E_{SS}=2K+4L\frac{1}{2}E_x\frac{1}{2}\sqrt{3}E_y`$. At a facet ridge these must be equal. This happens at $`E_y=\frac{4}{\sqrt{3}}L`$. However, the ZZ and SS facets roughen before this line is reached. Consider the zero temperature limit of the $`E_{ZZ}=E_{SS}`$ line. Line defects such as the one shown in Fig. 8 have the same energy as the two faceted configurations. The facets are unstable with respect to these defects, because at non-zero temperatures they gain entropy while the facets remain frozen. Two PT lines emerge from zero temperature and preempt the facet ridge.

As $`L`$ decreases, the ZZ facets shrink and disappear at $`L=0`$. At $`L0`$, the only remaining fully tilted facets are the three SS ones. Along $`L=0`$, they are still completely separated from each other, with rough, rounded orientations in between them, and PT transitions as borders. For $`L<0`$, facet ridges appear between the three SS facets. These facet ridges are capped by FRE points. The local structure at these FRE points changes with temperature in an exotic manner as shown in Fig. 9. The global crystal shape structure is less interesting. Initially it contains the central flat facet but this shrinks and vanishes across the roughening line.

Fig. 9(a) shows the local structure around the FRE point in region II. Two PT lines emerge from the FRE point and form the borders between the SS facets and the rough phase. Moreover, a FOR line emerges from the FRE point inside the rough phase. At this line the crystal orientation angle jumps. The skipped unstable surface orientations are associated with local ZZ facet type zig-zag step configurations. Those are unfavourable for $`L<0`$ and large $`K`$. The FOR line is definitely resolved in our numerical data, but remains extremely short (section 5). At its longest extension, it reaches only about two percent of the distance from the FRE point to the center. In some parts of region II it is so short that it can only barely be distinguished numerically, if at all.

The FOR lines disappear at the boundary between regions II and III. This boundary looks sharp in Fig. 6, but in reality it must be a smooth type of crossover. We discuss this in section 6.

In region III, see Fig. 6, the FOR line vanished. The border between the SS facets and the rough phase changed as well. Close to the FRE point the SS-rough boundaries are now first-order. They continue as PT lines beyond PTE-type critical endpoints. Fig. 9(b) shows this local structure.

The first-order segments become longer towards region IV. At the III-IV phase boundary, two PTE points from opposite FRE sides meet, and in region IV, the entire SS-rough phase boundary is first-order. In the mean time, the central rough region has been shrinking, and it disappears completely at the IV-V phase boundary.

## IV Numerical Methods

We determine the crystal shape by means of the transfer matrix method. This means that we calculate the exact free energies for semi-infinite lattices in terms of the largest and nearby eigenvalues of the transfer matrix. This is a standard method. Therefore we only need to comment here on the details concerning the surface tilt fields $`E_i`$, the amount of surface tilts $`Q_i`$ they create, and the Wulff construction (Legendre transformation) between these two types of quantities.

There are two natural choices for the directions in which the lattice is infinite. $`T`$ can be aligned with the $`y`$-axis or the $`x`$-axis, see Fig. 2. We refer to these as $`T_x`$ and $`T_y`$ type transfer matrices. In $`T_x`$, the infinite lattice direction coincides with the $`y`$-axis and the finite direction with the $`x`$-axis. In $`T_y`$ it is the other way around. Notice that for each choice there are three rotationally equivalent ones. For $`T_x`$, we are able to find the largest eigenvalues for $`4N18`$, where $`N`$ is the number of vertices in the time slice. For $`T_y`$, we can calculate the largest eigenvalues for $`2N10`$.

Consider the free energy as function of tilt angles $`f(Q_x,Q_y)`$. This is the free energy of a facet at a given orientation $`(Q_x,Q_y)`$. Some orientation ranges are thermodynamically unstable. These angles are not represented in the crystal shape. They are skipped, and result in sharp edges in the surface. To obtain the crystal shape one needs to minimize $`f(Q_x,Q_y)`$ under the constraint that the volume of the crystal (the amount of matter) is conserved. This leads to the famous Wulff construction  which is in essence a geometric construction for the Legendre transformation of $`f(Q_x,Q_y)`$. The tilt fields $`E_i`$ are the conjugate variables to the tilt angles $`Q_i`$. The crystal shape is a direct representation of the free energy function $`f(E_x,E_y)`$ at a certain temperature (at specific values of $`K`$ and $`L`$. The shape functions $`Q_x(E_x,E_y)`$ and $`Q_y(E_x,E_y)`$ tell us which surface orientations are represented in the crystal.

The transfer matrix provides us with a mixed version of this. It leads to $`f(Q_{},E_{})`$ and $`Q_{}(Q_{},E_{})`$. We control the tilt angle in the finite lattice direction, $`Q_{}`$, and the tilt field in the infinite direction, $`E_{}`$. $`Q_{}`$ sets itself. $`Q_{}`$ is set by the boundary condition in the finite lattice direction, $`h(x_{},x_{}+N)=h(x_{},x_{})+a`$. The total tilt of the surface can take the values $`a=0,\pm 2,\mathrm{}`$ up to $`N`$ when the transfer matrix is aligned with the $`y`$-axis. This means that an average slope $`Q_x=\sqrt{3}a/2N`$ is frozen into the surface in the $`x`$ direction. In the opposite set-up, with $`T_y`$, where the lattice is infinite in the $`x`$ direction, a tilt $`Q_y=a/2N`$ is frozen into the surface in the $`y`$-direction.

Consider the transfer matrix at zero tilt field $`E_{}=0`$ and a specific value of $`E_{}`$. We determine the largest eigenvalue, $`\lambda _0`$ of $`T`$ for each value of $`Q_{}`$. The free energy density is equal to $`f(E_{},0,Q_{})=\mathrm{log}(\lambda _0)/N`$. Fig. 10(a) shows a typical example, at $`K=0.24`$, $`L=0.29`$, and $`E_{}=0`$. Only the $`Q_{}>0`$ side is shown, because the curve is mirror symmetric. The discrete set of possible tilt angles $`Q_{}`$ increases with the strip width $`N`$, and we need to perform a finite size scaling (FSS) analysis to obtain the infinite-by-infinite lattice result. Fortunately these FSS corrections are reasonably small, and typically converge smoothly.

The tilt angle $`Q_{}`$ that is realized is the one that minimizes the free energy. This corresponds to the lowest point in Fig. 10(a). There is no need to repeat the calculation for other values $`E_{}0`$. They are related by $`f(E_{},E_{},Q_{})=f(E_{},0,Q_{})+E_{}Q_{}`$ because $`Q_{}`$ is conserved by the transfer matrix. In other words, the crystal shape $`f(E_{},E_{})`$ and the shape function $`Q_{}(E_{},E_{})`$ follow from the Legendre transform of $`f(E_{},Q_{})`$ with respect to $`Q_{}`$

It is easy to see from graphs like Fig. 10(a) how $`Q_{}`$ varies with $`E_{}`$. Simply rotate the page slowly and watch the minimum shift, and sometimes skip certain orientations. In Fig. 10(a) such skips do not take place. We follow the $`E_x`$-axis at $`E_{}=E_y=0`$ in Fig. 7 along which there are no orientational discontinuities. Notice that $`f(E_{}=0,Q_{})`$ develops a cusp at $`Q_x=0`$ with system size. This means that the central (001) facet is stable at small $`E_{}`$. The (001)facet boundary is a PT transition because the cusp is convex at $`Q_{}=0`$ and no orientations are being skipped for increasing $`E_{}`$ (rotate the page). The KT roughening transition takes place at values $`K`$ and $`L`$ where the cusp disappears. In the rough phase, $`f(E_{},Q_{})`$ has a quadratic minimum at $`Q_{}=0`$. Its curvature represents the roundness of the crystal shape at that orientation.

Next, consider the stability of the SS fully tilted facet. At maximum tilt angle the free energy $`f(E_{},Q_{})`$ has a definite slope. This means that the SS facet is stable for a specific range of $`E_{}=E_x`$. The absence of such a slope would imply that the SS facet is thermodynamically unstable, and would be absent in the crystal shape. The SS-rough phase boundary is a PT transition because $`f(E_{},Q_{})`$ is convex at $`Q_{}`$ near its maximum and no values of $`Q_{}`$ are being skipped by the Legendre transformation.

The above analysis is still incomplete because we did not check how $`Q_{}`$ is shifting simultaneously. The function $`Q_{}(E_{},Q_{})`$ is contained in the right and left eigenvector of the transfer matrix assocated with the largest eigenvalue. In the interpretation of $`T`$ as a quantum mechanical time evolution operator, $`Q_{}`$ is a velocity type off-diagonal expectation value.

Fig. 10(b) shows $`Q_{}(E_{},Q_{})=Q_y(E_y,Q_x)`$ at our example point. It shows that the orientation of the step excitations in the surface rotates as expected, from vertical just outside the $`Q_{}=Q_{}=0`$ flat facet to the correct slope just outside the SS fully tilted facet.

In the above example, we used the $`T_x`$ transfer matrix set-up to investigate the crystal shape in the horizontal direction. In that set-up $`Q_x`$ takes only discrete values, $`Q_x=\sqrt{3}a/2N`$. The opposite set-up with $`T_y`$ would have been better, since then $`Q_x`$ would have been continuous and $`Q_y`$ would have been rational. Fig. 11 illustrates this. We still use $`T_x`$ but now determine the crystal structure along the $`E_y`$ axis. The free energy curve lacks mirror symmetry. It shows the SS facet at the $`E_{}<0`$ side and the ZZ facet at the $`E_{}>0`$ side. It also shows the central facet around $`E_{}=0`$. The facet edges are smooth since they are PT-like.

From the above example it must be clear how we determined the development of the crystal shape with $`K`$ and $`L`$, described in section 3. It required numerous graphs of the type shown in Fig. 10 and Fig. 11. The results are quantitively represented in Fig. 6, Fig. 7, and Fig. 9. The latter two show the crystal shapes at actual values of $`K`$ and $`L`$ as specified in the figure captions. In the following two sections we limit the discussion to the most interesting features of the crystal shapes: the FOR lines, and the PTE points (see Fig.9).

## V region II: The FOR lines

In region II the crystal shape contains only SS-type, fully-tilted facets, with facet ridges in between them, capped by FRE points. The only structural difference with the KPZ-type FRE points in the square lattice BCSOS model is the presence of the FOR lines sticking into the rough phase, compare Fig. 9(a) with Fig. 1.

Such first-order ridges inside the rough phase have never been realized theoretically or experimentally before. Unfortunately, they remain very short in our model and are difficult to resolve numerically. Fig. 9(a) is drawn on scale for point $`K=1.43`$ and $`L=0.63`$ where the FOR line is at its longest, but remains very short compared to the distance of the FRE point to the center, $`\mathrm{\Delta }E/E0.02`$.

In the following we provide numerical evidence that FOR lines in our model are for real and not a fancy, transient, finite size scaling phenomenon. This is an important issue. Suppose the FOR lines vanish in the thermodynamic limit. Region II then represent a realization of the same type of FRE points as the ones in the square lattice BCSOS model, except that the transfer matrix is not stochastic anymore. We determined the anistropic scaling exponent $`z`$ from the leading, $`\lambda _0`$, and next leading, $`\lambda _1`$, eigenvalues of the transfer matrix in the sector with zero tilt $`Q_{}=0`$. The massgap $`m=\mathrm{log}(\lambda _0/\lambda _1)`$ represents the inverse of the time like correlation length, $`\xi _{}m^1N^z`$. We do this at our best estimates for the endpoint of the FOR line at several values of $`K`$ and $`L`$, and find $`z1.6\pm 0.1`$ (see Fig. 12). This is close to the KPZ exponent $`z=3/2`$. The assumption that the FOR lines vanish in the thermodynamic limit would confirm that KPZ scaling is valid beyond the stochastic subspace. However, in that interpretation, it remains unclear how far we are from the stochastic subspace. Crossover corrections to scaling might mask the “true” non-stochastic exponent $`z`$.

The existence of the FOR lines is much more likely, however. Fig. 13 shows the variation of the free energy $`f(Q_x)`$ at $`K=1.43`$ and $`L=0.63`$ for various $`E_y`$ close to the FRE point in Fig. 9(a). We need to use the $`T_x`$ set-up, the one with the infinite lattice direction aligned with the FOR line. The rational values of $`Q_y`$ in the $`T_y`$ set-up skip across the FOR line and make it invisible altogether (since it is so short). In the $`T_x`$ set-up the $`Q_x`$ are rational. In Fig. 13 we show only the data for $`N=18`$, otherwise the graph becomes too crowded.

Fig.14 shows the corresponding $`Q_y(Q_x)`$ behaviour at one value of $`E_y`$ (the curves for the other values of $`E_y`$ along the FOR line lie almost on top of this curve). Well inside the rough phase (below the FOR line in Fig. 9(a)), the $`f(Q_x)`$ curves are convex. The minimum is at $`Q_x=0`$, and it shifts smoothly with $`E_x`$ (rotate the page). Moreover, the curve remains convex all the way into the wings, implying PT-type boundaries with the SS facets.

In the opposite limit, well above the FRE point, the curves are concave. The two SS facets are coexisting at $`E_x=0`$, the entire $`f(Q_x)`$ curve except for its boundary values in the wings represent thermodynamically unstable orientations.

The intermediate behaviour of the $`f(Q_x)`$ curve determines the crystal shape close the FRE point. In the square lattice BCSOS model it flips at once from completely convex to completely concave (like changing the sign of $`a`$ in a polynomial like $`f(Q_x)=aQ_x^2`$). At that FRE point the transfer matrix is stochastic, implying that all surface orientations have the same free energy. In general, $`f(Q_x)`$ has more structure. There are many possible structures, but we observe in our model the simplest ones. In Fig. 13, the curvature of the central part of the $`f(Q_x)`$ curve changes before its SS facet wings come down. This creates the $`W`$-type shapes, and therefore gives rise to the FOR line. All of this happens in a very narrow $`E_y`$ interval. In Fig. 15 we show the location of $`E_y`$ where the curvature at the central part becomes zero, as function of $`N`$ for $`K=1.43`$ and $`L=0.63`$. To be more precise, we show the location of $`\lambda _0(Q_x)=\lambda _0(Q_x+\sqrt{3}/2N)`$ for $`Q_x=0`$, $`Q_x=\sqrt{3}/2N`$, and $`Q_x=2\sqrt{3}/2N`$. The curves pass the FRE point at small $`N`$, but then bend backward. Still, they do not converge to the FRE point (its location is known exactly and shown in the figure for reference). The curve for $`Q_x=0`$ is straightening out for system sizes $`12N18`$ and converges to a point below the FRE point at $`E_y=0.910\pm 0.005`$. We conclude from this that the FOR lines must be for real.

## VI Region III: Scaling at PTE points

In region III the $`f(Q_x)`$ curves near the FRE point are $`M`$ shaped (inverted $`W`$’s). This is shown in Fig.16 for $`K=0.37`$ and $`L=0.26`$. Fig.17 shows the corresponding $`Q_y(Q_x)`$ plots. The SS wings of the $`f(Q_x)`$ curve bend down and cross the central minimum before the curvature at $`Q_x=0`$ has a change to change sign. This behaviour implies that the $`f(Q_x)`$ curves become concave in the wings, and therefore that the PT transitions become first-order, see Fig. 9(b). Contrast this with the $`W`$ shapes in region II, Fig. 13, where the curvature of the central part of $`f(Q_x)`$ changes sign before its wings come down, and thus creates the FOR line.

The phase boundary between regions II and III, can not be sharp. It looks abrupt in Fig. 6, but that is an artifact of the shortness of the FOR lines. This changeover could take many forms. The FOR lines require the development of a concave segment in the center of $`f(Q_x,E_y)`$ while the PTE points require the development of concave segments in the wings. These two phenomena could be independent, and that would result in both a FOR line and first-order SS-rough segments, as shown in Fig.18(a). The FRE points become first-order as well. This means that subtle developments in the central part will be skipped and not be expressed in the crystal shape if they take place after the abrupt crossover from the central minimum to the wing minima. Our FOR lines are very short and move rapidly into the thermodynamically unstable region to vanish from the crystal shape. But this is only one of the possible scenarios.

We can not resolve accurately how the structure evolves from II to III. It happens too fast and the FOR lines are too short. We find some evidence that two concave bubbles move rapidly along the $`f(Q_x,E_y)`$ curves from the center to the wings in opposite directions. In that case the FOR line splits into two first-order lines, as shown in Fig. 18. Each line is still inside the rough phase, and each merges with the PT lines to form the first-order SS-rough segments.

The scaling behaviour at the PTE critical endpoints in region III deserves attention. The existence of such points has been anticipated  but not realized in solid-on-solid models before, and their scaling properties have not been confirmed numerically before.

It is advantageous to align the transfer matrix along the preferred direction of the steps. Therefore, we switch our attention to the SS-rough phase boundaries where, inside the SS facet, the steps run parallel to the $`x`$ axis, see Fig. 9(b). Inside this SS facet, the tilt vector is equal to $`\stackrel{}{Q}=(Q_x,Q_y)=(0,1)`$ We use the $`T_y`$ set-up, such that $`Q_x`$ is continuous.

Consider the crystal shape close to the PTE point. Define $`q_y=Q_y^{}Q_y0`$ and $`q_x=Q_xQ_x^{}=Q_x`$ as the deviation of the tilt with respect to the SS facet ($`Q_y^{}=1`$ and $`Q_x^{}=0`$). Let $`ϵ_y=E_yE_y^{}`$ and $`ϵ_x=E_x^{}E_x`$ be the distances from the PTE point The crystal is rough (faceted) when $`ϵ_y>0`$ ($`ϵ_y<0`$). The rough-to-faceted edge is PT-type (first-order) when $`ϵ_x>0`$ ($`ϵ_y<0`$). The free energy $`f(q_x,q_y)`$ as function of the tilt angles is an analytic function near PT transitions. This is a direct result of the dilute (free Fermion) nature of the “gas” of defects . Therefore we can expand it in terms of a polynomial in $`q_x`$ and $`q_y`$. Our presumption is that this remains true also at the PTE point.

$`f(q_x,q_y)=`$ $`a_1|q_x|+b_1ϵ_y|q_y|+`$ (10)   
$`{\displaystyle \frac{1}{3}}b_3ϵ_x|q_y|^3+{\displaystyle \frac{1}{4}}b_4|q_y|^4+\mathrm{}`$ (11)   
Notice the absence of the quadratic term. This is a well known property of the entropic hardcore repulsion of meandering defect lines (the missing steps in the SS structure) at low densities . Eq. 10 reproduces the usual PT scaling behaviour. The crystal shape follows from minimizing $`f(q_x,q_y)`$ with respect to $`q_x`$ and $`q_y`$ at constant $`ϵ_x>0`$ and $`ϵ_y>0`$. This gives $`q_x=0`$ and $`b_1ϵ_y+b_3ϵ_xq_y^2+b_4q_y^3=0`$. Along a path of constant $`ϵ_x>0`$ the tilt angle vanishes with a square root power, $`q_y|ϵ_y|^{\frac{1}{2}}`$ and the free energy has the famous PT power law singularity, $`f|ϵ_y|^{\frac{3}{2}}`$. This is the power with which the surface orientation connects to the SS facet. The anisotropic scaling exponent $`z=2`$, $`l_yl_x^z`$, follows from the $`fq_y^3`$ power with which the free energy scales at the PT points. This implies that the free energy of one single defect scales per unit length as $`ml_x[f(2\pi /l_x)f(0)]l_x^2`$ with $`l_x`$ a finite lattice length in the direction perpendicular to the defect. $`m`$ is inversely proportional to the correlation length in the $`y`$-direction, $`ml_y^1`$. Our numerical results are not shown here, since they are in complete agreement with this.

At the PTE point, $`ϵ_x=0`$, the same quantities scale according to Eq.(10) as $`q_y|ϵ_y|^{\frac{1}{3}}`$, $`fq_y^4`$, and $`z=3`$. Along the first-order line, $`ϵ_x<0`$, the jump in the orientation angle vanishes linearly, $`\mathrm{\Delta }q_y|ϵ_x|`$. $`q_x`$ develops a jump as well, but with a weaker power, because the rotation of the surface orientation inside the rough phase is an higher order effect.

It is not certain that $`q_y^4`$ is the lowest surviving power in the free energy at the PTE point. It could be the $`|q_y|^5`$ term instead. For example, in the conventional free fermion analysis for PT transitions the free energy is an odd function in $`|q_y|`$. In that case the critical exponents change to $`q_y|ϵ_y|^{\frac{1}{4}}`$, $`fq_y^5`$, and $`z=4`$.

In Fig. 19 we show the free energy as a function of $`q_y^4`$ at the PTE point for $`K=0.37`$, $`L=0.26`$ in region III for system sizes $`4N10`$. The curves are straight lines. This confirms that the free energy is going to zero as $`fq_y^4`$ and therefore that $`z=3`$.

The scaling behaviour $`q_y|ϵ_y|^{\frac{1}{3}}`$ across the PTE point and the power $`\mathrm{\Delta }q_y|ϵ_x|`$ with which the jump in $`q_y`$ along the first-order boundary vanishes are difficult to resolve numerically because of the discrete nature of $`q_y`$. We find evidence that the jump scales linearly, but the results are somewhat ambiguous.

At the III to IV phase boundary the PT segment along the rough-SS facet boundary shrinks to zero and the two PTE points merge. This can be described by a polynomial of the form

$`f(q_x,q_y)=`$ $`a_1|q_x|+b_1ϵ_y|q_y|+`$ (12)   
$`{\displaystyle \frac{1}{3}}b_3(cϵ_x^2)|q_y|^3+{\displaystyle \frac{1}{4}}b_4|q_y|^4+\mathrm{}`$ (13)   
$`ϵ_x`$ is now the distance along the rough-SS ridge from the point in the middle between the two PTE points, and $`c>0`$ ($`c<0`$) in region III (IV ). This polynomial, if correct, would imply that the critical exponents at the merging point, $`c=0`$, are the same as at the PTE points, $`z=3`$ and $`\beta =\frac{1}{3}`$, except that the jump in angle vanishes now as, as $`q_y|ϵ_x|^2`$. We find numerically that the free energy scales as $`fq_y^4`$ and, therefore, that indeed $`z=3`$.

## VII Conclusions

In this paper we obtained the crystal shape for a BCSOS model on a honeycomb lattice. This shape is quite complex. It contains two types of fully tilted facets, the SS and ZZ facets. The latter disappear for negative nearest neighbour interactions. In region III in that part of the phase diagram, we find that the rough-to-SS facet boundary becomes partially first-order with sharp edges where the surface orientation jumps. These first-order segments connect at PTE critical endpoints to conventional PT-type segments where the rough phase connects to the facet without a jump in orientation. The scaling properties at these PTE points follow from a simple higher-order polynomial free energy generalization of PT transitions, with anisotropic scaling with exponent $`z=3`$ instead of $`z=2`$.

In region II of the phase diagram the rough-SS phase boundary is PT-type everywhere, but a first-order line sticks out of the FRE point into the rough phase. These FOR lines remain very small in this model, but we resolved them numerically sufficiently to be sure they exist.

The main object of this study is to establish how general and universal the KPZ scaling properties of the FRE point in the square lattice BCSOS model are. Our conclusion is that the stochastic nature of their transfer matrix makes KPZ-type FRE points special and unstable. In general the local crystal shape around FRE points is more complex than simply one facet ridge splitting into two PT lines. Instead there are FOR lines sticking into the rough phase and/or the rough-facet phase boundaries become first-order. Even more complex structures (not observed in our model) are possible as well.

We are unable to resolve the scaling properties of the the critical point at the end of the FOR line (opposite to the FRE point). The FOR lines remain too short. It might still be KPZ-type scaling with $`z=\frac{3}{2}`$. In the exactly soluble square lattice BCSOS model, the KPZ nature of the scaling is closely linked to the degeneracy of the free energy with respect to all surface orientations (this is implied by the stochasticity of the transfer matrix). Our model lacks this special symmetry, but, at these critical points at the end of the FOR line, the free energies of small angle surface orientations become degenerate. Maybe this local version of that symmetry is sufficient for $`z=\frac{3}{2}`$ anisotropic scaling, or it could lead to a new, yet unknown value of $`z`$. This aspect needs further study.

This work is supported by NSF grant DMR-9700430.