Two-dimensional dilute Baxter-Wu model: Transition order and universality

We investigate the critical behavior of the two-dimensional spin-1 Baxter-Wu model in the presence of a crystal-ﬁeld coupling (cid:2) with the goal of determining the universality class of transitions along the second-order part of the transition line as one approaches the putative location of the multicritical point. We employ extensive Monte Carlo simulations using two different methodologies: (i) a study of the zeros of the energy probability distribution, closely related to the Fisher zeros of the partition function, and (ii) the well-established multicanonical approach employed to study the probability distribution of the crystal-ﬁeld energy. A detailed ﬁnite-size scaling analysis in the regime of second-order phase transitions in the ( (cid:2), T ) phase diagram supports previous claims that the transition belongs to the universality class of the four-state Potts model. For positive values of (cid:2) , we observe the presence of strong ﬁnite-size effects, indicative of crossover effects due to the proximity of the ﬁrst-order part of the transition line. Finally, we demonstrate how a combination of cluster and heat-bath updates allows one to equilibrate larger systems, and we demonstrate the potential of this approach for resolving the ambiguities observed in the regime of (cid:2) (cid:2) 0.


I. INTRODUCTION
Most of the commonly studied spin models of statistical mechanics such as the Ising and Potts or O(n) models are spin-inversion symmetric.A notable exception to this rule is the Baxter-Wu (BW) model [1,2] that was originally introduced by Wood and Griffiths [3,4].The commonly studied version is defined on the triangular lattice with N sites and the Hamiltonian function where J > 0 denotes a ferromagnetic exchange coupling.The sum ⟨ijk⟩ extends over all elementary triangles, and σ i = ±1 are Ising like spin-1/2 variables.The presence of three-spin interactions leads to the mentioned violation of spin-inversion symmetry, and it results in a fourfold degeneracy of the ground state: there is one ferromagnetic state with all spins up, and three ferrimagnetic states with down spins in two sublattices and up spins in the third.The triangular lattice can be decomposed into three sublattices, A, B, and C, as shown in Fig. 1 below.Note that the model of Eq. ( 1) is selfdual [3,4], resulting in the same critical temperature as that of the spin-1/2 Ising model on the square lattice, i.e., k B T c /J = 2/ ln ( √ 2 + 1) = 2.269185 • • • , where k B denotes Boltzmann's constant.
An exact solution of the model was provided early on by Baxter and Wu [1,2], supplying the critical exponents α = 2/3, ν = 2/3, and γ = 7/6.In the following, it was also shown that its critical behavior corresponds to a conformal field theory with central charge c = 1 [6,7].Due to the four-fold symmetry of the ground state, it is expected that the critical behavior of the q = 4 model Potts belongs to the universality class of the Baxter-Wu model [8].While, therefore, the critical exponents of the two models are identical, the same does not apply to the scaling corrections: the 4-state Potts model exhibits arXiv:2304.12379v2[cond-mat.stat-mech]7 Aug 2023 logarithmic corrections with the system size [9], whereas the Baxter-Wu model has power-law corrections with a correction-to-scaling exponent ω = 2 [6,7].Recently, the model has attracted renewed attention, and various aspects of its critical behavior have been studied in substantial detail [10][11][12][13][14][15][16][17][18].
A natural generalization of the Baxter-Wu model (1) results from the consideration of three spin orientations σ i = {−1, 0, 1} and the inclusion of an extra crystalfield (or single-ion anisotropy) coupling ∆.The resulting Hamiltonian then reads where E J and E ∆ denote the contributions of the exchange and the crystal field, respectively, to the total energy.Note that in the following we will use reduced units where J = 1 as well as k B = 1.This choice of units follows the practise of some of the present authors to implement multicanonical simulations on spin-1 Blume-Capel and Baxter-Wu models, see the discussion in Sec.III.Although still rather simple, for this spin-1 model there exists no exact solution, except for the case ∆ → −∞, where only configurations with σ i = ±1 are allowed and the pure spin-1/2 Baxter-Wu model is recovered, and also at zero temperature, where the four ordered phases coexist with the paramagnetic phase in a multiphase point at ∆/J = 2 (accordingly, no transition is observed for ∆/J > 2).
Based on the analogy between the Baxter-Wu and the diluted Potts model [19], but also on a series of more recent results [20][21][22], it is now well established that the phase diagram of the spin-1 Baxter-Wu model in the (∆, T )-plane includes a multicritical point separating first-from second-order transition regimes.This is in contrast to an earlier prediction by finite-size-scaling applied to transfer-matrix calculations, where a continuous transition only occurs in the limit ∆ → −∞ [23].In this respect, the model resembles the well-known Blume-Capel ferromagnet [24], which exhibits a phase diagram with ordered ferromagnetic and disordered paramagnetic phases separated by a transition line with first-and second-order segments (the latter in the Ising universality class) connected by a tricritical point, whose location is known with high accuracy [25][26][27][28].In contrast, there is no consensus on the precise location of the multicritical point for the spin-1 Baxter-Wu model -see Fig. 2 but also Fig. 5 of Ref. [22] for a summary regarding the phase diagram.Along the first-order transition line of Fig. 2 three ferrimagnetic phases and one ferromagnetic one coexist with the paramagnetic phase, forming a quintuple line that arrives at a pentacritical point where all five phases become identical.
In addition to the question of the location of the multicritical point, the reign of universality along the secondorder segment of the transition line has been put into question.While some earlier results based on the transfer matrix and conformal invariance suggested a contin- uous variation of critical exponents with the crystal field along the second-order transition line [21], more recent studies reported a match of the observed critical behavior with that of the 4-state Potts model [20].Some authors have also suggested the scenario of a mixed-order transition with both first-order and second-order properties [29].Recently, however, some clear-cut evidence for a simple, continuous transition in the universality class of the 4-state Potts model in the regime of ∆ < 0 could be provided based on a highly optimized combination of Wang-Landau simulations that cross the transition at constant ∆ and multicanonical simulations operating at constant temperature T [5,30], such that questions remain only in the regime ∆ ≳ 0.
In the present work we study the spin-1 Baxter-Wu model at several values of the crystal-field coupling that also reach into the regime ∆ > 0. To this end, we employ two complementary Monte Carlo schemes, a recently proposed variant of studying Fisher's partition function zeros [31] dubbed energy probability distribution zeros [32], and the multicanonical approach applied to the crystal-field energy [5,26,30,33].While the latter method is well established in the literature, the former was to date shown to be robust and useful in only a few cases, including some simple spin systems and polymer chains [32,34,35].In this respect, the purpose of the present work is twofold: firstly, to explore the scope and limitations of the method of energy probability distribution zeros for the more complicated spin-1 Baxter-Wu model that lacks the up-down symmetry and, secondly, to investigate the criticality and universality of the model specifically for ∆ ≥ 0, but still below the proposed location of the multicritical point.
The rest of the paper is organized as follows.In Sec.II we outline the method based on the energy probability distribution zeros and show results for both the pure spin-1/2 and for the spin-1 Baxter-Wu model, the latter for several values of the crystal-field coupling in the range −10 ≤ ∆ ≤ 0.5.In Sec.III we complement the outcomes of Sec.II via extensive multicanonical simulations at fixed values of the temperature in the regime where ∆ ≥ 0. Finally, in Sec.IV we summarize the main findings of the current work, comparing the implemented methodologies in the light of some additional preliminary results at ∆ = 0 obtained via an efficient numerical scheme that mixes cluster and heat-bath updates.

