Multicanonical analysis of the plaquette-only gonihedric Ising model and its dual

The three-dimensional purely plaquette gonihedric Ising model and its dual are investigated to resolve inconsistencies in the literature for the values of the inverse transition temperature of the very strong temperature-driven first-order phase transition that is apparent in the system. Multicanonical simulations of this model allow us to measure system configurations that are suppressed by more than 60 orders of magnitude compared to probable states. With the resulting high-precision data, we find excellent agreement with our recently proposed nonstandard finite-size scaling laws for models with a macroscopic degeneracy of the low-temperature phase by challenging the prefactors numerically. We find an overall consistent inverse transition temperature of 0.551334(8) from the simulations of the original model both with periodic and fixed boundary conditions, and the dual model with periodic boundary conditions. For the original model with periodic boundary conditions, we obtain the first reliable estimate of the interface tension, 0.12037(18), using the statistics of suppressed configurations.


Introduction
The gonihedric Ising model was originally formulated by Ambartzumian and Savvidy [1] as a possible lattice discretization of an alternative "linear" action for the string worldsheet in bosonic string theory. These early discretizations using triangulations were then translated to plaquette surfaces generated as the spin cluster boundaries of classical generalized Ising spin Hamiltonians by Savvidy and Wegner [2]. For a recent review, see [3].
The resulting gonihedric Ising model is a generalized three-dimensional Ising model, where spins σ, living on a three-dimensional cubic lattice, interact via nearest i, j , next-to-nearest i, j and plaquette interactions [i, j, k, l] and the weights of the different interactions are finetuned so that the area of spin cluster boundaries does not contribute to the partition function, but rather their edges and self-interactions. This leads to the family of Hamiltonians where κ parametrizes the self-avoidance of the spin cluster boundaries. The purely plaquette Hamiltonian with κ = 0, allows spin cluster boundaries to intersect without energetic penalty. It has attracted some attention recently, as it displays a strong first-order transition [4] and evidence of glass-like behaviour [5] at low temperatures in spite of the absence of quenched disorder. The strong first-order nature of the transition for the purely plaquette Hamiltonian has meant that it has proved difficult to obtain consistent values for the inverse transition temperature. Only recently, it was understood that the exponential degeneracy q = 2 3L [6] in the low-temperature phase, for the periodic system living on a cube with linear lattice size L, severely changes the finite-size corrections and leads to nonstandard finite-size scaling [7].
First estimates of the inverse transition temperature were given by a mean-field approximation that yielded β ∞ = 0.325 and early canonical Monte Carlo simulations gave β ∞ = 0.50(1) [8]. Later, detailed finite-size scaling analysis of Monte Carlo data with fixed boundary conditions found β ∞ = 0.54757(63) [9]. More recently, another value of 0.510 (2), that is apparently compatible with the first simulations [8], was suggested by analysing a dual representation of the model with periodic boundary conditions which turns out to be an anisotropic Ashkin-Teller model [10]. Here, two spins σ, τ live on each vertex of a cubic lattice, with nearest-neighbour interactions along the x, y and z-axes, To resolve these discrepancies, we have conducted multicanonical Monte Carlo simulations of the original model and its dual representation. In the remainder of the paper we present the results of these simulations and the underlying theory of finite-size scaling at first-order phase transitions used in extracting the conclusions.
In Section 2, after first discussing "standard" first-order finite-size scaling behaviour for models with periodic boundary conditions, we observe that nonstandard first-order finite-size scaling behaviour will occur when there is an exponentially large degeneracy of the low-temperature phase, which is the case for both the plaquette Hamiltonian and its dual. The scaling corrections for various observables are presented in detail in this case.
In Section 3 we discuss the simulations themselves, for the plaquette Hamiltonian with periodic boundary conditions, for the dual model with periodic boundary conditions and for the plaquette Hamiltonian with fixed boundary conditions. We determine characteristic quantities of these systems, such as the specific heat, Binder's energy parameter, the latent heat and interface tension, as well as the inverse temperature of the transition. We find for the first time, with the use of the nonstandard scaling relations, a consistent estimate of the inverse critical temperature from the plaquette Hamiltonian and its dual with periodic boundary conditions and the plaquette Hamiltonian with fixed boundary conditions. The self-consistency of the simulations is also confirmed by extracting various prefactors of scaling corrections both directly from fits and by calculation of energy moments.
Finally, in Section 4 we summarize our results for the measured physical quantities in the gonihedric model, make some observations on the apparent dependence of the latent heat and interface tension on the boundary conditions and conclude the paper with a brief discussion of potential further applications of nonstandard first-order finite-size scaling.

