Defining shapes of two-dimensional crystals with undefinable edge energies

The equilibrium shape of crystals is a fundamental property of both aesthetic appeal and practical importance: the shape and its facets control the catalytic, light-emitting, sensing, magnetic and plasmonic behaviors. It is also a visible macro-manifestation of the underlying atomic-scale forces and chemical makeup, most conspicuous in two-dimensional (2D) materials of keen current interest. If the crystal surface/edge energy is known for different directions, its shape can be obtained by the geometric Wulff construction, a tenet of crystal physics; however, if symmetry is lacking, the crystal edge energy cannot be defined or calculated and thus its shape becomes elusive, presenting an insurmountable problem for theory. Here we show how one can proceed with auxiliary edge energies towards a constructive prediction, through well-planned computations, of a unique crystal shape. We demonstrate it for challenging materials such as SnSe, which is of C2v symmetry, and even AgNO2 of C1, which has no symmetry at all.

is fixed by well defined (and computable) physical quantity E11`, but absolute positions of the striplongitudinal and transverse (depending on value α) shifts, are both arbitrary.So, indeterminacy is 2 and the auxiliary parameter only translates the Wulff shape, without changing its form (trivial in this case).
Next, two cuts in different directions, is a more realistic case, since this allows one to cut out a finite material piece.Then 4 unknown edge-energies ε1, ε1` and ε2, ε2` satisfy 2 equations only ε1 + ε1` = E11`, and ε2 + ε2` = E22`, leaving again edge-energies indeterminate, in need of two closure equations, e.g.
ε1 -ε1` = α, ε2 -ε2` = α`.Morphology from the Wulff construction now is a geometrical parallelepiped (Fig. 2b in main text an overlap of two stripes of widths E11` and E22`, whose positions shift with α and α` values.The indeterminacy is again 2, yielding 2 degrees of freedom in the construction translations, while its shape remains unaffected. Third, and sufficiently general to consider, is a material piece enveloped by the three cuts (1, 2 and 3) of infinite crystal.Total there are 6 unknown energies: ε1 and ε1`, ε2 and ε2`, ε3 and ε3`.There are 3 equations from the ribbons, ε1 + ε1` = E11`, ε2 + ε2` = E22`, ε3 + ε3` = E33`, and there is also a triangle with well defined energy, so ε1l1 + ε2l2 + ε3l3 = E123, that is 4 equations for 6 unknowns, requiring 6 -4 =2 closure equations; the indeterminacy being 2. In this case closure can be any pair of linear equations (or nonlinear but one better avoids complications), e. g. ε1 -ε1` = α, and perhaps ε1 + ε2 -ε3 = α`... Cannot be less than 2 (still underdetermined) and cannot be more than 2 making the system overdetermined, cannot be solved.So, again 2 degrees of freedom.
Further, one can prove by inspection that the arbitrary values of these auxiliaries do not change the Wulff construction shape but only its position.For this, at some choice of α, α`, consider a triangle with edges 1, 2, and 3, each at a distance ε1, ε2 and ε3 from the center (an asterisk in the figure), according to the W. c. rules.In principle, changing α, α` affects all ε's values, shifting the position of the edge-lines and possibly truncating the triangle at the corners, turning Wulff construction into irregular quadra-, penta-or hexagon.Assuming for simplicity the opposite edge-lines 2` and 3` to be far enough not to participate, we focus on 1 and 1` shift caused by some increase of auxiliary α.Distance between them being fixed, increase of α makes ε1 greater, shifting the strip 11` down, to the red horizontal lines in the figure, and possibly truncates the initial triangle at its top.However, greater ε1, must also change ε2 and ε3 to generally smaller values, to keep the sum ε1l1 + ε2l2 + ε3l3 constant, and thus shifting lines 2 and 3 to new positions, forming a new triangle (red in the drawing c above).The red triangle must have the same area, as a sum of three bases εi times heights li, line above, and therefore the same total vertical height as the original black triangle.The latter is important, indicating that the top vertex lies under the new 1`-edge (red dashed line), so the triangle is not truncated into trapezia.The change of auxiliary parameter only translates the Wulff construction without changing shape (here, black to red triangle), Q.E.D.
Examples above show that although auxiliary values do affect the edge energy densities, the Wulf construction shape remains unchanged, except translation.With increasing number of sides-edges, the complexity of drawings rises, but rendering the Wulff construction from the edge-energies can be coded, and this is what Figs. 3 and 4 on main text represent: computer-assisted proof by construction, common in geometry.In presence of a symmetry axis (as in C2) obviously the translation is only possible along the axis, as in main text Fig. 3, while in no-symmetry C1 case it occurs in two directions, as in Fig. 4.

Supplementary Section-2b Polygons redundancy
One can imagine many polygons based on set of basic edges, This makes obtaining more equations to resolve the indeterminacy and define the edge energies tempting, but one learns that it does not work, since any added polygon (within the same chosen set of edges/cuts) yields no nontrivial equations, as discussed below.In Supplementary Figure 2a, the edge energy of gray shaded triangle is Egrey = ε1 + ε2 + ε3'.If we add-join the remaining white inverted-triangle to the gray triangle, thus eliminating the edges 3 and 3', then the edge energy of combined gray-white rectangle is: Egrey+white = (ε1 + ε2 + ε3') + (ε1' + ε2 + ε3) -(ε3' + ε3) = (ε1 + ε1') + (ε2 + ε2) = E11'+E22, that is reduced to a combination of 2 ribbons edge energy.Therefore, the new (white) triangle brings no new independent equation.In other words, besides the basic-edge ribbons and the selected (grey) triangle (for "Y" crystal, there're 3 ribbons and 1 triangle), adding other polygons only results in repeating those basic-edge ribbons and single triangle equations, and won't introduce new information.
It is easy to sketch (for brevity) a general proof by mathematical induction.Let's assume that for an N-gon it is true that all perimeter edge energy is a linear combination of the energies of the three basic ribbons and one basic triangle.Then attaching one more basic triangle to it builds a new (N+1)-gon, whose energy obviously consists of N-gon plus the added triangle minus a ribbon (due to eliminated edges at the contact).Therefore, the (N+1)-gon energy is also a linear combination of the energies of the three basic ribbons and one basic triangle.And so on, by induction, any polygon adds no new independent equations, Q.E.D.Alternatively, instead of incremental attachment of basic triangles, any picked "new" polygon can be tessellated into basic triangles (plus/minus energies of the contacts-cuts), to see that such new polygon is redundant, its energy equation is merely a linear combination of already listed 4 equations (1.1-1.4) in the MS, for 3 ribbons and 1 triangle.This makes for indeterminacy 1 (degree of freedom α) and thus requires 1 closure equation.
In the lowest symmetry C1 case, as y-crystal in Supplementary Figure 2b, the analysis is exactly the same except that for 8 unknowns one has 4 basic ribbons plus 2 triangles, that is 6 equations, with no more added by any larger polygons or shapes.