II. ENERGY PROBABILITY DISTRIBUTION ZEROS
A. Description and finite-size scaling As was recently discussed in Refs.[32,34,35], the study of zeroes in the energy probability distribution (EPD) allows for a straightforward determination of critical temperatures and the shift exponent θ = 1/ν while avoiding the need of computing traditional thermodynamic quantities, such as the susceptibility or the specific heat.The method of EPD zeros is closely related to the Fisher zeros of the canonical partition function Z [31], expressed as where E is the energy of the system, g(E) is the number of states having energy E (degeneracy), and β = 1/T .In the last part of Eq. (3) we assume a discrete energy spectrum E = E n = ϵ 0 + nϵ, where ϵ 0 is the ground-state energy and ϵ denotes the level spacing, n = 0, 1, 2, . . ., N .
Fisher noted that since (3) is a polynomial in y = e −βϵ , it has N complex zeros and since g(E n ) ≥ 0 none of them are real.However, on approaching the thermodynamic limit N → ∞, some zeros might approach the real axis, thus leading to a non-analyticity at y c = e −βcϵ corresponding to the phase transition at the inverse critical temperature β c .Since the partition function is not so straightforward to sample in a Monte Carlo simulation, we consider a somewhat different formulation.To this end, we multiply the right hand side of Eq. ( 3) by 1 = e −β0ϵ0 e +β0ϵ0 to obtain where β 0 is the inverse of some reference temperature, , and x = e −△βϵ .Note that h β0 (n) is the unnormalized canonical energy probability distribution at β 0 , and it can be easily estimated from an energy histogram through Monte Carlo simulations.As is easily seen, when β 0 = β c the dominant zero of the EPD is located at the fixed value x c = (1, 0) in the complex plane, thus simplifying the analysis.
For the finite lattice systems of linear size L that are amenable to numerical simulation, one can systematically follow the behavior of the most dominant zero x * L that is approaching the real axis at x c = (1, 0) as L is increased.In this way, it is possible to use finite-size scaling arguments to retrieve the critical temperature as well as the critical exponent ν.Specifically, the algorithm proposed in Ref. [32] for this purpose is as follows.We first choose a starting guess β j=0 0 of the inverse transition temperature and then iterate through the following steps: (1) Simulate the system at β = β j 0 and construct a histogram h β j 0 .
(2) Find all the zeros of the polynomial with coefficients given by h (3) Find the dominant zero (x j ) * .Then: x * L = (x j ) * and stop; (ii) else, make and return to step (1).
After setting a convergence criterion, one ends up with estimates x * L for the dominant zeros and hence with pseudocritical temperatures T * L , resp.β * L .Previous numerical results for several spin systems of Ising, Potts, and Heisenberg type, as well as for homopolymeric models, indicate that [32,34,35]: (i) the choice of the starting temperature β 0 is largely irrelevant for arriving at the dominant zero x * L ; (ii) for β ≈ β c , only states with a non-vanishing probability are relevant to the transition, thus allowing to define a cutoff, h cut , affecting the leftand right-hand side margins of the energy distribution.Discarding configurations where h β j 0 (n) < h cut substantially reduces the degree of the polynomial, especially with increasing system size; and (iii) to further simplify the polynomial, one can rescale the histogram by setting its maximum value to unity so that max h β j 0 = 1, since an overall rescaling of the partition function does not affect the location of the zeros.
According to the well-established finite-size scaling theory, the shift of pseudocritical temperatures T * L is described by the power law [36,37] where T c is the critical temperature of the infinite system, b and b ′ are non-universal parameters, ν is the critical exponent of the correlation length, and ω denotes the correction-to-scaling (Wegner) exponent, fixed hereafter to the predicted value ω = 2 [6,7].On the same ground, one also expects that [32] x Since x * L ≈ (1, 0), the imaginary part Im(x * L ) should scale with the system size as [32] In this description, the standard process is to firstly compute the critical exponent ν via Eq.( 9), and then retrieve the critical temperature T c using the Eq. ( 7).

