Commensurate-incommensurate transition in the chiral Ashkin-Teller model

We investigate the classical chiral Ashkin-Teller model on a square lattice with the corner transfer matrix renormalisation group (CTMRG) algorithm. We show that the melting of the period-4 phase in the presence of a chiral perturbation takes different forms depending on the coefficient of the four-spin term in the Ashkin-Teller model. Close to the clock limit of two decoupled Ising models, the system undergoes a two-step commensurate-incommensurate transition as soon as the chirality is introduced, with an intermediate critical floating phase bounded by a Kosterlitz-Thouless transition at high temperature and a Pokrovsky-Talapov transition at low temperature. By contrast, close to the four-states Potts model, we argue for the existence of a unique commensurate-incommensurate transition that belongs to the chiral universality class, and for the presence of a Lifshitz point where the ordered, disordered and floating phases meet. Finally, we map the whole phase diagram, which turns out to be in qualitative agreement with the 40 year old prediction by Huse and Fisher.


I. INTRODUCTION
Recent experiments on Rydberg atoms [1,2] have brought back the problem of commensurateincommensurate (C-IC) transitions, and of the nature of the melting of an ordered period-p phase into an incommensurate one. This transition was originally discussed in two dimensional classical systems in the context of the melting of an adsorbed layer on a substrate lattice, and simultaneously introduced by Huse [3] and Ostlund [4] as a series of chiral models exhibiting p − 1 types of domain walls, such that domain walls between ordered domains A | B and B | A have different energy, introducing a chiral into the problem. For p = 2, the system is commensurate even in the disordered phase, and there is no commensurate-incommensurate transition. The transition, if continuous, generically belongs to the Ising universality class. For p ≥ 5 the transition is always a two-step one, with a critical phase in between. In that phase, the configurations are made of permutations of stripes of the p different orders with domain walls at an average distance of l ∝ 1/q, where q is the incommensurate wave-vector of the correlations. As the temperature approaches the ordered phase, the domain walls repel each other and the wave-vector converges to a commensurate value with an exponentβ (q − q 0 ∼ tβ, where t = (T − T c )/T c is the reduced temperature). It is commonly agreed that the transition is in the Pokrovsky-Talapov (PT) [5] universality class, characterised by critical exponents : ν x = 1/2, ν y = 1,β = 1/2, where ν x/y describes the divergence of the correlation length in the directions x and y respectively. The high temperature transition between the critical and the disordered phases is a Kosterlitz-Thouless [6] one driven by the unbinding of the dislocations [7] which become relevant at η = 1/4, where η is the exponent of the decay of the spin-spin correlation function.
For the p = 3, 4 cases, the picture might be different. It was first argued by Huse and Fisher [8] that a floating critical phase bounded by two transitions might not appear right away when introducing the chirality. Instead the commensurate-incommensurate transition could be unique and would belong to a new chiral universality class characterised by the product of the correlation length ξ and the wave-vector q − q 0 converging to a constant at criticality (β = ν), and by a dynamical (or anisotropy) exponent z ≡ ν y /ν x different from 1 (ν y = ν x ).
By contrast to the p = 3 case, the p = 4 literature is less rich. We note however that recently an investigation of a quantum version of the model has suggested the presence of a chiral universality class as well [34]. In this article, we address the p = 4 case by studying the classical version of the chiral Ashkin-Teller model, and we show that both possibilities (chiral transition or two-step transition) are realized depending on how close one is to the Potts or clock limit of the model. The paper is organised as follows. In Section II, we describe the model and the different possibilities for the phase diagram. In Section III, we describe the methodology and the CTMRG algorithm. In Section IV, we benchmark the algorithm and discuss its power and limitation. In Section V we present our results. We compare and study in detail the phase diagram of the chiral clock-Ising model, where a floating phase opens right away, and of the four-state chiral Potts model, where we observe a chiral transition, and we propose a full phase diagram of the model.

