{'fraction_non_alphanumeric': 0.050398857593066775, 'fraction_numerical': 0.011596415205830215, 'mean_word_length': 3.804471256210078, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Conical Point in the Ferroelectric Six-Vertex Model

We examine the last unexplored regime of the asymmetric six-vertex model: the low-temperature phase of the so-called ferroelectric model. The original publication of the exact solution, by Sutherland, Yang, and Yang, and various derivations and reviews published afterwards, do not contain many details about this regime. We study the exact solution for this model, by numerical and analytical methods. In particular, we examine the behavior of the model in the vicinity of an unusual coexistence point that we call the ``conical'' point. This point corresponds to additional singularities in the free energy that were not discussed in the original solution. We show analytically that in this point many polarizations coexist, and that unusual scaling properties hold in its vicinity.


Introduction
The six-vertex model is one of the few exactly solved models in statistical mechanics. Its origin lies in the problem of counting the number of configurations of the hydrogen atoms in hydrogen-bonded crystals, such as ice (1) and KH 2 PO 4 . (2,3,4) In 1967, Lieb, Yang, and Sutherland (5) obtained exact solutions for various versions of the model, and the most general solution was described by Sutherland, Yang, and Yang. (6) These publications contain virtually no details, and later more explicit accounts were published by Lieb and Wu, (3) and Nolden. (7,8) The paper by Nolden contains a detailed description of the exact solution, and explicit results for some special cases where these can be obtained analytically. However, the low-temperature ferroelectric regime is not treated there.
A more recent application of the six-vertex model is as a model of crystal surfaces. Various restricted solid-on-solid models can be mapped onto the six-vertex model, as was shown by van Beijeren for the body-centered-cubic lattice, (9) and by Jayaprakash and Saam for the face-centered-cubic lattice. (10) This last mapping in particular focuses attention on the low-temperature ferroelectric regime of the six-vertex model. In ref. 10 it was conjectured that there would be a discontinuity of slope at a particular point in the crystal shape. Since in this mapping the crystal shape is, up to a scaling factor, given by the shape of the free energy surface as a function of applied fields, (11) this implies the presence of a singularity in the free energy at this point. However, such a singularity was not mentioned in the general solution of the sixvertex model. (6) Our re-examination of the exact solution in this regime, following the method presented in ref. 7 shows that there is indeed a singularity at this point, and that it has new and unusual properties. In this point, a one-parameter family of polarizations coexists, and near it interesting scaling properties are found.
The organization of the rest of the paper is as follows. In Section 2, we give a short, general discussion of the six-vertex model, and briefly review the expressions from the exact solution that we use further on. Then, in Section 3, we apply those expressions to a particular point in the phase diagram in the low-temperature regime, where a completely analytical solution is possible. The behavior of the free energy and the order parameters at and around this point is discussed. In Section 4 the same is done, in less detail, at the transition between the low-and high-temperature regimes. Conclusions and a summary of the results are given in Section 5. The appendices contain calculational details that were omitted from the main text of the paper. A brief account of this work was published in ref. 12. 2 The ferroelectric six-vertex model A configuration of the six-vertex model is given by a covering of the bonds of a square lattice with arrows, satisfying the ice rule: Every lattice point must have two incoming and two outgoing arrows. There are six possible vertices that satisfy this rule (see Fig. 1). To make the model amenable to a solution by transfer matrix techniques it is necessary to impose periodic boundary conditions in the horizontal direction. As a consequence, the number of down arrows is equal in each row of bonds; this conservation law is essential for the exact solution of the model. The transfer matrix method used in the exact solution is discussed in, e.g., refs. 5, 3, 4. After assigning an energy to each of the vertices, the partition function of the model is simply given by the sum over all allowed configurations of the Boltzmann factors of the configurations. The usual way to write the energies of the vertices is given in Fig. 1. The quantities h and v are to be thought of as fields, acting on the horizontal and vertical arrows respectively, favoring one orientation over the other. (The name "asymmetric six-vertex model" refers to the fact that these fields create an asymmetry between left and right and between up and down arrows.) The fact that vertices 5 and 6 have the same energy is no restriction: since they are sinks and sources of horizontal arrows they have to occur in pairs because of the periodic boundary conditions, so only the sum of their energies is relevant.
The variables that are conjugate to the fields h and v are the polarizations x and y. These are defined as x = 1 − 2f ← and y = 1 − 2f ↓ , where, e.g., f ← is the fraction of the horizontal arrows pointing to the left. x and y are order parameters for the various ordered phases. The free energy per vertex can be viewed as either a function of the fields, F (h, v), or of the polarizations, F (x, y). The two are related by a two-dimensional Legendre transformation or conversely The free energy that is calculated in the exact solution is actually F (h, y), because of the conservation law for the number of down arrows. This amounts to solving the model for fixed values of h and y. The relevant free energies F (x, y) and F (h, v) are then obtained from F (h, y) through one-dimensional Legendre transformations.
Most of the analysis in ref. 7 is valid for any values of the vertex energies. However, certain transformations that are applied do depend on them; more precisely, they depend on the parameter ∆ ≡ 1 2 (η + η −1 − e 2βε ), where η ≡ e βδ , and β = 1/k B T is the inverse temperature. The discussion in ref. 7 is restricted to the case ∆ < 1; we will focus on ∆ ≥ 1. There is a different transformation for ∆ = 1 and for ∆ > 1; both are given in ref. 6. As a consequence, various quantities and functions that occur in the exact solution have different forms in these two cases; they are given in table I (the expressions for the other ranges of ∆ are given in ref. 3,7). ∆ > 1 corresponds to the low-temperature phase of models that are "ferroelectric", i.e. that have ε < δ/2, so that the vertices 1 and 2, with a non-zero net polarization, have lower energies (for h = v = 0) than the vertices 5 and 6 with zero net polarization (we choose δ ≥ 0 without loss of generality). By varying T at fixed δ and ε, a transition to a high-temperature phase (with ∆ < 1) takes place at a critical temperature T = T c (corresponding to ∆ = 1).
We now give a brief review of the main results in Nolden's paper (7) that we will use; for derivations and more details the reader should consult that paper. For fixed ε, δ, and T , the exact solution provides expressions for the free energy, and various other quantities of interest, like the fields and the polarizations, as a function of two parameters, a and b. Different values of a and b then correspond to different points in the phase diagram, i.e. different values of h, v, x, and y. The first step is to find the function R(u) from the integral equation R is a function of u, a real parameter in the interval [−a, a]; it also depends on the values of a and b. 1 This equation results from the consistency relations for the wavenumbers in the Bethe ansatz. Once R is known, the free energy per site can be obtained as the largest eigenvalue of the transfer matrix, with upper or lower sign for the right (R) or left (L) eigenvalue, respectively. Expressions for h and y can be obtained from the "generalized normalization" To find v, x, and the relevant free energies, F (h, y) must be differentiated with respect to h and y, on which variables it depends only implicitly through a and b. All this is described in ref. 7, and that treatment applies to the ferroelectric case without any modification.
For general values of a and b, the only way to solve these equations is numerically; the integral equation can be solved with standard techniques (See ref. 13,Chap. 18), and all quantities can be calculated. The ranges of a and b in table I lead to that part of the phase diagram where y ≥ 0. The other half can be obtained through a simple symmetry operation. An analytic solution is possible for a = 0 and a = −π, and an expansion around these points can also be made. The case a = 0 was examined in ref. 10; it corresponds to second order boundaries between completely and incompletely polarized phases in the phase diagram. We will focus on the case a = −π , which corresponds to the two "conical" points.

Analytic solution for ∆ > 1
In the case that a = −π, the integral equations in the previous section can be solved analytically by means of Fourier series, as was done for a = π in the case ∆ < −1 in ref. 7. It is then also possible to calculate an explicit expression for the free energy, and to examine its behavior around the point a = −π.

Finding R(u)
The first step in the solution is to solve Eq. (2.4). To find the behavior both at and around a = −π, we write the solution R as an expansion in 2 ǫ = π + a: All integrals of the form −a a f (v)dv can be expanded in ǫ as follows: Using this in Eq. (2.4), and substituting the expansion (3.1) everywhere, we get a set of integral equations, one for each order in ǫ, containing R 0 and the δR m . All functions R 0 and δR m are defined on the interval u ∈ [−π, π], and all integrals are over the same interval. Since K and ξ have period 2π, the solutions R 0 and δR m also have this periodicity. (Note that by this periodicity terms with even m in the expansion (3.2) vanish. The terms with odd m reduce to twice the value at v = π.) This periodicity makes it possible to solve for the functions R 0 and δR m by applying a Fourier expansion to the equation. We will now examine the first five orders in ǫ, and solve for R 0 and the first four δR m . The Fourier components of, e.g., R 0 (u) are defined by and similarly for δR m , K and ξ. The Fourier coefficients K n are the same as those given in ref. 7, and ξ n is calculated in Appendix A: This ǫ should not be confused with the energy ε introduced earlier.
Eq. (2.4) to order ǫ 0 is Fourier transforming this gives so, for n = 0, R n is given by For n = 0, (3.6) reduces to 0 = 0, so R 0 is not yet determined. The Fourier series (3.3) for R 0 can be summed, and we find with the upper (lower) signs for b > ν (b < −ν). R 0 will be determined by the next order in ǫ.
To recapitulate, the total expression for R(u) is where δR 2 and δR 4 are real and independent of u, and δR 3 (u) is imaginary.

Calculating y and h
Having determined R, we can now determine y and h using (2.6). We will first calculate the values y 0 and h 0 , the first terms in the expansion in ǫ, and then find terms up to order ǫ 4 from a slightly simplified equation. From table I, p 0 (−π) is given by The second term on the right hand side of (2.6) is The integrals I n ≡ − π −π Θ(−π − v)e −inv dv are calculated in Appendix B. They turn out to be For the aforementioned term we now have with the upper (lower) signs for b > ν (b < −ν). The final result then is So we find, using Eq. (3.11), Now we use Eq. (2.28) from ref. 7 to find the next orders in ǫ. This equation is Substituting the expansions for R, ∂ a R, expanding all occurrences of a to order ǫ 3 , and using Appendix B, we find 0 are real and R ′ 0 imaginary, we have Putting all this together, we now have power series in ǫ for y and h: with the upper (lower) signs for b > ν (b < −ν).

Calculating x and v
To rigorously calculate expressions for x and v, and the Legendre transforms F (x, y) and F (h, v), one has to construct the derivatives of F (h, y) by applying the chain rule. This is what is done in ref. 7. However, by making a few assumptions, it is possible to derive the same expressions from the symmetry properties of the exact solution. The assumptions can be checked against known exact results, in which case they turn out to hold. The first step is to realize that, instead of calculating the free energy as a function of h and y, one might just as well calculate it as a function of another pair of variables, for example v and x. The results will be exactly the same as before, except for the change of variables, and the region of the phase diagram that one finds the solution for. The preceding calculation gives h, y, and F (h, y) expressed in terms of the parameters a and b. The solution holds for y ≥ 0, which corresponds to regions II and III (for b ≥ ν), and region IV (for b ≤ −ν) in Fig. 2. 3 The solution for y is of the form where Y 1 and Y 2 are given functions. If, on the other hand, one calculates F (v, x), one finds exactly the same solution, but now v and x play the roles of h and y, and the parameters are now calledã andb to distinguish them from the other set of parameters. The solution is now for that part of the phase diagram where x ≥ 0, regions V and IV (forb ≥ ν), and region III (forb ≤ −ν). The expression for x is just the same as the one found previously for y, with the same functions Y 1 and Y 2 . We now have expressions for both y and x in regions III and IV, only they are given in terms of two different sets of parameters, a, b andã,b. We need a relation between them to find x and y expressed in a and b. A simple ansatz is to assume that a =ã, and that b andb are related linearly.
If we further assume that the expression for x in region III also holds in region II, as is true for y, we now have expressions for x throughout the region y ≥ 0. The expression for x for b ≥ ν can simply be obtained by taking y for b ≤ −ν and Similarly, x for b ≤ −ν follows from the expression for y for b ≥ ν with the same substitution. The same prescription can be used to obtain v from h. It is possible, to a certain extent, to check the assumptions we made. As a first check, we have calculated x and v to zeroth order in ǫ the rigorous way, obtaining the same results as above. Second, from the symmetry properties we can also derive consistency relations for the expressions for y in various regions. For example, we can repeat the above reasoning for F (−v, −x). This gives the region of the phase diagram where x ≤ 0, regions I and II (forb ≥ ν), and VI (forb ≤ −ν). The expressions for x are now If we make the same assumptions as before,ā = a andb = ν . So to be consistent the expressions for y have to satisfy . Checking the expressions for the first few orders in ǫ (3.25) shows that this is indeed true, giving confidence that the assumptions we made are justified. In addition, the numerical results confirm the expressions for x and v found in this way.

Finding the free energy
Now we are in a position to calculate F (h, y) from (2.5). We need to know the integrals These can be expressed in terms of integrals of the form for n ≥ 0. These are calculated in Appendix C. The details of the calculation of the zeroth order part F 0 of the free energy as a function of y and h can be found in Appendix D. There it is shown that the right eigenvalue is largest for b < −ν and b > φ 0 , and the left eigenvalue is largest for ν < b < φ 0 . Even though at b = φ 0 the two eigenvalues cross, the analytic form of the expression for the free energy F 0 is the same on either side of b = φ 0 . (In fact, this remains true at least up to order ǫ 3 ). In Appendix D we find To find the next two orders of F (h, y) in ǫ we expand Eq. (2.5) as usual, to order ǫ 3 , which gives The coefficients F n can be expressed in terms of the coefficients in the expansions of h, v, x, and y (see Eq. (3.25)): Since the Legendre transforms of the free energy are given by

Discussion
The overall phase diagram as a function of h and v can be obtained numerically; its central region is shown as the base plane in Fig. 3. It contains four phases where the polarizations (x, y) are "frozen in" at their extremal values (±1, ±1). Two of these phases, (1, 1) and (−1, −1), meet along the line h = −v. x and y are discontinuous across this line, so it represents a first order phase transition. In between the four frozen regions are two incompletely polarized phases, where x and y change continuously with h and v. They are separated from the frozen phases by continuous transitions, which have been shown to be of the Pokrovsky-Talapov (PT) type. (10,14) The points where two PT lines meet the first order line are the "conical" points.
To show the relation between the values of the parameters a and b, and the phase diagram in terms of h and v, several lines of constant a and b have been drawn in Fig. 4. This relation is shown only for the part of the phase diagram with b ≥ (ν + φ 0 )/2. This gives one quarter of the phase diagram; the rest of the diagram is then obtained by symmetry. We restricted our computation to this region to avoid certain numerical difficulties that arise near b = ±ν: At b = ±ν, R(u) has a pole and this makes the solution of the integral equations increasingly difficult as b → ±ν. An additional difficulty is that at b = φ 0 , Φ R,L (u) have poles. This does not affect the integral equations themselves, but does make the integrals involved in the calculation of F (h, y), x, and v increasingly difficult as b → φ 0 . We have dealt with this problem by using an integration routine with adaptive stepsize control (See ref. 13,Chap. 4), as opposed to the simple Gaussian quadrature used for solving the integral equations and determining the less-sensitive integrals. However, near b = φ 0 the integration does become quite slow. Fig. 3 also shows the free energy surface F (h, v) superimposed on the phase diagram. F (h, v) is purely linear in h and v in the frozen phases, reflecting the fact that (x, y) = −(∂F/∂h, ∂F/∂v) is fixed. In the incompletely polarized phases, x and y change continuously, so F (h, v) is smoothly curved. Along the PT lines the curved part of F joins smoothly onto the linear part, since the singularity in F (h, v) corresponding to the (second order) PT transition leaves the first derivative continuous. Along the first order transition line there is a jump in slope, so that a ridge is formed. At the two conical points the slope of F depends on the way in which the point is approached, as will be shown later on.
We now turn to the phase diagram in x-y space, Fig. 5. Here the four frozen phases correspond to the corner points (x, y) = (±1, ±1). The incompletely polarized phases take up the rest of the phase diagram outside the lines ℓ + and ℓ − . The two conical points now correspond to the two lines ℓ + and ℓ − , which are given by (x 0 (b), y 0 (b)) for b > ν and b < −ν, respectively (see Eq. (3.25)). By eliminating b in Eq. (3.25) we find that the lines ℓ ± are given by From Eq. (3.36) we see that on, e.g., ℓ + , where a = −π and b > ν, the free energy is F (x, y) = −δ/2 + v 0 y 0 + h 0 x 0 = −δ/2 + (ν/2β)(y 0 − x 0 ). This means that the slope of F (x, y), which is equal to (h, v) = (−ν/2β, ν/2β), is the same everywhere along ℓ + , and that every (x 0 , y 0 ) on ℓ + corresponds to the same point in (h, v), the conical point. So, all the polarizations (x 0 , y 0 ) are stable in the conical point. This is similar to the coexistence of a one-parameter family of magnetizations m with |m| = m 0 (T ) in a zero-field XY model for T < T c . However, there the coexistence is due to an obvious symmetry in the Hamiltonian, while in the six-vertex model the symmetry that allows the polarizations to coexist is not an obvious symmetry in the Hamiltonian and therefore appears to be generated spontaneously. It can be shown that different values of b correspond to entering the conical point with different angles from the incompletely polarized region; if b is kept fixed as a → −π, the conical point is approached from an angle with respect to the line h = −v. Since every different b leads to a different value of (x, y), and (x, y) gives the slope of F (h, v), the slope in the conical point depends on the angle at which it is entered. This gives F (h, v) a geometry similar to the tip of a cone at this point, hence the name "conical" point. The exact solution does not give any states with polarizations in the regions C ± between ℓ + and ℓ − . Thus, no pure equilibrium states with those polarizations exist, and strictly speaking the free energy F (x, y) is not defined in this region. However, a state with an average polarization (x, y) in, e.g., C + can be formed as a mixture of states on ℓ + , properly weighted to give the right average polarization. This leads to a free energy that is also an average of the free energies of the states in the mixture, which means that it is linear in x and y. In this way, the regions C ± are interpreted as describing mixtures of the coexisting phases on ℓ ± , similar to states that are a mixture of coexisting gas and liquid phases at a gas-liquid phase transition. The slope of the free energy surface F (x, y) is then given by (h, v) = (−ν/2β, ν/2β) in C + , and (h, v) = (ν/2β, −ν/2β) in C − . This matches the slopes on the lines ℓ + and ℓ − , so that the first derivative of the free energy F (x, y) is continuous across these lines. Its higher derivatives are not continuous, however, and the free energy has singularities at the lines ℓ ± , since an analytic function cannot be matched to a purely linear function in a regular way. This form of the free energy also means that the fields corresponding to the regions C ± are the same as those for the lines ℓ ± , i.e. the mixtures described by C ± occur in the conical points. The two linear parts of F (x, y) meet along the line x = y, forming a groove. This groove corresponds to the coexistence of the two frozen phases along the first order transition in Fig. 3, and it was already mentioned in ref. 6. However, it is stated there that the groove is the only set of points in x, y ∈ (−1, 1) where F (x, y) is non-analytic, while in fact there are additional non-analyticities along the lines ℓ ± .
To examine the nature of the singularities in the free energy we now specialize to the case b = (ν + φ 0 )/2, so that we are on the lines y = −x and v = −h, moving from the incompletely polarized phase towards the conical point as ǫ → 0. The behavior along this line should be qualitatively the same as along any other line emanating from the conical point; however, it is much easier to deal with just two variables instead of four. Solving for ǫ in δv ≡ v − v 0 = ǫ 3 v 3 + O(ǫ 5 ), and substituting this in the expression for y (Eq. 3.25), we find that So, as δv goes to zero, y approaches its value of y 0 in the conical point with a slope that diverges like δv −1/3 . This value is in agreement with recent work (15,16) which shows that certain (1+1)-dimensional surface growth models map onto the conical point in the (h, v) phase diagram, and that Kardar-Parisi-Zhang (KPZ) (17) scaling holds. In particular, the exponent of 2/3 in Eq. (3.40) is given by 1/z, (16) where z = 3/2 is the dynamic exponent for the KPZ universality class. Relation (3.40) can be integrated to give F (v), since ∂F (v)/∂v = −2y (the factor 2 comes from the fact that x = −y and h = −v). This gives The first few terms of this expression reproduce the expansion of F (v) given in terms of ǫ earlier (Eq. 3.37). Similarly, an expression for F (y) can be derived by writing v in terms of δy ≡ y − y 0 and integrating the relation ∂F (y)/∂y = 2v. It is Again, this confirms the expansion given earlier (Eq. 3.36). The appearance of the terms δv 5/3 and δy 5/2 means that both F (x, y) and F (h, v) have non-analyticities at the first-order transition, in addition to the expected discontinuities in their derivatives. Thus, the traditional mean-field picture of a first-order phase transition as two branches of the free energy simply crossing is not valid. As the transition is approached from the incompletely polarized phase, the free energy becomes singular, and there is no obvious way to extend this branch of the free energy into the regions C ± in Fig. 5 to find metastable states. It is thought that the existence of non-analyticities at coexistence boundaries is a generic feature; (18) the two-dimensional Ising model, e.g., shows a weak essential singularity in F (h) at h = 0 and T < T c . Vector spin models (where, as here, a whole continuous family of order parameter values coexists at a single point) show power law singularities similar to the ones found here, albeit with different values for the exponents (we are not aware of any previous observation of, e.g., a "susceptibility" ∂ 2 F (v)/∂v 2 with an exponent of −1/3). As a further comparison we mention that at the PT transition, where y goes to 1, the free energies behave like (10) F (y) = A + B(1 − y) + C(1 − y) 3 + · · · F (v) = α + βδv 3/2 + γδv 2 + · · · (3.43) Note that the conical points are reached either by letting a → −π or by letting b → ±∞ (or some combination of both), as was remarked earlier. Taking a → −π brings us into the point along some nonzero angle with respect to the PT line in (h, v) space, whereas b → ±∞ takes us in tangent to the PT line. The scaling properties derived above [Eqs. (3.40-3.42)] correspond to the former case. In the special case of a trajectory tangent to the PT line, the scaling would be determined by the limit b → ±∞. In particular, we know that the free energy F (h, v) is analytic in h and v along the PT line itself (b → ∞ with a = 0). By coming into the conical point along other curves which are tangent to the PT line, intermediate scaling forms of various kinds could be obtained, but such scaling behavior is not generic and therefore somewhat contrived.
This brings us to a closely-related point which is: In which region is the behavior near the conical point dominated by the proximity to the conical point, and in which region is it dominated by the proximity to the PT lines? For example, we might ask how to determine the region over which we can observe the scaling (e.g. the Pokrovsky-Talapov behavior) associated with the boundary of the completely polarized region. The preceding discussion suggests that this region is very narrow. Indeed, by considering when the second term in the denominator of (A19) of ref. 10 is of order 1, we find that a curve of the form d f ∼ d 2 c separates the two regions. Here, d c is the distance in the h-v plane to the conical point and d f is the perpendicular distance to the PT line. Note that this curve is tangent to the PT line, deviating only quadratically as one moves away from the conical point. Application of (A19) in ref. 10 also tells us that the coefficients of the (1 − y) and (1 − y) 3 terms of the expansion of the free energy about the completely polarized region [as in Eq. (3.43)] are unchanged by the presence of the conical point, and thus that the coefficient of the PT scaling term is not modified. However, higher-order terms in the expansion diverge as one approaches the conical point, representing the fact the the PT scaling will be observable over an ever-smaller angular region near the PT line.
Finally, we might ask whether there are any analytic approximations we can make which are valid over the entire phase diagram. Indeed, we can perform an expansion of Eqs. (2.4-2.6) to lowest nontrivial order in the temperature. (19) To do this, we ignore terms in ξ(u), K(u − v), p 0 (u), and Φ L,R (u) which are down by factors of order e −2ν (or higher powers) from the dominant terms. Since to this order the kernel K(u − v) is identically 1, the integral equation can then be solved analytically and we can write down explicit expressions for h, v, x, and y in terms of a and b. The expressions for h and v, however, involve an integral which apparently must be evaluated numerically. Furthermore, even the expressions for x and y cannot be inverted to eliminate the parameters a and b from the equations entirely (leaving us, as we had hoped, with direct relations between h, v, x, and y). Therefore, these expressions, while providing us with a simpler numerical task than solving integral equations, do not lead to any profound simplification.
Since we have ignored terms down by orders of e −2ν , the low-temperature approximation is quite accurate even up to a large fraction of the critical temperature. One particularly interesting fact about the accuracy of this approximation should be noted: Although the general features, such as the locations of the PT lines (corresponding to a → 0) are not obtained exactly, everything involving the conical points (i.e., the limit a → −π and the limit b → ±∞) is exactly correct. This includes the location of these points in the (h, v) plane, the expressions for l ± , and even the expansions for x, y, h, and v about the conical point [Eqs. (3.24-3.25)] to the orders we know them. (In particular, this means that the transition temperature T c is obtained exactly-a curious fact for the first term of a low-temperature expansion!) We hypothesize that the fact that the behavior at the conical points is obtained exactly could be evidence that this behavior is controlled by a zero-temperature fixed point.

Analytic results for ∆ = 1
The case ∆ = 1 corresponds to the critical point T = T c , where the first-order line in Fig. 3 has shrunk to zero length (see Fig. 6), and the two conical points at the end of this line have now coalesced into one point, which will disappear for ∆ < 1.
Here, it is possible to obtain some analytic results for a → −∞, which corresponds to this limiting point. It turns out to be much harder to make an expansion about a = −∞ than about a = −π, and we were only able to obtain a limited amount of information about this case.
We can go through the same steps as in section 3, using the relevant expressions from table I. Now equation (2.4) can be solved by means of a Fourier transform, defined as The Fourier transform of K is again the same as in ref. 7, and the one for ξ can easily be calculated: As usual, R 0 (0) is not determined. However, that does not prevent us from Fourier transforming back to find R 0 (u), It is, to a certain extent, possible to expand R(u) in 1/a around a = −∞. It can be shown that the leading corrections are of the form with f and g real functions, f even and g odd. These two functions are determined by two integral equations which still contain a. Even though numerical evidence shows that f and g have well-defined limits for a → −∞, there is no obvious way to take this limit in the integral equations, and we have not been able to obtain closed-form expressions for the two functions.
To get y 0 and h 0 we again apply (2.6). It turns out that we have to be slightly careful in taking the limit a → −∞, because we have to consider the possibility that b diverges along with a, so we cannot just set a = −∞ as we did before. From table I it is apparent that p 0 (a) → 0 no matter how a and b go to infinity. For the integral term we find with the upper (lower) signs for b > 1/2 (b < −1/2). So we find that for a → −∞ and h 0 = 0 in both cases. So by letting b go to infinity along with a, any positive value of y 0 can be selected. It follows from the expansion (4.5) that the corrections to both h and y will be of order 1/a; however, the coefficients of these terms cannot be calculated without knowledge of the functions f and g.
To find expressions for x and v we can again use the symmetries discussed in section 3.3, if every occurrence of ν is replaced by 1/2. In particular, we immediately find that v 0 = 0.
It is also possible to find the free energy in the limit that a → −∞. As before, it turns out that for b < −1/2 and b > φ 0 the right eigenvalue is largest, and for 1/2 < b < φ 0 the left eigenvalue is largest. In all cases we find that F 0 (h, y) = −δ/2, independent of b. This confirms the above conclusion that v 0 = ∂F 0 /∂y 0 = 0.
We can again examine the scaling behavior of the free energies. As in the case ∆ > 1, F (h, y) is a regular power series in the expansion variable, 1/a. However, in contrast to that case, the expansions for y and h both contain linear terms in the expansion variable. This implies that the polarizations are non-singular functions of the fields. Doing the Legendre transform of F (h, y) to obtain F (x, y), we find that the coexistence regions C ± have shrunk to the line x = y, and the behavior around this line is now quadratic in x − y. So both the groove at x = y and the singularities at the lines ℓ ± (which are now collapsed to x = y) have disappeared. The free energy can thus be written where the coefficient A = A(x+y) is a function of the position along x = y. Similarly, the free energy as a function of the fields is given by where y 0 depends on the ratio a/b according to Eq. (4.7), i.e., it depends on the details of how the point h = v = 0 is approached. Note that both the jump in the slope of F (h, v) on entering this point, and the singularities on approaching it, have disappeared; however, the slope in the h + v direction (given by −y 0 ) still depends on the way in which this point is approached. This reflects the fact that, as T → T c from above, the two PT boundaries approach each other, and the curvature of F (h, v) in the h + v direction goes to infinity.

Conclusion
To summarize, we have studied the phase diagram of the asymmetric, ferroelectric six-vertex model in the low-temperature phase. From numerical results we get a global picture of the behavior of the free energy as a function of the polarizations and fields, F (x, y) and F (h, v), and by analytical methods we found several new and interesting features. For T < T c (∆ > 1), in addition to the groove that was already known to exist in the free energy surface F (x, y), there are additional singularities along the lines ℓ ± , corresponding to two conical points in F (h, v). In these points, which lie at the ends of a ridge in F (h, v), all polarizations on ℓ ± are stable for one value of the fields h and v. Also associated with these points are the two coexistence regions C ± , where mixtures of the coexisting states are stable.
As T → T c (∆ → 1), the size of the regions C ± decreases, as does the length of the ridge in F (h, v). At T = T c , the two conical points merge as the ridge is reduced to a point, and the regions C ± vanish. There is still a family of polarizations, x = y, that is stable in the (one) conical point. If there are any non-analyticities left in the free energy they are of higher order than the terms we have been able to study. Above T c , the conical point disappears completely.
The behavior of the free energies when approaching the conical points is rather unexpected. Both the fact that there are power law divergences, while the phase transition is discontinuous and thus first-order, and the specific exponents of the power laws, are interesting. It would be useful to find metastable states in the regions C ± with a uniform polarization, and examine how such a state decays into a mixture of stable states on the lines ℓ ± . However, because of the singularities in F there is no obvious way to analytically continue the free energy into the unstable region, in order to find such metastable states in a straightforward manner. Also, there does not seem to be any mechanism to select one of the many possible mixtures that such a state might break up into. Thus one can ask how the process of phase separation, or spinodal decomposition, would actually take place.
A particular system in which this issue needs to be addressed is the phase separation of crystal surfaces. (20,21) Through the mapping of the six-vertex model onto restricted solid-on-solid models (9,10) the polarizations x and y are found to be equivalent to surface orientations, and the free energy F (h, v) turns out to be a replica of the crystal shape itself, (10,11) while F (x, y) is the surface tension. Thus, the conical points represent actual points on the crystal surface that have a conical geometry. 4 In these points, various surface orientations all coexist. Also, there is a discontinuity in orientation as one crosses a conical point, since the orientations that correspond to states in the regions C ± are not stable anywhere. Furthermore, if a crystal is cleaved to form a surface having one of these unstable orientations (say in C + ), it should decompose into a mixture of stable orientations selected from the set ℓ + , which then coarsens over time. This is the phase separation process mentioned in the previous paragraph. These issues are discussed in more detail in ref. 21.
So for n ≥ 0 For n < 0, we substitute z = e −iu instead, and in a similar way we find

Appendix B The integrals I n
The integrals I n are defined as The function Θ(−π − v) can be rewritten as where the ln has a branch cut along [0, ∞), so that its imaginary part runs from 0 to 2πi. This is because Θ(−π − v) must run from 0 for v = −π to −2π for v = π.
For n ≤ 0 we use (B.2), define z = −e iv , so that with A = e −2ν , B = e 2ν , C = 2ν. The contour is given in Fig. 7. So we must examine integrals of the type for n ≥ 0. For n = 0 where now the ln's are the principal branch (i.e. with a branch cut along (−∞, 0]), the γ's are the imaginary terms that result from approaching the cut from above or below. Since the integral around the little loop goes to zero, and the line integrals give (1 − A n )(γ 1 − γ 2 )/n, we find For n = 0 there is an extra pole at z = 0, so we must replace the right hand side of (B.5) with The two line integrals give −(γ 1 − γ 2 ) ln A = 2πi ln A. So now For n < 0 a similar calculation is performed by defining z = −e −iv . The end result for I n is Substituting this in Eq. (2.5) gives Now, defining z = e iu , and using Appendix C, Here we find that − βF R (h, y) = 1 2 (ln η + 2βh) + ν R 0 (D.10) Now we have to collect all information on −βF R,L (h, y) and find the maximum over R, L for the various intervals. The simplest case is b < −ν: Here, − βF R (h, y) = 1 2 (ln η + 2βh) + ν R 0 = βF L (h, y) (D.12) Also, R 0 = − e ν+b 1 + e ν+b and 2βh = ν (D.13) We also know that ln η > ν, say ln η = ν + ζ with ζ ≥ 0. Then So −βF R (h, y) corresponds to the larger eigenvalue in this regime. Then, For b > ν similar reasoning shows that for b ∈ [ν, φ 0 ] the left eigenvalue is largest, and for b > φ 0 the right eigenvalue is largest, and we find Table I. The definitions of various quantities and functions for the two transformations for ∆ ≥ 1. Note that some of the functions depend on b, but the explicit dependence has been suppressed for brevity. We use the definitions of Θ and K as they are given in ref. 7; these definitions differ from the ones used in refs. 3, 5, 6 by a sign in Θ and a factor of 2π in K.      13. Note that the phase boundaries themselves correspond to the limit a → 0, with the value of b determining the location along these boundaries. The conical point corresponds to either a → −π or b → ∞. In the former limit, the angle of approach to the conical point is determined by the value of b, while in the latter limit, the approach is tangent to the phase boundary for all values of a.  Fig. 3. The two coexistence regions C ± correspond to the two conical points (h, v) = ±(− ν 2β , ν 2β ) at each end of the first-order line in Fig. 3. The regions are bounded by the lines l ± and y = x (shown as bold lines). The line y = x is also a coexistence region, corresponding to the entire first-order line in Fig. 3. [Although it is not obvious from this perspective, F (x, y) is constant along y = x.]  A branch cut runs along the real axis between A and B. The labels 0 and 2π above and below this branch cut give the value of the imaginary part of the logarithm. For the case n = 0, there is also a pole at z = 0. Fig. 8. Contour of integration in the complex z-plane for the integrals J n of Eq. (C.1) for the case B < A < 1. A branch cut runs along the real axis between A and B. The labels ±π above and below this branch cut give the value of the imaginary part of the logarithm. For the case n = 0, there is also a pole at z = 0.