B. Results
For the application of the EPD zeros method to the Baxter-Wu model, histograms were accumulated using the standard single-spin-flip Metropolis algorithm.To accommodate for all ground states, periodic boundary conditions must be considered and the allowed values   Finite-size scaling analysis of the pseudocritical temperatures, see Eq. (7), where the critical exponent ν is fixed to the value from panel (a).
of the linear size of the lattice L must be a multiple of three [5,20,21,30].(Note that a triangular lattice on the torus is tripartite when its linear dimensions are multiples of three).In the course of our simulations we considered linear sizes in the range 12 ≤ L ≤ 120, respecting this constraint (a practise followed also in the multicanonical and hybrid simulations described below).During thermalization, 10 5 Monte Carlo steps per spin (sweeps) were discarded for L ≤ 45 and 3×10 5 sweeps for the larger sizes.An additional 10 8 sweeps were performed to accumulate the energy histograms, leading to a quite precise estimate of the dominant root.The iterative process of finding the dominant EPD zero terminated when the temperature difference between two consecutive steps became smaller than a predefined accuracy of ε = 10 −4 .Note that one may also look at the real part of the dominant zero and halt the process when | Re(x * L )−1| ≤ ε, and also that smaller values of ε may be considered, without any significant consequences in the results.Regarding the cutoff, the value h cut = 10 −4 was used throughout the simulations.Errors have been computed by averaging over 10 different independent runs.Finally, for all fits performed throughout this paper we restricted ourselves to data with L ≥ L min , adopting the standard χ 2 test for goodness of the fit.Specifically, we considered a fit as being acceptable only if 10% < Q < 90%, where Q is the quality-of-fit parameter [38].
It is clear that finding the zeros of a high-order polynomial is far from trivial, in general.In the present case, this difficulty is being added to by the need to determine the cutoff of the EPD while monitoring the precision of coefficients h β j 0 (n) in order to obtain a sensible accuracy for the zeros.The precision of the coefficients h β j 0 (n) strongly depends on the length of the Monte Carlo time series.On the other hand, even in case of rather accurate values of the coefficients results still depend on the cutoff threshold of the EPD.However, as it has been recently shown for the two-dimensional Ising and 6-state Potts models [39], the EPD method is indeed quite robust against the cutoff threshold and the number of Monte Carlo sweeps used, giving accurate results for the critical parameters.
Since the method has not yet been checked on the spin-1/2 Baxter-Wu model, our first port of call is to test it against the well-known exact results.Figure 3 depicts typical results for the zeros of a system L = 75 at a temperature close to its pseudocritical.Note that since the density-of-states factors g(E n ) are real, the zeros all come in conjugate pairs.The upper panel shows a global view of all zeros located around a unit circle, and the bottom panel depicts an enlargement of the dominant zero (x j ) * near the area x c = (1, 0).Changing the temperature ac- cording to Eq. ( 6) will furnish a different new dominant zero that converges to the desired x * L after just a few iterations.The finite-size scaling analysis of the imaginary part of the dominant root, as given by Eq. ( 9), is shown in Fig. 4(a).A linear fit in a log-log scale gives an estimate of ν = 0.668 (6) for the critical exponent of the correlation length, in very good agreement with the exact result ν = 2/3 [1,2].Fixing ν to this value and ω = 2 [40], Eq. ( 7) gives T c = 2.2692(4) in accordance with the exact result 2.269185 • • • [2], see Fig. 4(b).
We proceed now with the study of the spin-1 Baxter-Wu model at ∆ = {−10, −1, 0, 0.5}.This selection allows a direct comparison with results already reported in the literature by other approaches [5,20,29,30].For brevity, we chose to show here in Fig. 5 the case ∆ = 0.5, which is the largest positive value of ∆ considered in this work.The scaling analysis in both panels of Fig. 5 is in direct analogy with that of Fig. 4, giving ν = 0.623 (11) and T c = 1.5301 (3).Although the estimate for T c is in excellent agreement with conformal invariance, see Tab.I, the value of ν appears to deviate from the expected 2/3 result.A similar but slighter deviation was observed also for the case ∆ = 0, see again Tab.I.This trend is a note of warning indicating the presence of strong finite-size effects as ∆ approaches the location ∆ pp of the multicritical point, suggesting the need of studying larger system sizes.Finally, in Fig. 6 we provide a summary concerning the finite-size scaling behavior of the imaginary part of the dominant zero for all values of ∆ considered, including the case of the spin-1/2 model.Inspecting Fig. 6 one may observe that as we lower ∆ from 0.5 to −10 the trend of the numerical data follows the expected passage to the spin-1/2 model (∆ = −∞).However, this approach appears to be rather slow, and it could instructive to study even more negative values of ∆.The gathered results for T c and ν are listed in Tab.I and are critically discussed in Sec.IV.Overall, we may deduce that the EPD zeros method appears to be a promising alternative for deter- mining critical aspects of the transition in the Baxter-Wu model.