II. CHIRAL ASHKIN-TELLER MODEL
The 2D Ashkin-Teller model (AT) [35] is defined with two Ising spins τ, σ ∈ {±1} on each site with energy: where the sum runs over pairs of nearest neighbors. It has been shown [36] that there is an exact critical line known from duality where the correlation length exponent varies continuously from the four-state Potts model (λ = 1) to the clock model (λ = 0) equivalent to two decoupled Ising models as: To our knowledge there is no known theoretical value known forβ, but it is believed to be larger than 1, implying that the product ξq → 0 along this line [9]. The chiral Ashkin-Teller model can be defined by adding a chiral perturbation along the x-direction: This model was originally introduced by Schulz [11], who showed that the chiral perturbation is irrelevant for λ > λ c1 ∼ 0.9779. For large enough values of ∆, the C-IC transition is expected to be a two-step transition separated by a floating phase. For λ < λ c1 , for small enough values of ∆ two possibilities arise, either the transition is two-step and the floating phase opens right away or the transition is unique and belong to the chiral universality class with ν x =β. The two different scenarios are depicted in Fig. 1. On the other hand, for λ > λ c1 , the perturbation is irrelevant, and for small enough values of ∆, the transition is expected to be in the AT universality class. The parameter range between λ c1 and the Potts limit is very small, actually too small to be observed numerically as we shall see. For larger ∆, beyond this small range of AT transition, the question of a unique chiral transition or a floating phase remains unchanged.
For the rest of the paper, we will refer to the λ = 1 case as the four-state chiral Potts model and to the λ = 0 case as the chiral clock-Ising model in reference to their respective universality classes at ∆ = 0.

III. METHODOLOGY
We use the same methodology already developed in the previous work on the three-state chiral Potts model [37] where more details are available.

A. CTMRG
The CTMRG algorithm was first introduced by Okunishi and Nishino [38]. It is a numerical method which combines Baxter's corner matrices [39] and Steve White's renormalisation group density matrix method [40]. It is most commonly used for two dimensional quantum systems as a contraction algorithm of wave functions [41], but has recently shown promising results on classical systems as well [42][43][44].
The investigation of classical systems with tensor networks is done by expressing the partition function in the thermodynamic limit as an infinite tensor network and contracting it. CTMRG sets a way to contract such an infinite square tensor network made of local tensors a on each vertex as shown in Fig. 2. In particular, it allows one to express the Gibbs measure observables as a contraction made of eight tensors denoted by T = {C 1 , T 1 , C 2 , T 2 , C 3 , T 3 , C 4 , T 4 }, with corner tensors C i of dimension χ × χ and row/ column tensors T i of dimensions χ×4×χ. The parameter χ is what controls the numerical approximation and is usually referred to as the bond dimension. The local tensors a can be expressed in many different ways, we choose the most common one given by with Q x/y ij the Boltzmann weight matrices between two spins on the horizontal/vertical axis. In the (σ, τ ) = {(1, 1), (1, −1), (−1, 1), (−1, −1)} basis, Q x is given by: e x0 e x1 e x2 e x3 e x2 e x0 e x3 e x1 e x1 e x3 e x0 e x2 e x3 e x2 e x1 e x0    with Q y is defined similarly with ∆ = 0. CTMRG therefore provides an easy way to measure directly local observable and two points functions. The algorithm can be summarised in two steps [45]: (i) Update: each corner matrix is extended by adding a row, a column and a local a tensor. And each row/ column tensor is extended by adding a local tensor a. We illustrate this step for the tensors C 1 and T 1 in Fig.  3. The other tensors are updated similarly.
(ii)Truncation: the update increases the bond dimension of the environment tensors by a factor 4. In order to avoid an exponentially growing bond dimension one needs to project the tensors into a subspace. The isometries are given by doing the singular value decomposition on some density matrices and keeping only the χ largest singular values. The choice of the density matrices influences greatly the convergence of the algorithm. We choose the ones originally proposed by Nishino and Oknishi and given by: with C i the extended corner matrices. A full iteration for one corner and one column tensor is illustrated in Fig. 3. Both steps are repeated until the energy between two iterations converges to some threshold, at which point we consider the thermodynamic limit reached. We note that the initial tensors represent the boundary conditions on the infinite systems and influence the convergence of the algorithm. We choose to set the boundary conditions to be open as it gives the best convergence.