Supplementary Section-2c
Adding cuts/edges expands the basic set, yet the same method works The rank of MS increases but does not change the level of indeterminacy or degrees of freedom α's, which do not affect the shape, so the latter can still be predicted.
One substantial motivation to introduce extra cut(s) can be when empirical knowledge of a particular crystal suggests a specific facet, but which is not among the a priori basic edges.
Let us consider a general no-symmetry (C1) case.A fragment of a polygon from Supplementary Figure 2b (Fig. 1a in the main text) is shown below in Supplementary Figure 3. Including a new cut (edges marked red n and blue n`) adds 2 new unknown variables εn and εn`.In the MS set, one new equation to add obviously follows from a ribbon: εn + εn` = Enn`.One more equation arises from a triangle formed by this cut and two edges from the prior set, for instance ε1 + ε3 + εn` = E13n`, as shown, shaded pink.Any other extra-triangles would not add new nontrivial equations.A geometrical proof, although a bit tedious, is: such extra-triangle also containing n (or n') edge can be presented as a tiling combination (+/-) of the first-chosen triangle 13n` and a few from the prior set.That is by including an additional cut we add 2 unknown edge-energies, but also gain exactly 2 more equations in the MS, still in need of 2 closure equations containing auxiliary α and α` --same 2 degrees of freedom.
Supplementary Figure 3: A fragment of a polygon from Supplementary Figure 2b (Fig. 1a in the main text) including new cut, n`.
Another proof is as follows.Let us assume that for M cuts, that is a polygon with 2M sides-edges, the MS set contains M equations from the ribbons for each cut-direction, and M-2 equations from triangles, so 2 closure-equations make it fully determined (M + M-2 + 2 = 2M equations for 2M unknowns); for a selected fixed α, α`values one has exactly 2M equations and the linear system is fully determined.Now, adding one new "suggested" cut brings 2 variables (its edge energies) and also 2 equations (one ribbon, one triangle).Any other extra-triangle cannot add a nontrivial equation since this would make our system overdetermined, which is unphysical.So, we again conclude that even if the basic edges set is made larger to M → M+1, it merely expands the MS rank while preserving the number of closure equations 2 (or 1, for the C 2 symmetry case).Since we know it is true for small M =1, 2, 3, by mathematical induction it is true for any M, Q.E.D.
An important corollary: extra edges can be added to the polygonal construction, for approximating the edge-energy function ε(a) to any desired accuracy (that is discretizing ε(a) curve, in polar coordinates).This merely makes the algebraic MS larger, yet with the same couple of closure equations, so that solving such a linear system presents no principal difficulty, although is much more taxing and less elegant than using the IA.Importantly is shows that our interpolation anzats IA is not essential for solving the lowsymmetry problem; it is definitely a convenient shortcut in going from small number of basis edges (8 in our approach to no-symmetry case C1) to continuum representation of edge energy for arbitrary directions ε(a), which is customary to use in discussion of Wulff construction.Moreover, an anzats (attempt, "trial answer") is often not expected to have a rigorous proof but is validated by practical use; in our case, all 2D materials considered in literature (graphene, hBN, several TMD, as cited in main text) and in this work produce rather satisfactory agreement with other theoretical predictions (done for higher symmetry cases) of with experimental images shown here.
Supplementary Section-3.Symbol correspondence between schematic and real crystals.Wulff construction requires ε(a) for arbitrary direction (angle a) and for this we invoke an interpolation ansatz (IA): any slanted, vicinal edge is viewed as a sequence of small segments of the basic edges and, accordingly, its energy is decomposed into a sum of the basic energies, in proper proportions [25][26][27] , like c1ε1 + c3ε3, etc.After straightforward trigonometry transformations, one can express the general ε(a) through basic EEs only and the geometry, i.e. the Bravais lattice parameters and angles 25 : Note that at basic edges directions it is exact, while at mid-angles the precision of this interpolation is not essential, usually having no effect on the resultant Wulff shape mostly defined by the minima of this function, which falls on the basic edges.The values of c1 to c6 parameters for the SnSe (and SnS) and AgNO2 are given in Supplementary Tables 1 and 2 respectively.
IA holds as long as the basic edges chosen are lowest in energy.In case if the energy of any other edge due to specific reconstruction or passivation is found to be low, in computation or empirically observed, close or below the reconstructed basic edges and IA, such "new" edge should be added as the basic edge for IA, see also Supplementary Section-2c and the corollary.As we noted in the main text, this ever-possible issue is not specific to low-symmetry indeterminacy that we focus on, but is universal also for any high symmetry exploration.Further, as mentioned above in Supplementary Section-2c, for all 2D materials analyzed in literature (graphene, hBN, several TMD, cited in main text) and in this work, IA produces shapes in rather satisfactory agreement with other theoretical predictions (possible for higher symmetry cases, see hBN in Supplementary Section-10) or with experimental images, shown here.Supplementary Section-6.Master system (MS), its right-hand side (RHS) parameters, and example edge structures.

Supplementary
We estimate the RHS of the master system of equations like (1.1-1.4) by DFT total energy calculations of the 2D shapes, with edges, taken relative to that of the corresponding bulk crystal phase, or the exact same numbers Na, Nb, … of each constituent element-atom, with their appropriate thermodynamic chemical potentials μa, μb, ... .Note that, such an atom-precise subtraction corresponds to a proper Gibbs surface choice in the treatment of continuum-gradual interfaces 38,39 .Subtracting exact same number of atoms, i.e. moles, of each element is also fully in accord with the equimolar cut requirement 5,40 in the atomistic models comparing a cluster with a geometrical volume cut out of a bulk phase, to determine energy excess contribution from the surface (perimeter, in case of 2D).Consider a bi-elemental 2D layer with elements 'a' and 'b', of 1:1 stoichiometry.The energy of the perimeter of any shape of this 2D layer can be expressed as: Eperimeter = Etotal(2D shape) -μaNa -μbNb, where μa and μb are the chemical potentials of 'a' and 'b'.Note that the sum, μab = μa + μb is fixed to the energy of the 'ab'-phase present (3D or 2D form), which is directly computable with standard DFT (or empirical-atomistic potentials) as μab = Ebulk/cell.Remaining degree of freedom is either one of the elemental chemical potentials, or for the convenience of keeping things symmetric, their difference, μ = ½(μa -μb).Dictated by equilibrium with the environment, a certain assumed value of μ can change the RHS and accordingly affect the crystal shape, as we illustrate in case of SnS and SnSe, below and in the main text, Fig. 3.For practical use, it is convenient to rewrite the RHS in the MS, that is the perimeter energy, as Eperimeter = Etotal(shape) -μaNa -μbNb = Etotal(shape) -½μab(Na+Nb) -μ(Na-Nb).In the latter, the first two terms are exactly computable, while the last one reflects the elemental imbalance from exact stoichiometry (Na=Nb) and depends on the environment's chemical potential.All the energies on the RHS are estimated using DFT of the appropriate geometries (i.e.triangles and ribbons) whose perimeter is composed of the edges, while μ can be chosen according to the conditions (e.g.high for Se-rich or low for Sn-rich conditions, in case of SnSe crystal).For example, 1) In the case of ribbons with edges 1 and 1', the master equation is estimated as follows.l1ε1 +l1' ε1' = Etotal(ribbon) -μaNa -μbNb, where l1 = l1' are the lengths of the edges.
2) In the case of a triangle with edges 1, 2 and 3, with lengths l1, l2 and l3 the master equation is as below.
In addition to calculating the energies of ribbons with different edge terminations, we have also considered edge-reconstruction depending on the chemical potential (μ).In doing so, we have constructed several different structures for a given edge.E.g. for Sn rich conditions, the edge has Se vacancies (in varying concentrations), and for Se rich conditions we simulate excess Se-adatoms (in varying concentrations).
Then we utilize only those reconstructed edges that yielded the lowest energies at a given μ.
From Fig. 3 b,c (main text) and Supplementary Figure 9 a,b, we can see that both SnSe and SnS have element-balanced edges as the most stable ones in a wide chemical potential μ range (it's -0.6 to 0.7 eV for SnSe, and -1.16 to 0.44 eV for SnS), where the equilibrium shape is a rhombus.It means that the shape most stable for a wide range of the chemical potential, and hence most probable for these two materials is a rhombus.The element-balanced edges result in constant energies of the ribbons along basic edge directions (Z, A, and A').While under extreme chemical potential conditions (Sn-rich or Se-/S-rich), the element-unbalanced edges reach lower energy in both materials, which introduces chemical potential dependent and competing Z/A/A' edges with lower energies, and hence the rhombus shape is lost.Due to different edge reconstructions, SnSe has a rectangle, a two-corners-truncated rectangle, a non-regular octagon, and a one-corner-truncated rhombus as equilibrium shapes under extreme chemical potential conditions, shown in Fig. 3; while SnS has 2 more shapes than SnSe, a two-corners-truncated rhombus and an octagon, shown in Supplementary Figure 9. Supplementary Section-9.Symmetry classification.Apart from the more popular 2D materials with well defined edge energies, such as graphene, hBN or metal dichalcogenides, there exist several less studied 2D van der Waals materials that lack the sufficient symmetry to define their edge energies, making their shapes theoretically unpredictable.Some such materials are SnS, SnSe, GeAs2, SiP2, AsO2, SnTe, GaTeCl, VOBr2, VOCl2, As2O3, I2O5, In2Se3, LiBH4, YbP5 and AgNO2, all of them predicted to be exfoliable into 2D layers.In particular, SnSe which exists as 2D flakes 39 is especially important, for its ability to exhibit exceptional thermoelectric properties with very high ZT merit in bulk [19][20][21][22][23][24] .Hence, formulating a strategy to predict the shapes of these important materials is not only of interest for fundamental physics, but also for predicting the change in structural and electronic properties of the material that accompany the change in shape.
Here we show that all 2D materials can be grouped into 4 types according to the determinacy of edge energy as Supplementary Table 3, and classify 230 space groups into the 4 types as Supplementary Table 4. 3. The types of 2D materials according to the determinacy of ε.

Supplementary
Figure 2. (a) Y-crystal and (b) y-crystal geometry, with gray shaded triangles.

Table 1 .
The parameters of ε(a) expression for Y-crystal (and SnSe, SnS).

Table 2 .
The parameters of ε(a) expression for y-crystal (and AgNO2).