III. MULTICANONICAL SIMULATIONS A. Method and observables
The multicanonical (MUCA) method [33] consists of a substitution of the Boltzmann factor e −βE with weights that are iteratively modified to produce a flat histogram, usually in energy space.This ensures that suppressed states such as those in the co-existence region in an (effectively) first-order transition scenario can be reliably sampled, and a continuous reweighting to arbitrary values of the external control parameter becomes possible [41,42].Due to the two-parametric nature of the density of states, g(E J , E ∆ ), in the spin-1 Baxter-Wu model, the process was applied only to the crystal-field part E ∆ of the energy.This allowed us to reweight to arbitrary values of ∆ while keeping the temperature fixed.Starting from the partition function of Eq. ( 3) we can write where the Boltzmann weight associated with the crystalfield part of the energy has been generalized to W (E ∆ ).
For a flat marginal distribution in E ∆ , it should hold that In order to iteratively approximate the generalized weights W (E ∆ ), we sampled histograms of the crystalfield energy.Supposing that at the n th iteration a histogram H (n) (E ∆ ) was sampled, then its average should depend on the weight of the iteration W (n) (E ∆ ) as From Eqs. (11) and ( 12) it follows that Hence, in order to approximate the W (E ∆ ) that produces a flat histogram a weight modification scheme of the form The simulations can terminate when a flat-enough histogram has been sampled, based on a suitable flatness criterion.For our purposes we used the Kullback-Leibler divergence to test the flatness [42,43].After this initial preparatory part, the final fixed weights can be used for production runs.
As has been shown in detail in Refs.[42,44], the multicanonical method can be adapted for the use on parallel machines by performing the sampling of histograms in parallel, with each parallel worker using the same weights but a different (independent) pseudorandom number sequence.The accumulated histogram can then be used to update the weights, keeping communication between the parallel parts of the code minimal.This scheme has been successfully applied for the study of spin systems in the past, including the spin-1 Blume-Capel and Baxter-Wu models [5,26,28,30,45].Here we performed our simulations on an Nvidia Tesla K80 GPU, using a total of 26 624 workers assigned to independent copies of the system.At each time, a subset of these threads are actually running in parallel on the 4 992 cores of the device, while the excess in the number of parallel tasks is employed to hide the latencies due to memory accesses [46].In the course of the multicanonical simulations the sampled observables include estimates of the mean energy ⟨E⟩, the order parameter ⟨m⟩ which is estimated from the root mean square average of the magnetization per site of the three sublattices A, B, and C [29,47,48], and the magnetic susceptibility As the multicanonical method allows for continuously reweighting to any value of ∆, canonical expectation values for an observable O = O({σ}) at a fixed temperature can be attained by estimating the expressions In this framework, it is natural to compute ∆-derivatives of observables rather than the usual T -ones.For instance, in place of the usual specific heat one may define a specific-heat-like quantity [26] which shows the shift behavior expected from the usual specific heat [5,26,30].Additionally, in order to obtain direct estimates of the critical exponent ν from finite-size scaling, one may compute the logarithmic derivatives of nth power of the order parameter [36,49,50]