B. Effective exponents and critical temperatures
Our analysis is based on effective exponents and some general hypotheses regarding the nature of the transitions. The effective exponent of some physical quantity A which behaves algebraically A ∼ t θ is defined as: with t =| T − T c |, the reduced temperature. Due to crossovers and various corrections, the effective exponent might not be equal to the critical exponent if evaluated outside the critical regime but should converge to it as t → 0. If the C-IC transition is unique, effective exponents, and in particular ν eff from both sides of the transition, should converge to the same value at t = 0. It turns out that this requirement is enough to fix T c . On the other hand, if the transition is two-step, due to the exponential divergence of the correlation length from the KT transition at high temperatures, asking for a unique limit of ν eff from both sides of the transition will lead to ν > 1. Such an exponent is unphysical and one can conclude that the hypothesis is wrong. Thus, the transition must be two-step. In that case, fixing T KT by looking for the regime where η = 1/4 turned out to be more accurate than fitting the correlation with an exponential, where η is the exponent of the decay of the spin-spin correlation function: T P T is fixed by the condition lim t→0 ν y eff = 1 from the ordered phase. From now on, we will use ν HT and ν LT to denote the limit of the effective exponent from the respective high and low temperature regimes.

C. Extrapolation of the correlation length and wave-vector
The correlation length and wave-vector can be easily computed through partial diagonalisation of transfer matrices. If one denotes the normalised ordered eigenvalues of the transfer matrix with λ i = j e iφj , j = 1, 2, . . ., then the correlation length and wave-vector are given by: This estimation, which corresponds to a given finite bond dimension, can be improved using extrapolations [46] with respect to gaps (δ) within the transfer matrix spectrum. Rams et al suggested that the inverse correlation length scales linearly such that:  Various gaps δ can be used. In particular, in the disordered phase we systematically used δ = 4 − 2 . In some cases the energy levels of the transfer matrices cross and in order to follow a given gap closing, we extrapolate using different eigenvalues. An example of such a case is given in Fig. 4. As for the wave-vector extrapolation, we used differences in the argument of the eigenvalues, usually we used δ = φ 4 − φ 2 .

D. Error bars
The errors on the effective exponents come from the uncertainty on the critical temperature δT c and the error in the correlation length from the linear fit δξ. A first order Taylor expansion gives: Error bars on the critical exponents are computed by fitting error bars of the effective exponents as shown in

Ising point
At the λ = ∆ = 0 point, the model decouples into two independent Ising models. Therefore, the transition belongs to the Ising universality class characterised by an exponent ν = 1 and a specific heat that diverges logarithmically, i.e. an exponent α = 0 with a multiplicative logarithmic correction. The critical temperature is known from duality and is given by T c = 2/ log( √ 2 + 1). The numerical results are summarised in Fig. 5. A linear extrapolation gives ν HT = 0.999 ± 0.01 and ν LT = 0.992 ± 0.008 in excellent agreement with the theoretical value ν = 1. The wave-vector exponent extrapolates tō β = 1.262 ± 0.004, implying that the product ξq → 0.

Four-state Potts point
We now turn to the ∆ = 0 and λ = 1 point, which belongs to the four-State Potts universality class. It is characterised by the exponents ν = 2/3 and α = 2/3 with multiplicative logarithm corrections [47] both in the specific heat and in the correlation length. The critical temperature T c = 4/ log 3 is also known from duality .
Due to these corrections, the measurement of ν is known to be difficult. In particular, previous studies with the Monte Carlo Renormalization group method [48]  1 What about > 1?
We know that the there are two Ising transitions at = 0. We sho  1 What about > 1?
We know that the there are two Ising transitions at = 0. We show t  1 What about > 1?
We know that the there are two Ising transitions at = 0. We  sults. A linear extrapolations of the effective exponent on Fig.6 gives ν = 0.713. This is quite far from the exact value. Indeed, due to the logarithmic correction, extrapolating the results with a linear fit is wrong and one expects the effective exponent to converge to 2/3 with an infinite slope. We are not close enough to the critical temperature to observe this change of behaviour. For the incommensurability, we observeβ = 1.71±0.02. Due to the poor accuracy of ν, we can ask how much one can trust this result. However, we notice that the product ξq has exponentβ−ν = 0.99±0.06. Thus, if the product cancels the corrections in the wave-vector and in the correlation length, one would haveβ = 2/3 + 0.99 = 1.66 ± 0.06 which we note to be close to the three-state Potts critical exponent. Results are summarised in Fig. 6. In any case, we established that, as expectedβ > ν and ξq → 0.
We now look at the C-IC transitions. In particular, we study in detail the phase diagram of the chiral Ising and the four-state chiral Potts models. For both systems, we observe a KT and a PT transition enclosing a floating phase. In addition, for the four-state chiral Potts model, we found a unique chiral transition.