Finite-size scaling for first-order phase transitions
In spite of the ubiquity of first-order phase transitions [11] it was only relatively recently that the initial studies of finite-size scaling for first-order transitions were carried out [12] and subsequently pursued further in [13]. Rigorous results for periodic boundary conditions were derived in [14,15]. For the discussion of scaling laws under periodic boundary conditions here, we will first employ the two-state model of [16] which is capable of correctly reproducing the prefactors of the leading contributions. To have this paper reasonably self-contained, we will recall the principles and main results in the following. In this model we assume the system spends some time in either one of the q ordered phases or in the disordered phase, where transitions between the phases are instantaneous and fluctuations within the phases are neglected. Let W o denote the fraction of the total time spent in one of the q ordered phases and W d = 1 − W o the fraction spent in the disordered phase. We associate energies e o and e d with the phases. Under these assumptions, the energy moments are simply given by e n = W o e n o + (1 − W o )e n d . The specific heat C V (β, L) = −β 2 ∂e(β, L)/∂β as an expression in terms of these moments becomes with e = e d − e o . Varying W o , we find the maximum for W o = W d = 0.5, where the ordered and disordered peaks of the energy probability density have equal weight. Here, we have assumed that β, e o and e d deviate from the infinite-volume limit β ∞ , ê o = e o (β ∞ ) and ê d , respectively, only by terms of order 1/L d .
In close analogy, we can write the Binder parameter in terms of the two-state moments from which we can calculate the minimum with respect to the weights at The free-energy densities, f o and f d , of any one of the ordered phases or the disordered phase govern their probability and the fraction of time spent in the ordered phases is proportional to qp o . Neglecting exponentially small corrections in the linear lattice size L [14][15][16], the ratio of the fraction of time spent in the respective phases is given by Expanding the logarithm ln( which after truncation can be solved for the inverse temperatures β. For W o = W d = 0.5 we find the inverse temperature β eqw , where both peaks of the energy probability density have equal weight, and, to leading order, the location β C max V of the specific heat maximum, The minimum of Binder's parameter at W o /W d =ê 2 d /ê 2 o is located at the inverse temperature In spite of its simplicity the model captures the essential features of first-order phase transitions and, importantly for our purposes, correctly predicts the prefactors of the leading finite-size scaling corrections for a class of models with a contour representation, such as the q-state Potts model, where a completely rigorous theory of scaling also exists [15]. This rigorous theory enables the calculation of the coefficients of higher-order terms in a systematic asymptotic expansion in powers of 1/L d [16,17]. In addition, there are further corrections that decay exponentially fast with growing system size [18].
These models for periodic boundary conditions have, up to exponentially small corrections in L, a canonical partition function of the form [15] allowing a more rigorous derivation of inverse transition temperatures. The inverse temperature of equal peak weight then reads [16] where ŝ = β ∞ ê is the transition entropy and Ĉ =Ĉ d −Ĉ o . For the location of the specificheat maximum and the minimum of the Binder parameter one finds [16,17] where the expression for a 2 is explicitly known but very complicated [16], and will simplify in our special case (see below). The leading correction to the inverse temperature of equal peak height, β eqh , comes more heuristically from a double gaussian approximation of the energy probability density [16], Usually, the low-temperature degeneracy q is a constant (like in a q-state Potts model) and standard finite-size scaling behaviour at a first-order transition displays a leading contribution proportional to the inverse volume L −d . However, for the three-dimensional gonihedric Ising model (1), the degeneracy q is exponentially dependent on the linear system size. By construction, the model shows a highly degenerate ground-state for all parameters κ. In the special case of vanishing energetic penalty for self-intersecting spin cluster boundaries, κ = 0, the degeneracy, is apparent even for finite temperatures [6]. The usual 1/L 3 behaviour is therefore transmuted to a 1/L 2 behaviour [7]. Furthermore, the finite-size scaling corrections to β eqw (L) in (14) and in the scaling law (15) for β C max V (L) now coincide up to order O(L −4 ), The scaling law for the peak location of Binder's parameter (16) becomes where we have used the fact that only the contribution to a 2 with the highest power of ln q, a 2 = (ln q/ ŝ) 2 Ĉ /2 ê + . . . , contributes to the order given. The leading contribution to the finite-size correction is thus also proportional to L −2 , and the prefactor of the contribution O(L −4 ) becomes the same as that found for the inverse temperatures of the equal peak weight and the peak location of the specific heat. Note that there is, however, an additional correction term of O(L −3 ). The leading correction to the inverse temperature of equal peak height, β eqh , is now also of order O(L −2 ), The extremal values of the specific heat and Binder's parameter change with the system size according to and The prefactor in the first correction for B min (L) reads which comes from an even more complicated expression of O(L −3 ) in the general case [16].
Here we have already taken the degeneracy q = 2 3L into account. The Taylor series of the energy in the ordered and disordered phases, e o/d , around β = β ∞ reads where the specific heat of the ordered and disordered phase enters the leading correction. Calculating the energies at inverse temperature β eqw , the scaling of the energy fulfils The difference β eqw − β ∞ is known from (19), therefore the expression shows the finite-size corrections to the energy. The same change of leading contributions is apparent in the dual model (3), where a similar low-temperature phase degeneracy is expected (but, in contrast to the original model, not proven). These considerations will also apply to other models with periodic boundary conditions which have a degeneracy that depends exponentially on the system size.
For fixed boundary conditions, surface effects play an important role [19]. Here, the inverse transition temperature of the gonihedric Ising model is shifted by a leading term of order O(L −1 ) for finite lattices of linear lattice size L. Thus in this case we expect for the peak locations of both the specific heat and Binder's parameter.

Simulation results
An effective way of combating supercritical slowing down near first-order phase transitions, where canonical simulations tend to get trapped in either the disordered or ordered phases, is to use the multicanonical Monte Carlo algorithm [20,21]. At such first-order transitions cooling a system down or heating it up makes estimation of the transition temperature inherently difficult using standard algorithms due to hysteresis effects. In a multicanonical simulation it is possible to systematically improve guesses of the energy probability distribution before the actual production run by using recursive estimates [22].
The rare states lying between the disordered and ordered phases in the energy histogram are then promoted artificially, decreasing the autocorrelation time and allowing the system to oscillate more rapidly between phases. Canonical estimators can then be retrieved by weighting the multicanonical data to yield Boltzmann-distributed energies. Such reweighting techniques [23] are very powerful when combined with the multicanonical simulations, allowing the calculation of observables over a broad range of temperatures.

Observables
Standard observables such as the specific heat (4) and Binder's energy parameter (6) were calculated at different temperatures from our data for both the gonihedric Ising model (2) and its dual (3). The positions of their peaks, β C max V (L) and β B min (L) were then determined by the reweighting techniques.
Other estimates of the inverse critical temperature are given by β eqw (L) and β eqh (L), where the ordered and disordered peaks of the energy probability density p(e) have the same weight or height, respectively. In practice, we use reweighting techniques to get an estimator of the energy probability densities at different temperatures. β eqw is chosen systematically to minimize (29) where the energy of the minimum between the two peaks, e min , is determined beforehand to distinguish between the different phases. The location of the minimum, β eqw , is then used to calculate the energy moments of the ordered and disordered phases, where e o/d (L) = e 1 o/d (L) is the energy in the respective phases, and their difference is an estimator of the latent heat e(L) = e d (L) − e o (L). Also, the second and first moments combine to give the specific heat of the ordered and disordered phases, To find the inverse transition temperature where both phases have equal height we minimize as a function of β. The probability density p(e, β eqh ) itself at β eqh is also of interest since one can make use of it to extract the reduced interface tension for periodic boundary conditions. This characteristic quantity of first-order phase transitions is almost impossible to extract reliably from canonical Metropolis simulations, since getting reasonable statistics on the suppressed states is very hard. Multicanonical simulations, on the other hand, are perfectly tailored for measurements of such rare events.

Original plaquette model with periodic boundary conditions
The quality of the simulations for the original plaquette model (2) with periodic boundary conditions can be judged by the integrated autocorrelation time τ int and the number of sweeps in Fig. 1. Here, τ int has been determined with a self-consistent cutoff at 6τ int and the error comes from the known formula for this algorithm [24], τ int = (τ int ) 3/2 ( 24 n ) 1/2 , where n is the number of measurements (= number of sweeps when performing measurements every sweep). We would expect that the integrated autocorrelation time with perfect multicanonical weights should in principle increase linearly with the volume, τ int ∝ L 3 .
Error-weighted nonlinear least-squares fits of a power law, τ int ∝ L α , to the actual measured integrated autocorrelation times yield much larger exponents α ≈ 6.40(17) that vary a bit around 6.0 depending on the lattice sizes included in the fits. Also for those fits with acceptable χ 2 ≈ 1 that only include lattices with size L ≥ 16, fits to an exponential growth with L are of comparable quality, see Fig. 1. With least-squares fits and no proper model testing, we cannot distinguish between the two alternatives [25]. In any case, we find that the autocorrelation time grows significantly faster than expected, an effect that may be attributed to free-energy barriers. Such hidden barriers appear for instance in droplet condensation [26], whose analog with the gonihedric Ising model is, however, still unclear. For each lattice, we performed a maximum number of n max = 2 17 L 3 = 131 072 L 3 sweeps with an upper bound on the computer time of around 500 hours of real time for each lattice size. All lattices with size L ≤ 20 finished the desired number of sweeps, the larger lattices were aborted after 500 hours and collected correspondingly less statistics. The ratio n/τ int is a measure for the number of effectively uncorrelated data. Although the autocorrelation time increases dramatically with the system size, the simulation of the largest lattice of V = 27 3 spins still effectively flipped more than 250 times between the two phases during the simulation. This is remarkable, given that rare states were suppressed by more than 60 orders of magnitude compared to the most probable states (see the inset of Fig. 2).
From our multicanonical data, we have extracted the resulting quantities of interest for different lattice sizes and listed them in Table 1, where errors have been extracted by jackknife analysis using 20 blocks for each lattice size [27].
We find that the inverse temperatures of the specific-heat maximum β C max V and of equal peak weights β eqw fall nearly together for all lattice sizes. This is accounted for by the fact that the higher-order corrections of order O(L −4 ) in the scaling law (19) for these quantities collapse due to the exponential degeneracy of the low-temperature phase to induce the same prefactor.
Least-squares fits to the data in Table 1 according to the laws in Section 2 have been conducted. We have left out the smaller lattices systematically for each fit, until a goodness-of-fit value of at least Q = 0.5 was found for each observable individually. The same protocol was employed earlier [7] for a reduced time series, where only every L 3 -th measurement was used. There we were not challenging the prefactors of higher-order corrections so the reduced time series was sufficient. We list all the fit parameters obtained for both the full time series and the reduced one in Table 2 along with the quality-of-fit parameters Q and the degrees of freedom left. The inverse transition temperatures in the thermodynamic limit are effectively identical and do not depend on whether we use the reduced or the full dataset or on the precise final averaging procedure. A graphical representation of the best fits is given in Fig. 2. Table 1 Resulting quantities for the gonihedric Ising model (2) with periodic boundary conditions: extremal values for the specific heat C max V , Binder's energy parameter B min , with their respective pseudo-critical inverse temperatures β, and inverse temperatures where peaks of the energy probability density have equal heights and weights for finite lattices with linear length L. The finite-lattice interface tension is listed as σ , the energy of the ordered and disordered phases and their difference as e o , e d and e. The infinite lattice limits are listed as parameters of fits whose goodness-of-fit value Q is given. Highlighted in light grey are the data points used for fits with only leading order correction (lo). Additional data points used for fitting with subleading corrections (so) up to and including order O(L −4 ) are marked in dark grey.  (7) 0.88265 (7)  11 0.532894(13) 73.766(7) 0.531286(13) 0.27998 (7) 0.532897 (13) (7) 513.571 (12) 0.545821 (7) 0.326601(29) 0.546044 (7) 0.545952 (7) 0.104440 (20) Table 2 Resulting parameters of the best fits to the extremal points β for the specific heat C max V , Binder's energy parameter B min ; or inverse temperatures β eqw and β eqh to laws of the form β(L) = β ∞ + p 2 /L 2 + p 3 /L 3 + p 4 /L 4 . Parameters p i not used in the specific fit are marked with -. The error-weighted average over all four inverse temperatures are listed as β av , whereas β av w/o eqw is the average, where β eqw is left out because it is strongly correlated with β C max V and would effectively weight this value twice.  (7) 1.313493 (12) full time series, linear fit  Since the inverse temperatures β eqw and β C max V are obviously strongly correlated, we leave out the former and average over β C max V , β B min , and β eqh , neglecting cross-correlations [28] between those and taking the smallest contributing error for the final estimate. Our best estimate of the inverse transition temperature is then given by which accounts for the higher-order scaling corrections up to O(L −4 ).
Although the inverse transition temperatures do not change, we employ the full data set. The reason is that the error on β eqh becomes smaller for the time series that uses the full, correlated data set. This is attributed to the fact that the observable relies on the statistics in single bins of the energy histogram, which in total becomes smoother when using more, correlated measurements. The same argument is valid for the calculation of the interface tension (33), for which the best fit with corrections of order O(L −2 ) yields a value of Table 3 Quality and resulting quantities of the canonical simulations. In the time series of n meas = 1 048 576 measurements we found autocorrelation times τ < 7 (in units of measurements), leading to approximately n uncorrelated measurements. The energy ê o/d and the specific heat Ĉ o/d in the different phases have been measured by preparing two independent systems for each respective phase at inverse temperature β = 0.5513. We give the error-weighted average over lattices L ≥ 32, where the dependence on the lattice size is smaller than the statistical error. The last line gives infinite-volume limits from the multicanonical simulations for comparison.
Moments of the energy in the pure ordered and disordered phases are also expected to become more accurate using the full data set, since autocorrelation times in the pure phases are then significantly smaller than τ int for the full energy range (see below). By fitting the scaling law (27) to these moments, one obtains the latent heat in the infinite-volume limit, ê = 0.85148(5) original model, periodic bc.
Taking a careful look at the scaling laws in Section 2, we find that the prefactors of the scaling corrections only depend on the moments of the energy or their differences. We have two methods at hand to test the self-consistency of our simulations. Firstly, since the statistics of the observables are very high, we can retrieve the prefactors of the corrections as parameters of (nonlinear) least-squares-fits with all corrections up to and including O(L −4 ). Secondly, from multicanonical simulations we get estimators (30) of the energy moments, allowing a direct computation of those prefactors.
In addition, we carried out independent canonical simulations for the original model under periodic boundary conditions for very large lattices. The goal was to get independent measurements of the moments of the energy in the ordered and disordered phases. Here we prepared the system in the appropriate phase and performed the simulations at a fixed temperature β = 0.5513, near the transition temperature, exploiting the fact that in canonical simulations, for large lattices, flips between the two phases are extremely unlikely. Of course, this was only possible after having determined the transition temperature with high accuracy by the multicanonical simulations.
The quality of the canonical measurements and estimators on the energy and the specific heat are summarized in Table 3. The autocorrelation times within the phases are significantly smaller, because the system does not traverse suppressed, improbable states between the phases. The statistical error has again been retrieved by jackknife analysis. For lattices with size L ≥ 32, physical quantities indicate no further dependence on the lattice size within the statistical error. Therefore we can safely take the error-weighted averages over energy moments and their differences for those lattices. The multicanonical and the canonical estimates of energetic moments agree astonishingly well.
With use of the energy moments from both simulations, we can challenge the prefactors in the finite-size scaling laws numerically by comparing the numerical values for the fit parame- Table 4 Resulting prefactors of the finite-size scaling corrections of the original model, retrieved by fitting the ansatz, compared to direct calculations from estimators for the energy ê o/d and specific heat Ĉ o/d of the ordered and disordered phases. In the multicanonical simulations these moments were determined by finite-size scaling of e o/d (β eqw (L)), C o/d (β eqw (L)); and in the canonical case by measuring time series directly at β = 0.5513 β ∞ . (7) fit on β C max V (L) 0.8771 (14) 2.371 (4 (16) 0.34723 (9) ters to the theoretically expected prefactors in terms of the energy moments. The results of this cross-check are collected in Table 4. Employing the scaling relation for the specific-heat maximum (22), we can calculate ê from the fit parameter of the leading contribution. Using our estimate of β ∞ = 0.551332(8), we get ê = 0.85130 (7), very close to our estimate 0.850968 (18) from the moments of the canonical simulations. The leading correction to the specific-heat ansatz (22) has a prefactor which computes to 0.2197 (17) from the canonical moments. The fits find 0.17 (6), which is compatible, if not quite accurate.
The minimum of the Binder parameter (23) for the infinite lattice is found to be 0.34729(7) from the direct fit of our multicanonical data which agrees within error bars with 0.34723(9) from the canonical and 0.347(4) from the multicanonical energy moments. The first correction in Eq. (24) yields a value of −9.195 (14) when inserting the energy moments from the multicanonical simulation. The fits find −9.12(4) which is very close.
The coefficient of the leading correction apparent for all inverse temperatures, p 2 = 3 ln(2)/ ê, agrees reasonably well for the fits on β B min and β eqh (cf. Table 4). The fits on β eqw and β C max V yield a slope of 2.371(4) which within error bars is slightly off from the value of our best estimate of 2.44362(6) from the energy moments. The relative error between the two values is very small though, around 3%, which is acceptable given that the leading contribution probably accounts for the omitted higher-order contributions with different sign and that we have neglected all exponential corrections.
The second leading correction of order O(L −3 ) of β B min has a prefactor of the form 2 ln(ê o /ê d )/ ê, which we expect to have a value of 2.03625(6) from the energy moments. The corrections of fourth order to β B min , β C max V and β eqw are supposed to be identical from the analytical expansion in Section 2. The fits of the inverse temperatures β C max V and β eqw in Table 2 suggest a value around 17 for the O(L −4 ) contribution. For the lattice sizes accessible to the multicanonical approach, this is of the same absolute order of magnitude but with different sign compared to the third-order contribution. Therefore they should, in principle, compensate each other. This is reflected by the fact that the second-order contribution p 2 of β B min is closest to the expected one. In accordance, fitting the observable to the law β ∞ + p 2 /L 2 + p 3 /L 3 gives a prefactor p 3 = 0.65 (12) with the wrong sign compared to the theoretical prediction (20), compensating the next contribution. We therefore also looked at the fit including the fourth term of order O(L −4 ). Not taking the explicit values too seriously, we find p 3 = −1.6(8), p 4 = 13(4), which reflects qualitatively the compensation of those contributions for our lattice sizes at hand. Overall, we must conclude that least-squares fitting cannot be pushed any further given our statistics.
The observation that β C max V and β B min have the same O(L −4 ) contribution can be exploited (implicitly also making use of the cross-correlations) by looking at their difference. Here, we expect from (19) and (20) (7), in excellent agreement with 2.03625(6) from the formula, where the relative error between the two is less than 0.1%.
Finally, we can also compare the numerical values for the correction to the energies via (27). For ê o , we find a prefactor of 1.329(5) from the specific heat Ĉ o , compared to the value of 1.397(5), for ê d a value of 6.80(3) compared to 6.09(4), which is roughly 10% off.
The overall consistency of our results is very good, given that we neglected all exponential corrections. No estimates for the prefactors differ by more than 10%, and the various estimates of the inverse transition temperature are insensitive to the actual fitting protocol we use. This clearly demonstrates that the first correction terms are properly predicted by the simple two-state model even in the case of models with an exponential degeneracy of the low-temperature phase.
The earlier canonical Monte Carlo simulations of the original plaquette model yielded values of β ∞ = 0.50(1) [8] and more recently canonical simulations of the dual model (3) gave β ∞ = 0.510(2) [10]. Another previous estimate for the infinite-lattice inverse transition temperature, reported by Baig et al. [9] from canonical simulations using fixed boundary conditions, β ∞ = 0.54757(63), is much closer to the results here.
We therefore measured the inverse transition temperature again using multicanonical simulations for both the dual model (3) under periodic boundary conditions and the original model (2) with fixed boundary conditions. We resolve those inconsistencies, as we show in the following.

Dual model with periodic boundary conditions
For the dual model, we performed a number of n max = 4 × 10 6 sweeps and took measurements every sweep for even lattice sizes up to L = 24. The inverse temperatures of the dual model were fitted to laws with the leading correction of order O(L −2 ), which should be well covered by our data. The best fits on the inverse temperatures are shown in Fig. 3, where we used the data in Table 5 that also lists the other quantities of interest.
Since the inverse temperatures β eqw and β C max V are again obviously strongly correlated, we leave out the former and average over β C max V , β B min , and β eqh , neglecting cross-correlations [28] between those. We then find the error weighted average, for the inverse transition temperature. The error is again taken as the smallest error of the contributing estimates. The temperature β ∞ dual of the dual model is related to the temperature in the original model, β ∞ , by the duality transformation Applying standard error propagation, we retrieve a value of β ∞ = 0.55143(7) from duality, periodic bc (39) for the original model. This agrees very well with 0.551 332(8), obtained from the direct simulation, considering that higher-order corrections in the dual model are omitted and additional exponential corrections [15,16,18] in the finite-size scaling were ignored completely for both models.
We argue that the earlier estimates of the transition temperature were clearly hampered by strong hysteresis effects. Apart from the locations of the hysteresis branches being cooling-rate dependent, it is hard to estimate the transition temperature reliably from their locations. This is illustrated in Fig. 4, where the multicanonical data of the dual model is located between the two hysteresis branches. Such effects are very difficult to tackle using canonical Monte Carlo data, as already remarked by the authors of Ref. [10].
For the interface tension (33) of the dual model we find σ = 0.1214 (13). This value agrees very well with the interface tension of the original model, σ = 0.12037 (18), which raises interesting questions about the duality of the model.

Original plaquette model with fixed boundary conditions
The remaining open question is the difference between our inverse transition temperature compared to the value in the same model under fixed boundary conditions. In the thermodynamic limit, we expect a system to be independent of its boundary conditions, since the boundaries grow Table 5 Simulation results for the dual model (3): extremal values for the specific heat C max V , Binder's energy parameter B min , with their respective pseudo-critical inverse temperatures β, and inverse temperatures where peaks of the energy probability density have equal heights and weights for finite lattices with linear length L. The finite-lattice interface tension is listed as σ , the energy of the ordered and disordered phases and their difference as e o , e d and e. Light grey cells mark the values used for fitting, so that the goodness-of-fit parameter Q > 0.5. If Q < 0.5 for all fits, we took that one with the largest Q.  Fig. 10 of Ref. [10]. One can see that heating the system up (decreasing the inverse temperature β) or cooling it down (increasing β) can lead to strong hysteresis effects in the energy. Our multicanonical data lies in between both branches of the hysteresis curve, but not in the centre, as one may heuristically assume.
like a surface, whereas the system size grows with the volume. We therefore reinvestigated the model (2) using the multicanonical algorithm for such fixed boundary conditions. For our simulations we enclosed L 3 free spins in a larger cube with frozen outer planes, so the whole system contained (L + 2) 3 spins. Our results are listed in Table 6. We performed linear regression on the peak locations β(L) of the specific heat and Binder's parameter according to the law that was also used by Baig et al. [9], and fitted the inverse temperatures. The statistical errors of the constant a 2 turned out to be of the same order as the value itself, therefore we set a 2 = 0 for the fits in Table 6, intentionally neglecting the contribution O(L −2 ). The best fits and the energy probability density are shown in Fig. 5 and the weighted average of inverse transition temperatures is given by: This estimate of the inverse transition temperature is thus in excellent agreement with the other results obtained here from multicanonical simulations with periodic boundary conditions for the gonihedric Ising model β ∞ = 0.551332 (8) and the dual model β ∞ = 0.55143 (7).
Direct comparison to Ref. [9] shows that while inverse transition temperatures are reproduced, the extremal values of observables are not. The following observations may help to clarify the deviations. The authors of Ref. [9] simulated the system by applying periodic boundary conditions and fixing one plane parallel to the xy-, yzand zx-planes each. If their simulation box consisted of a total number of L 3 spins, they simulated (L − 1) 3 free spins. Thus our data with lattices of linear length L has to be compared to their data with L = L + 1. Also, their specific heat and magnetic susceptibility χ = βL −d ( M 2 − M 2 ) have to be multiplied by a factor of (L + 1) 3 /L 3 to be comparable with our normalization, since these quantities are proportional to the inverse system volume. Here, M = i σ i is the total magnetization, which for fixed boundary conditions is a well defined order parameter. Table 6 Extremal points of the gonihedric Ising model (2) for the specific heat C max V , Binder's energy parameter B min , and the pseudo-critical inverse temperatures β for finite lattices with linear length L contained in a box with fixed boundary conditions. For each lattice size, a number of n prod = 10 6 measurements were taken. Errors have been calculated with jackknife error estimation using 20 blocks.  Table 7 Values for comparison with Ref. [9], with the hat denoting the observables as calculated by Baig et al. Linear lattice lengths are L = L + 1, with L being the number of free spins in each direction. The specific heat Ĉ max The inverse temperatures β are the same for their data and ours. Magnetic susceptibilities χ max = L 3 /(L + 1) 3 χ max are listed as well. Binder's energy parameter has no explicit volume-dependence by design, but it is sensitive to offsets in the energy, which cancel in the specific heat. Our values of the Binder parameter minima differ significantly from [9]. However, if we shift our measured energies E to get Ê = E − 1.5L 2 = E − 1.5(L + 1) 2 and calculate Binder's parameter (6) with Ê instead of E, our measurements reproduce those of [9] very well. The energy E of the system can be written in terms of the number of plaquettes with an even or odd number of parallel aligned spins, n + or n − , Since we measure the same cumulant values for shifted energies Ê , in Ref. [9] an additional number of n + = n + + 3L 2 plaquettes contribute to the system's energy because energetic contributions from the fixed planes, where all spins are aligned, were included. For direct comparison, the resulting quantities are listed in Table 7 after applying all corrections, showing that our data is then in very good agreement with Ref. [9]. For completeness we include here also our data for the magnetic susceptibility χ. The deviation from the fitting results  (4) 0.0281 (7) in Ref. [9] simply stems from the fact that our simulations are performed with the multicanonical method that allows a finite-size scaling analysis with more and significantly larger lattice sizes. The interface tension (33) as a function of the linear lattice size was also extracted and its infinite-volume limit yields a value of σ = 0.0281 (7) where we allowed corrections of order O(L −2 ) in the fits. Note that this value is about four times smaller than that for the same model with periodic boundary conditions, which may point to unexpectedly large finite-size effects for this quantity.

Summary
We simulated the plaquette gonihedric Ising model and its dual to shed some light on discrepancies in the recent literature on the reported value(s) of the first-order phase transition temperature. High-precision results from multicanonical simulations forced us to review the traditional scaling ansatz for first-order finite-size corrections by taking the exponential low-temperature phase degeneracy of the model into account. The leading correction in such circumstances is then no longer proportional to the inverse volume of the system, O(L −3 ), but is rather O(L −2 ). With this finite-size scaling ansatz, our simulations with periodic boundary conditions produced consistent results for both the original formulation of the model as well as its dual representation. Since our results also differed from earlier simulations using fixed boundary conditions, where the leading corrections are now O(L −1 ), we carried out multicanonical simulations of the gonihedric Ising model with these boundary conditions too. The resulting inverse transition temperature was fully consistent with the value found using periodic boundary conditions when larger lattices were included, and hopefully settle once and for all enduring inconsistencies.
The main resulting physical quantities that characterize the first-order phase transition are summarized in Table 8 for the different models and boundary conditions in our simulations. We find an overall consistent value for the inverse transition temperature of β ∞ = 0.551334 (8) ( and we measure the interface tension of the original model and its dual for the first time. We find values of σ = 0.12037 (18) and σ = 0.1214 (13) for the original and dual model with periodic boundary conditions, respectively. The interface tension of the original model with fixed boundary conditions is found to be much smaller, σ = 0.0281 (7). Interestingly, we also estimate different values for the latent heat with different boundary conditions, ê = 0.694(4) in the case of fixed boundary conditions compared to ê = 0.850968(18) for periodic boundaries.
That the latent heat may apparently depend on the boundary conditions has been observed earlier in simulations of the q-state Potts model [19]. We suspect, however, that the estimates for fixed boundary conditions are flawed by extreme finite-size effects which are difficult to quantify. Any model with an exponentially degenerate low-temperature phase will display the modified scaling at a first-order phase transition we have delineated for the three-dimensional gonihedric model and its dual here. Apart from higher-dimensional variants of the gonihedric model or its dual, there are numerous other models in which the scenario could be realized. Examples range from ANNNI models [29] to topological "orbital" models in the context of quantum computing [30] which all share an extensive ground-state degeneracy. Among the orbital models for transition metal compounds, a particularly promising candidate is the three-dimensional classical compass or t 2g orbital model where a highly degenerate ground state is well known and the signature of a first-order transition into the disordered phase has recently been found numerically [31].
Other systems, such as the three-dimensional Ising antiferromagnet on an FCC lattice, have an exponentially degenerate number of ground states but a small number (6 in the case of the FCC Ising antiferromagnet) of true low-temperature phases. Nonetheless, they do possess an exponentially degenerate number of low-energy excitations so, depending on the nature of the growth of energy barriers with system size, an effective modified scaling could still be seen at a first-order transition for the lattice sizes accessible in typical simulations. The crossover to the true asymptotic (standard) scaling would then only appear for very large lattices. Indeed, previous simulations appear to have found nonstandard scaling for the first-order transition in the three-dimensional Ising antiferromagnet on an FCC lattice [32].