B. Results
We performed simulations at T = 1.6606 and T = 1.5301 which approximate the critical points at ∆ = 0 and ∆ = 0.5 respectively, see Table I, using system sizes 12 ≤ L ≤ 96, again with periodic boundary conditions.For T = 1.6606, 4 × 10 6 sweeps were used in the production run for the smallest system and 3 × 10 8 sweeps for the largest.For the lower temperature T = 1.5301, due to its proximity to the proposed multicritical point (see Fig. 2), sampling was increased to 2.5 × 10 7 sweeps in the production run for the smallest system and 10 9 sweeps for the largest.After the initial iterations for the calculation of the generalized weights, an additional 10% of the total production sweeps were discarded by each worker for thermalization.Preliminary tests indicated that distributing the production sweeps equally among the workers results in sampling the equivalent of ∼ 5 autocorrelation times worth of data points per worker for the larger temperature and ∼ 20 for the smaller.The results were analyzed using the jackknife resampling method [51] and the location of pseudocritical points was estimated via reweighting and bisecting in ∆.
As discussed in Refs.[5,22,29], there have been recent reports of first-order transition features even along the presumed continuous part of the transition line.In relation to such claims, we put forward here some additional evidence for the clarification of the nature of the phase transition at ∆ < ∆ pp .Following the prescription of Ref. [5] we studied the reweighted probability density function P (E ∆ ).It is well known that a double-peak structure in the density function in finite systems is an expected precursor of the two δ-peak behavior in the thermodynamic limit that is expected for a first-order phase TABLE I: Representative critical-point estimates (Tc or ∆c) of the phase diagram of the spin-1 Baxter-Wu model from the present work and as well as previous studies, including estimates of the critical exponent ν.In the first column we either indicate the value of ∆ for simulations that vary T or the value of T for simulations that vary ∆.Columns 2 -5 feature results obtained in the current work from the EPD zeros method (columns 2 and 3) and multicanonical simulations (columns 4 and 5).Columns 6 -9 append earlier estimates from Wang-Landau (WL) simulations (∆ = −10 and −1 [5] and ∆ = 0 [29]) and conformal invariance (CI) [20].The first line of results relates to the pure spin-1/2 model (∆ = −∞).b This estimate (and the one at T = 1.5301) corresponds to the average value of ν obtained from the fits of Fig. 9. Cross-correlations were not taken into account, but see Ref. [60].