A. Pokrovsky-Talapov transition
We first consider the λ = 0 model at which we observe a PT transition for ∆ as small as 0.05. As we will see later, in contrast to larger values of λ, the floating phase extends over a broad parameter range for relatively small values of ∆ and we are able to locate T KT using the criterion η = 1/4 for ∆ > 0.1. Looking at the phase diagram of Fig. 7, the KT and the PT transitions do not seem to cross before the Ising point and we conclude that, as expected for this model, the floating phase opens up as soon as the chirality is introduced.
In Fig.8, we present a more careful study of the transition at ∆ = 0.4. We have plotted the inverse correlation length in the y-direction and the square of the wave vector with respect to the temperature. We were not able to measure ξ x from the ordered phase due to poor extrapolations. Because in the floating phase the algorithm converges very slowly, the wave-vector values shown have not been extrapolated but were made with fixed bond dimension χ = 200 and χ = 150. One first notice the accuracy of both linear fits on a large temperature interval, in agreement with critical exponent ν y = 1 andβ = 1/2 with no or very little corrections. If these are indeed the correct critical exponents, the critical temperature can be determined by the intersection with the temperature axis either using the inverse correlation lenght or with the squared wave-vector. This gives respective values of T P T = 1.806 and T P T = 1.807, in agreement with each other up to 10 −3 , and provides a self consistent picture in agreement with a PT transition. Furthermore, in Fig. 9, we observe that the specific heat does not diverge from the ordered phase but has a clear divergence from the floating phase, as expected at a PT transition.
For λ = 1, we observe a PT transition only at ∆ > 0.7 with at first, a narrow floating phase. We illustrate the results at ∆ = 0.8 in Fig. 11. For λ = 0 and ∆ = 0.4, we deduced T P T from the critical exponents. Here, we work the other way around, we fix T P T and deduce the exponents. Setting T P T with ν LT y = 1 at the transition fixesβ = 0.53 ± 0.04 and ν LT x = 0.51 ± 0.05, once again providing a self consistent picture with exponents in good agreement with the PT universality class. In general, the larger the floating phase, the slower the exponent η varies inside it. Therefore, we have best estimates of the critical temperature for narrow floating phases rather than larger ones.

B. Kosterlitz-Thouless transition
The main limitation of the algorithm is its convergence. It turns out that close to the PT transition, the algorithm does not converge to sufficiently low precision to give sensible results. This forbids us to use the η = 1/4 criteria to determine T KT when the floating phase is too narrow. In that case, an exponential fit of the correlation length will not give enough accuracy to distinguish the two critical temperatures.
As already mentioned, for λ = 0, the floating phase extends over a large parameter range so that even at ∆ = 0.1 the η = 1/4 criterion can be used to locate the KT transition. We now illustrate such log-log fit with different temperatures in Fig. 10. We notice that the log-log plot is linear only above some threshold distance that increases as the temperature decreases. This is due to the fact that the true power law decay is only visible above the average distance between domain walls l = 2π 4q which diverges as one approaches the PT transition. For T = 2.22, we evaluate log l 3.99 which is in good quantitative agreement of where the curve starts to be linear in Fig. 10. This makes the study of the power law decay close to the transition extremely difficult.
For λ = 1, we can use the η = 1/4 criterion to determine T KT only above ∆ > 1. In contrast, for 0.7 < ∆ < 1 the floating phase is too narrow to distinguish the two critical temperature transitions. In any case we can identify the transition to be in the KT universality class from the divergence of ν x and ν y from the high temperature regime at the transition as shown in Fig. 11.

C. Chiral transition and Lifshitz point
We now move to λ = 1 and small values of ∆. We recall that at this point, the chiral perturbation is irrelevant and the transition should initially be a unique one in the AT universality class. Yet, we do not see any incommensurate melting in the AT universality class. This is probably due to the fact that the parameter range of the AT transition is too small to be detected by our approach. The results are summarised in Fig.12 and 13. Based on those, we exclude the possibility of a two-step transition. Indeed, fixing the critical temperature up to 10 −4 with a unique exponent ν y fixes a unique exponent ν x with good accuracy. We note that both exponents are smaller than one and the transition must then be unique.  Fig. 13. This is probably due to strong crossovers. However, in contrast to the Potts point we extrapolateβ = 0.70 ± 0.02. Such drastic change of behaviour allows us to conclude that the nature of the transition must be different. Once again, we note that ν HT x =β within 1.4% and in agreement with a chiral transition.
As already mentioned, just beyond the Lifshitz point, we will obtain ν y > 1 if we try to locate a unique transition due to the KT nature of the upper transition. We can therefore locate the Lifshitz point with the condition ν y = 1, leading to a Lifshitz point at ∆ L = 0.705 ± 0.01.

D. Phase Diagram
By applying the same methodology systematically we have been able to map the whole phase diagram, and to identify a line of Lifshitz points. This leads to three different regimes depending on the value of λ. For 0 ≤ λ < λ c2 , the chiral perturbation opens right away a floating phase. By contrast, for λ c2 < λ < λ c1 , the chiral perturbation leads to a unique transition in the chiral universality class, and a floating phase only opens for a non-zero value of ∆. Finally, for λ c1 < λ ≤ 1, the transition should first be in the AT universality class, a regime too small to be detected numerically. For small but-non zero value of ∆, the transition should become a unique transition in the chiral universality class, and finally, for larger ∆, a floating phase should open. These results are summarised in Fig.14, where we show the three dimensional phase diagram, and in Fig.15, where we show a projection of the phase diagram on the (∆, λ) plane. The similarity between the phase diagram of Fig.14 and Fig. 12 of Ref. [9] is striking. The transition is clearly unique for a range of λ below the Potts point λ = 1, and a two-step one for a range of λ above the clock-Ising point λ = 0. A linear extrapolation of the Lifshitz point at λ = 0.6 and λ = 0.7 towards ∆ = 0 leads to the estimate λ c2 0.42, but the curvature of the line of Lifshitz points in Fig. 15 suggests that this value should be considered as an upper bound for λ c2 . This result for λ c2 is in agreement with recent DMRG simulations on the quantum version of the model [34], for which a chiral transition has been identified in a model with λ ∼ 0.57.

E. Specific heat and Kibble-Zurek exponent
We now discuss and analyse the specific heat exponent α obtained along the line λ = 1 and its relation with the Kibble-Zurek exponent. We recall the hyper-scaling relation ν x + ν y = 2 − α. If such a relation applies, there are two different ways to determine if, as reported earlier for the three-state Potts model [37], α is constant along the transition: Either through the direct computation of α itself, or through the computation of the sum ν x + ν y .
For the Potts point, the exact value is known to be α = 2/3. Due to logarithmic corrections in the correlation length, we have found 2ν = 1.426 which gives α = 0.57. If we make the hypothesis of a chiral transition, we can use the identityβ = ν x and for ∆ = 0.1 and ∆ = 0.3, we respectively found ν HT y +β = 1.42 and 1.427. We also directly studied α eff (see Fig. 16) and found evidence of a constant exponent along the transition at α 0.55, which establishes the validity of the hyper-scaling relation. Thus, we observe that α seems to be more or less constant along the chiral transition, as already observed in the p = 3 [37] case and in a quantum version of the model at λ 0.57 [34]. We note that a constant value of α = 2/3 along the transition and in particular at ∆ = 0 would imply that there are still logarithmic corrections in the correlation length along the chiral transition. Thus, we are probably overestimating the sum ν x + ν y , and in fact most certainly ν x since, in order to keep ν x + ν y = 4/3 with ν y > ν x , ν x should be smaller than 2/3.
The Kibble-Zurek exponent can be measured in arrays of Rydberg atoms [2]. It is defined as µ = ν/(1 + νz), with ν the correlation length exponent along the chain and z the dynamical exponent. In our case, z plays the role of the anisotropy exponent z = ν y /ν x , and ν is the correlation length exponent in the x-direction. Therefore, we can compute the Kibble-Zurek exponent as µ = ν x /(1 + ν y ). For reasons discussed above, we choose to express µ with respect to ν y . Then, using the hyper-scaling relation and the hypothesis of a constant α = 2/3 along the transition, the Kibble-Zurek exponent becomes µ = (4/3 − ν y )/(1 + ν y ). Using this expression of the exponent and the critical exponents obtained from the disordered phase for λ = 1, we found µ = 0.36 and µ = 0.34 for ∆ = 0.1 and ∆ = 0.3 respectively. These results indicate that the Kibble-Zurek exponent is smaller along the chiral transition than at the Potts point, where µ = 0.4. These results are in qualitative agreement with recent experiments [2] on chains of Rydberg atoms, which have been argued to be described by the AT model with λ close to 1 in the vicinity of the p = 4 ordered phase [34], and for which values µ 0.25 have been reported [2] .
Since our main purpose was to identify a possible chiral transition we have concentrated in the preceding section on the interval λ ∈ [0, 1] since it contains the relevant parameter range λ c2 0.32 ≤ λ < λ c1 = 0.9779. For completeness we now comment on the effect of a chiral perturbation outside this interval.
In this range, the chiral perturbation stays relevant and the floating phase open up as soon as the chirality is introduced. The zero temperature ground state stays four times degenerate and we do not expect any qualitative change from the λ < λ c2 case.