EPD
c Private communication by the authors of Ref. [20].
transition [52,53].However, this observation must be taken with a grain of salt, since there have been many cases reported in the literature, for which this two-peak structure tends to a unique peak in the thermodynamic limit.A warning example is the two-dimensional 4-state Potts model [54].
We start the presentation of our results with Fig. 7(a) where we show the probability density function P (E ∆ ) for selected system sizes at the temperatures T = 1.6606 and 1.5301.A double-peak structure is observed in both cases, in agreement with the evidence in Ref. [29] for ∆ = 0.As is clearly visible, stronger first-order-like characteristics are present for the lower-T (higher-∆) example that is closer to the multicritical point.
The multicanonical method is optimal for studying these phenomena in the framework of the method proposed by Lee and Kosterlitz [55], as it allows the direct estimation of the barrier associated with the suppression of states during a first-order phase transition.Considering distributions with two peaks of equal height (eqh) [56], as the ones shown in Fig. 7(a), allows one to extract the surface tension in the E ∆ -space, where P max and P min are the maximum and local minimum of the distribution P (E ∆ ), respectively.This parameter is expected to scale in two dimensions as possibly with higher-order corrections [57][58][59].For the system under investigation here, the scaling behavior of the surface tension Σ(L) is depicted in Fig. 7(b) for all temperatures studied.The dashed lines show fits of the form (19) with L ≥ L min = 30 leading to a practically zero value of Σ ∞ in all cases.In particular, we obtain the extrapolated values Σ ∞ = −0.00005(11), −0.0003( 9), and 0.0003(3), for T = 1.8503, 1.6606, and 1.5301, respectively.This analysis suggests a continuous transition in the thermodynamic limit for the regime of ∆ < ∆ pp , in favor of the scenario originally discussed in Ref. [5].
In order to extract critical crystal fields ∆ c (T ) as well as a first estimate of the correlation-length exponent ν, we present in Fig. 8 the shift behavior of suitable pseudocritical fields ∆ * L .These are defined as the peak locations of ∆-dependent curves, such as the specific heat C ∆ , the magnetic susceptibility χ, and the logarithmic derivative of the order parameter ∂ ln ⟨m⟩/∂∆.For each of the two temperatures studied the dashed lines show joint fits to the expected power-law behavior [26,28] where ∆ c and ν are common parameters and ω = 2 [6,7].Using L min = 24 and 48 for T = 1.6606 and 1.5301, respectively, the evaluated critical points ∆ c (T = 1.5301) = 0.4999(2) and ∆ c (T = 1.6606) = 0.0008 (7) are in good agreement with the results of Sec.II B but also with those reported in Tab.I from Wang-Landau simulations [29] and conformal invariance [20].More importantly, our estimates ν = 0.656 (28) and 0.632 (45) for T = 1.5301 and 1.6606, respectively, confirm to a good accuracy the q = 4 Potts model universality class [9].Additional estimates for the critical exponent ν can be extracted from the maxima of the logarithmic derivatives of the order parameter according to Eq. ( 17), which are expected to scale as [36,49,50] Figure 9 shows our data for n = 1 (main panel) and n = 2 (inset) at the two temperatures under study.The dashed lines are power-law fits of the form (21) using L min = 36, providing an average of ν = 0.654 (19) and 0.652 (10) for T = 1.5301 and 1.6606, respectively, thus reinforcing the scenario of the q = 4 Potts model universality class [8].Finally, we turn to the finite-size scaling behavior of the maxima of the specific heat, C * ∆ , and magnetic susceptibility, χ * , in order to probe the critical-exponent ratios α/ν and γ/ν, respectively.Figure 10 contains the relevant numerical data at the two temperatures considered.The dashed lines are fits of the expected form [5,28] and with L min = 36.These led to the estimates α/ν = 1.13(5) and 1.06(3), and γ/ν = 1.74(4) and 1.75(3) for T = 1.5301 and 1.6606, respectively.Here, at the lower temperature T = 1.5301 we had to include a secondorder correction term (∼ L −2ω ) in our fitting attempts to improve the quality-of-fit.All of the above results are clearly compatible with the exact values α/ν = 1 and γ/ν = 7/4 of the 4-state Potts model universality class [9].

IV. SUMMARY AND OUTLOOK
In closing, we return to the question of the current understanding of the behavior of the model along the II: Summary of results for the critical exponent ν of the spin-1 Baxter-Wu model at ∆ = 0 obtained via conformal invariance (2 nd row) [20], Wang-Landau simulations (3 rd row) [29], EPD zeros (4 th row), and multicanonical simulations (5 th row).The 6 th row showcases the estimate of ν via the hybrid approach [61,65], see also Fig. 11.For all methods the maximum accessible system size Lmax used in the simulations is also given in brackets.The 3 rd column highlights the deviation δν of each estimate from the exact value [9], which is included in the last row for reference.Exact solution 2/3 0 a Note that Lmax denotes the maximum strip width considered in Ref. [20].

Method
b Average value of ν obtained from the fits of Fig. 9.
phase boundary.Our results as well as some reference estimates from the recent literature are summarized in Table I.On inspecting these values, the following comments are in order: (i) A very good agreement between different methods of estimating the location of points (∆, T ) along the phase boundary of the model is observed, cross-validating the different numerical approaches used in the present but also in previous works.(ii) The values of the critical exponent ν at ∆ < 0 are fully compatible with the value 2/3 of the 4-state Potts universality class [8,9] for all methods.However, with increasing ∆, a slight decrease in the value of ν is observed and may be attributed to the presence of finite-size effects that become more pronounced as one approaches the pentacritical point ∆ pp ≈ 0.89 − 1.68 [20,22].(iii) Although the multicanonical simulations allowed us to significantly improve the limited capability of the Metropolis algorithm to reduce correlations, much larger system sizes are required for a safe determination of critical exponents, in particular in the regime 0 ≤ ∆ < ∆ pp .
In light of the above discussion, it would be very valuable to have at one's disposal some simulation method that allows to equilibrate significantly larger systems than those considered here.This holds especially for the scaling at the pentacritical point itself, where one may need to take into account possible multiplicative and additive logarithmic corrections, similar to those present in the 4-state Potts model.A suitable cluster update for the spin-1/2 Baxter-Wu model was proposed by Novotny and Evertz [61].Its basic idea is as follows: for each update step one of the sublattices is chosen at random and its spins kept fixed, resulting in an effective Ising model on the other two sublattices with non-frustrating couplings.Hence, the Swendsen-Wang [62] algorithm can be applied for simulations of these embedded models.Our preliminary tests indicate that a combination of this cluster formalism that improves the decorrelation of configurations but is not ergodic as it does not affects the diluted spins σ i = 0 with the heat bath algorithm [63,64] results in an efficient algorihtm capable of thermalizing rather large systems.A detailed analysis of the critical dynamical behavior of this hybrid scheme will be presented elsewhere [65].
Here, we confine ourselves to an exemplary application of this technique to the case ∆ = 0 beyond which the deviation in the estimates of ν from the expected value 2/3 appears to grow, cf. also Table I.In Fig. 11 we present the results of a test calculation of the critical exponent ν from hybrid simulations at the critical point (∆, T ) = (0, 1.6606), studying systems up to linear size L max = 240.The finite-size scaling analysis of the logarithmic derivatives of the order parameter (for both n = 1 and 2) produces the estimate ν = 0.669(3), in excellent agreement with the value 2/3 [1,8,9].A comparative set of results for the critical exponent ν of the spin-1 Baxter-Wu model at ∆ = 0 is given in Tab.II, where one may notice the superior accuracy of the hybrid approach.
To conclude, this work complements previous results that map the universality class of the spin-1 Baxter-Wu model to that of the 4-state Potts model, a nontrivial task obscured by the presence of strong finite-size effects as revealed by our analysis.Clearly, it would be very instructive to add data for additional values of ∆ in the regime ∆ > 0.5.Yet this would require a huge computational effort given that crossover phenomena become more pronounced as we move towards the expected pentacritical point.This indicates that in order to perform a safe finite-size scaling analysis much larger system sizes would be needed with increasing values of ∆.For future work, we propose the following two-stage process: (i) identify with good numerical accuracy the location of the pentacritical point (∆ pp , T pp ) and (ii) perform extensive simulations around this point using the hybrid approach in order to quantify all these interesting phenomena outlined above, including crossover effects and possible logarithmic corrections to scaling.A possible tool for such an endeavor could be the field-mixing technique [66] in combination with the numerical methods reported in this paper.Such attempts are the subject of ongoing investigations.