B.
λ < −1 In this case, at ∆ = 0, the nature of the zero temperature ground state has changed and so has the transition. It is expected to be described by Ising critical exponents. Indeed, in the λ → ∞ limit, the model is described by Ising variables u i = τ i σ i and the hamiltonian reduces to H = i∼j u i u j . The zero temperature ground state is infinitely degenerate and we recover a log(2) residual entropy. Introducing the chirality will lift the residual entropy but the zero temperature ground state stays 2 L degenerate for a L × M lattice size, and we expect the transition to stay a C-C one. We thus also expect the transition to stay in the Ising universality class. Other values of λ within the interval should be described by the same physics. We show results for λ = −2 at finite chirality ∆ = 0.1 in Fig. 17. And as expected we found a unique π commensurateπ commensurate transition with a correlation length exponent consistent with the Ising value ν = 1. For this simulation, we have changed the isometries and used those in [45], which we noticed to perform better in the anisotropic case.

C. λ > 1
At ∆ = 0, for λ > 1, the self dual line does not represent the critical line anymore and the Ashkin-Teller transition splits up into two Ising transitions with in between a nematic phase characterised by τ 1 σ 1 = 0. When introducing the chirality, the two Ising transitions first stay in the Ising universality class but as the chirality increases, the nematic phase closes and the two Ising transitions merge into a unique one bounded by commensurate phases, as illustrated in Fig.18 for λ = 1.7. Upon further increasing the chirality ∆, the high temperature phase becomes incommensurate, and for large enough ∆ a floating phase bounded by KT and PT transitions appears. For this value of λ, the range of chiral transition (if any) is too small to be detected. By contrast, for not too large λ, we found evidence for the existence of a chiral transition to follow up on the AT one. This result is illustrated in Fig. 19. In the high chirality limit, one recovers a two-step transition separated by a floating phase (not shown). Based on the results illustrated in Fig. 18  and 19, a sketch of the expected generic phase diagram as a function of ∆ is shown in Fig. 20. Whether the chiral transition persists for all λ is an open question that would require further investigation.

VII. SUMMARY
Based on the qualitative study of effective exponents, we have been able to give a self consistent picture of the critical exponents and universality classes, and to identify two different types of commensurate-incommensurate transitions, either a two-step one through a KT and a PT transition, or a unique chiral transition. If we take for granted that for λ c1 0.9779 ≤ λ ≤ 1, the transition should remain AT for a while, something we have not been able to show numerically because this range is too small, we are led to the conclusion that the Ashkin-Teller family of models can be classified into three regimes according to the way they react to a chiral perturbation: (i) A floating phase opens right away for 0 ≤ λ < λ c2 ; (ii) The transition is unique in the chiral universality class for λ c2 ≤ λ < λ c1 ; (iii) The transition remains in the AT universality class before becoming chiral for λ c1 ≤ λ ≤ 1. The critical value for λ c2 could not be pinned down very precisely, but it can be expected to satisfy λ c2 0.42. We further found evidence for the chiral transition to be characterised by a specific heat exponent that keeps its AT value. In many models, the melting of an ordered period-4 phase into incommensurate phase takes place and the nature of the transition is still debated. If an Ashkin-Teller point is present along the transition, by identifying its universality class and its exponent ν, our prediction can be used to determine the nature of the melting close to that point. Further study on such a model and in particular on the quantum chiral Ashkin-Teller model would be of particular interest to test these predictions.