FIG. 3 :
FIG.3: Zeros of the EPD in the complex plane for the spin-1/2 Baxter-Wu model and a system with linear size L = 75 at β = 0.43991, close to the corresponding pseudocritical temperature; see also Eq. (5).Panel (a) illustrates a global view of the zeros around a unit-radius circle and panel (b) a zoom in on the area of the dominant zero near (1, 0).As discussed in the main text the roots appear in conjugate pairs.

FIG. 4 :
FIG. 4: (a) Log-log plot of the imaginary part of the dominant zero as a function of the lattice size for the spin-1/2 Baxter-Wu model.The solid line shows a fit of the form (9). (b)Finite-size scaling analysis of the pseudocritical temperatures, see Eq.(7), where the critical exponent ν is fixed to the value from panel (a).

5 FIG. 6 :
FIG.6: Log-log plot of the imaginary part of the dominant zero as a function of the lattice size for both the spin-1/2 and spin-1 Baxter-Wu models.For the spin-1 case results at various values of ∆ are shown.The solid lines are fits of the form (9) as described in the text.

FIG. 8 :
FIG.8: Shift behavior of several pseudocritical fields as a function of the inverse linear system size at the temperatures T = 1.5301 and T = 1.6606, corresponding roughly to ∆ ≈ 0.5 and ∆ ≈ 0, respectively.

FIG. 11 :
FIG.11: Finite-size scaling of the logarithmic derivatives of powers n = 1 and 2 of the order parameter at the critical point (∆, T ) = (0, 1.6606), as indicated by the EPD zeros method; see also Tab.I.The dashed lines are fits of the form(21).Data generated via the hybrid approach.