Planar ordering in the plaquette-only gonihedric Ising model

In this paper we conduct a careful multicanonical simulation of the isotropic $3d$ plaquette ("gonihedric") Ising model and confirm that a planar, fuki-nuke type order characterises the low-temperature phase of the model. From consideration of the anisotropic limit of the model we define a class of order parameters which can distinguish the low- and high-temperature phases in both the anisotropic and isotropic cases. We also verify the recently voiced suspicion that the order parameter like behaviour of the standard magnetic susceptibility $\chi_m$ seen in previous Metropolis simulations was an artefact of the algorithm failing to explore the phase space of the macroscopically degenerate low-temperature phase. $\chi_m$ is therefore not a suitable order parameter for the model.


Introduction
The 3d plaquette ("gonihedric") Ising Hamiltonian displays an unusual planar flip symmetry, leading to an exponentially degenerate low-temperature phase and non-standard scaling at its first-order phase transition point [1,2]. The nature of the order parameter for the plaquette Hamiltonian has not been fully clarified, although simulations using a standard Metropolis algorithm [3] have indicated that the magnetic ordering remains a fuki-nuke type planar layered ordering, which can be shown rigorously to occur in the extreme anisotropic limit when the plaquette coupling in one direction is set to zero [4,5,6] by mapping the model onto a stack of 2d Ising models.
The simple 3d plaquette Ising Hamiltonian, where the Ising spins σ i = ±1 reside on the vertices of a 3d cubic lattice, can be considered as the κ = 0 limit of a family of 3d "gonihedric" Ising Hamiltonians [7], which contain nearest neighbour i, j , next-to-nearest neighbour i, j and plaquette interactions [i, j, k, l], j,k,l] σ i σ j σ k σ l .
For κ = 0 parallel, non-intersecting planes of spins may be flipped in the ground state at zero energy cost, leading to a 3 × 2 2L ground-state degeneracy on an L × L × L cubic lattice which is broken at finite temperature [8]. For the κ = 0 plaquette Hamiltonian of Eq. (1), on the other hand, the planar flip symmetry persists throughout the low-temperature phase and extends to intersecting planes of spins. This results in a macroscopic low-temperature phase degeneracy of 2 3L and non-standard corrections to finite-size scaling at the first-order transition displayed by the model [1,2]. The planar flip symmetry of the low-temperature phase of the plaquette Hamiltonian is intermediate between the global Z 2 symmetry of the nearest-neighbour Ising model and the local gauge symmetry of a Z 2 lattice gauge theory H gauge = − [i,j,k,l] U ij U jk U kl U li (4) and naturally poses the question of how to define a magnetic order parameter that is sensitive to the first-order transition in the model. The standard magnetization will clearly remain zero with periodic boundary conditions, even at lower temperatures, because of the freedom to flip arbitrary planes of spins. Similarly, the absence of a local gauge-like symmetry means that observing the behaviour of Wilson-loop type observables, as in a gauge theory, is also not appropriate.

Fuki-Nuke Like Order Parameters
Following the earlier work of Suzuki et al. [4,5], a suggestion for the correct choice of the order parameter for the isotropic plaquette Hamiltonian comes from consideration of the J z = 0 limit of an anisotropic plaquette model where we now indicate each site and directional sum explicitly, assuming we are on a cubic L × L × L lattice with periodic boundary conditions σ L+1,y,z = σ 1,y,z , σ x,L+1,z = σ x,1,z , σ x,y,L+1 = σ x,y,1 . This will prove to be convenient in the sequel when discussing candidate order parameters. When J z = 0 the horizontal, "ceiling" plaquettes have zero coupling, which Hashizume and Suzuki denoted the "fuki-nuke" ("no-ceiling" in Japanese) model [5]. The anisotropic 3d plaquette Hamiltonian at J z = 0, may be rewritten as a stack of 2d nearest-neighbour Ising models by defining bond spin variables τ x,y,z = σ x,y,z σ x,y,z+1 at each vertical lattice bond which satisfy L z=1 τ x,y,z = 1 trivially for periodic boundary conditions. The τ and σ spins are related by an inverse relation of the form σ x,y,z = σ x,y,1 × τ x,y,1 τ x,y,2 τ x,y,3 · · · τ x,y,z−1 , and the partition function acquires an additional factor of 2 L×L arising from the transformation. The resulting Hamiltonian with J x = J y = 1 is then simply that of a stack of decoupled 2d Ising layers with the standard nearest-neighbour in-layer interactions in the horizontal planes, (τ x,y,z τ x+1,y,z + τ x,y,z τ x,y+1,z ) , apart from the L 2 constraints L z=1 τ x,y,z = 1 whose contribution should vanish in the thermodynamic limit [6]. Each 2d Ising layer in Eq. (9) will magnetize independently at the 2d Ising model transition temperature. A suitable order parameter in a single layer is the standard Ising magnetization which when translated back to the original σ spins gives which will behave as ± |β − β c | 1 8 near the critical point β c = 1 2 ln(1 + √ 2). More generally, since the different τ x,y,z layers are decoupled in the vertical direction we could define Two possible options for constructing a pseudo-3d order parameter suggest themselves in the fuki-nuke case. One is to take the absolute value of the magnetization in each independent layer the other is to square the magnetization of each plane, to avoid inter-plane cancellations when the contributions from each Ising layer are summed up. We have explicitly retained the various normalizing factors in Eqs. (13) and (14) for a cubic lattice with L 3 sites. The suggestion in [3,5] was that similar order parameters could still be viable for the isotropic plaquette action, namely where we again assume periodic boundary conditions. Similarly, for the case of the squared magnetizations one can define with obvious analogous definitions for the other directions, m y sq and m z sq , which also appear to be viable candidate order parameters. In the isotropic case the system should be agnostic to the direction so we would expect m x abs = m y abs = m z abs and similarly for the squared quantities. The possibility of using an order parameter akin to the m sq in Eq. (14) had also been suggested by Lipowski [9], who confirmed that it appeared to possess the correct behaviour in a small simulation.
In Ref. [3] Metropolis simulations gave a strong indication that m x,y,z abs and m x,y,z sq as defined above were indeed suitable order parameters for the isotropic plaquette model, but these were subject to the difficulties of simulating a strong first-order phase transition with such techniques and also produced the possibly spurious result that the standard magnetic susceptibility χ behaved like an order parameter. The aforementioned difficulties precluded a serious scaling analysis of the behaviour of the order parameter with Metropolis simulations, including an accurate estimation of the transition point via this route.
With the use of the multicanonical Monte Carlo algorithm [10,11] in which rare states lying between the disordered and ordered phases in the energy histogram are promoted artificially to decrease autocorrelation times and allow more rapid oscillations between ordered and disordered phases, combined with reweighting techniques [12], we are able to carry out much more accurate measurements of m x,y,z abs and m x,y,z sq . This allows us to confirm the suitability of the proposed order parameters and to examine their scaling properties near the first-order transition point.
The results presented here can also be regarded as the magnetic counterpart of the highaccuracy investigation of the scaling of energetic quantities (such as the energy, specific heat and Binder's energetic parameter) for the plaquette-only gonihedric Ising model and its dual carried out in [2] and provide further confirmation of the estimates of the critical temperature determined there, along with the observed non-standard finite-size scaling.

Simulation Results
We now discuss in detail our measurements of the proposed fuki-nuke observables [3] defined above, using multicanonical simulation techniques. The algorithm used is a twostep process, where we iteratively improve guesses to an a priori unknown weight function W (E) for configurations {σ} with system energy E = H({σ}) which replaces the Boltzmann weights e −βE in the acceptance rate of traditional Metropolis Monte Carlo simulations. In the first step the weights are adjusted so that the transition probabilities between configurations with different energies become constant, giving a flat energy histogram [13]. The second step is the actual production run using the fixed weights produced iteratively in step one. This yields the time series of the energy, magnetization and the two different fukinuke observables m sq and m abs in their three different spatial orientations. With sufficient statistics such time series together with the weights can provide the 8-dimensional density of states Ω(E, m, m x abs , . . .), or by taking the logarithm the coarse-grained free-energy landscape, by simply counting the occurrences of E, m, m x abs , . . . and weighting them with the inverse W −1 (E) of the weights fixed prior to the production run. Practically, estimators of the microcanonical expectation values of observables are used, where higher dimensions are integrated over in favour of reducing the amount of statistics required.
Although the actual production run consisted of N = (100−1000)×10 6 sweeps depending on the lattice sizes and is therefore quite long, the statistics for the fuki-nuke order parameters is sparser. For these, we carried out measurements every V = L 3 sweeps, because one has to traverse the lattice once to measure the order parameters in all spatial orientations and this has a considerable impact on simulation times. With skipping intermediate sweeps we end up with less statistics, but the resulting measurements are less correlated in the final time series. More details of the statistics of the simulations are given in [2].
In Fig. 1 we show the full time series for the magnetization m and the fuki-nuke parameter m x abs of the multicanonical measurements for an intermediate (L = 20) and the largest (L = 27) lattice size in the simulations along with the system energy. The time series of the energy per system volume e = E/L 3 of the larger lattice can be seen to be reflected numerous times at e −0.9 (coming from the disordered phase) and e −1.2 (coming from the ordered phase) which shows qualitatively that additional, athermal and non-trivial barriers may be apparent in the system [14]. It is clear that the standard magnetization is not a suitable order parameter, since it continues to fluctuate around zero even though the system transits many times between ordered and disordered phases in the course of the simulation. The fuki-nuke order parameter, on the other hand, shows a clear signal for the transition, tracking the jumps which are visible in the energy time series. This is also reflected in Fig. 2, where we show the estimators of the microcanonical expectation values with Ω i being the occurrences of pairs (E, O) in a reduced time series where the i-th block of length b has been omitted. The variance of the Jackknife estimators then is proportional to the squared statistical error of their mean where the prefactor (N B − 1)/N B accounts for the trivial correlations caused by reusing each data point in N B − 1 estimators [15]. Aside from reducing a systematic bias in derived quantities, the Jackknife error automatically takes care of temporal correlations as long as the block length b is greater than the autocorrelation time which was measured and discussed in great detail in [2]. We additionally confirmed that N B = 10 and N B = 40 for selected lattice sizes yield the same magnitude for the statistical error.
That the fuki-nuke parameters are, indeed, capable of distinguishing ordered and disordered states is depicted in Figs. 1 and 2. In the microcanonical picture we can clearly see that the different orientations of the fuki-nuke parameters are equal for the isotropic gonihedric Ising model, which we show for L = 20 in Fig. 3(a). This confirms that the sampling is at least consistent in the simulation. We also collect the microcanonical estimators for m x abs for several lattice sizes in Fig. 3(b), where a region of approximately linear increase between e −0.9 and e −1.3 can be seen for the larger lattices. This interval corresponds to the energies of the transitional, unlikely states between the ordered and disordered phases. Plaquettes successively become satisfied towards the ordered phase and thus the estimators that measure intra-and inter-planar correlations must increase, too.
As we stored the full time series along with its weight function, we are able to measure the microcanonical estimators for arbitrary functions of the measured observables f (O), which can be exploited to give a convenient way of calculating higher-order moments as well.
For canonical simulations reweighting techniques [12] allow the inference of system properties in a narrow range around the simulation temperature. That range and the accuracy are then determined by the available statistics of the typical configurations for the temperature of interest. Since multicanonical simulations yield histograms with statistics covering a broad range of energies, which is their most appealing feature and common to flat-histogram techniques, it is possible to reweight to a broad range of temperatures. The canonical Figure 4: Canonical curves for the fuki-nuke parameters m x abs and m x sq over a broad range of inverse temperature β for several lattice sizes L (compare with the canonical data in [3]). estimator at finite inverse temperature β > 0 is thus obtained by and again Jackknife error analysis is employed for an estimate of the statistical error, where we form individual Jackknife estimators by inserting the i-th microcanonical estimator of Eq. (18) into Eq. (21) and apply an analogue to Eq. (19). Since the microcanonical estimators for different orientations agree within error bars the canonical values will also be the same. Therefore, in Fig. 4 we only show one orientation for the two different fuki-nuke parameters. The overall behaviour of m x abs and m x sq from the Metropolis simulations of Ref. [3] is recaptured by the multicanonical data here. Namely, sharp jumps are found near the inverse transition temperature, as expected for an order parameter at a first-order phase transition. The transition temperature that was determined in the earlier simulations where energetic observables were measured under different boundary conditions and from a duality relation was β = 0.551 334(8) [2]. The positions of the jumps seen here (and in the earlier simulations) depend on the lattice size, and finite-size scaling can be applied to estimate the transition temperature under the assumption that the fuki-nuke parameters are indeed suitable order parameters.
To carry out such a finite-size scaling analysis, it is advantageous to look at the canonical curves of the susceptibilities, , since their peak positions provide an accurate measure of the finite-lattice inverse transition temperature β χ O (L). As an example Fig. 5 shows the peaks of the susceptibilities belonging to m x abs and m x sq for several lattice sizes. Qualitatively the behaviour of both susceptibilities is similar and as for the specific heat [2,16] their maxima scale proportional to the system volume L 3 but they differ in their magnitudes.
Empirically, the peak locations for the different lattice sizes L can be fitted according to the modified first-order scaling laws appropriate for macroscopically degenerate systems discussed in detail in [1,2], with free parameters a, b for the available 24 lattice sizes. Smaller lattices are systematically omitted until a fit with quality-of-fit parameter Q bigger than 0.5 is found. This gives for the estimate of the inverse critical temperature β χ (L) from the fuki-nuke susceptibility χ m with a goodness-of-fit parameter Q = 0.64 and 12 degrees of freedom left. Fits to the other directions m y,z abs and fits to the peak location of the susceptibilities of m x,y,z sq give the same parameters within error bars and are of comparable quality. The estimate of the phase transition temperature obtained here from the finite-size scaling of the fuki-nuke order parameter(s), β ∞ = 0.551 37 (3), is thus in good agreement with the earlier estimate β ∞ = 0.551 334 (8) reported in [2] using fits to the peak location of Binder's energy cumulant, the specific heat and the value of β where the energy probability density has two peaks of the same height or same weight. Interestingly, we find that the value of the coefficient for the leading correction also coincides. Assuming that the coefficient a = −2.46(3) of the leading correction is related to the inverse latent heat by a = −3 ln(2)/∆ê, as with the previous estimates [2], we find from Eq. (23) that ∆ê = 0.845 (8), in good agreement with the latent heat ∆ê = 0.850 968 (18) reported earlier in [2]. For visual confirmation of the finite-size scaling, the fuki-nuke magnetizations m x abs,sq along with their susceptibilities divided by the system volume are plotted in Fig. 6 by shifting the x-axis according to the scaling law, incorporating the fit parameters. The peak locations of the susceptibilities then all fall on the inverse transition temperature.
In the earlier work [3] the susceptibility χ m of the standard magnetization m unexpectedly behaved like an order parameter and it continues to behave idiosyncratically in the multicanonical simulations, but in a different manner. For compatibility with [3], the magnetic susceptibility divided by the inverse temperature, χ m = L 3 ( m 2 − m 2 ), is plotted in Fig. 7 along with the standard magnetization on a very small vertical scale (note that m should be between −1 and +1). χ m = 1 in the high-temperature phase but for the ordered, low-temperature phase the error rapidly increases below the transition temperature, though it is clear that the susceptibility is non-zero in this case too. Since m = 0, the behaviour of m 2 can provide insight into this behaviour of the susceptibility. Above the transition temperature in the high-temperature phase the sum over the free spin variables behaves like a random walk with unit step-size, therefore the expectation value of the squared total magnetization is given by M 2 = L 3 . Taking the normalization m = M/L 3 into account gives χ m = 1 in this region, as seen in Fig. 7.
Below the transition temperature, it is plausible that simulations in general get trapped in the vicinity of one of the degenerate low-temperature phases, each of which will have a different magnetization. A canonical simulation cannot overcome the huge barriers in the system and "freezes" with the same magnetization that the system had when entering the ordered phase. This accounts for the zero variance seen in the Metropolis simulations of [3]  below the transition temperature, since m is frozen. In multicanonical simulations, on the other hand, the system travels between ordered and disordered phases, thus picking one of the possible magnetizations each time it transits to an ordered phase which it sticks with until it decorrelates again in the disordered phase. Therefore, what is seen in the low-temperature region of Fig. 7 for χ m is that the variance of m is taking on rather arbitrary values due to the low statistics obtained compared to the large number, q = 2 3L , of degenerate phases one would have to visit to sample m 2 properly. Even with multicanonical simulations it is not possible to visit all of these macroscopically degenerate phases, and the increasing error bars reflects this. In the canonical case one gets stuck with one magnetization and would not notice the different values, leading to much more severe ergodicity problems in the finite simulation runs.
We investigate further the behaviour of the standard magnetization and susceptibility and fuki-nuke magnetizations in the model by preparing several fixed configurations with a given magnetization for a lattice with 10 3 spins and then only flipping complete planes of spins (a "flip-only" update), measuring the running average of the magnetization and fukinuke parameters. An example with the first thousand out of a total of 10 6 measurements is shown in Fig. 8 for three configurations picked at random from the ordered (e = −1.46), intermediate (e = −1.29) and disordered (e = −0.80) regions, respectively, along with the histograms for the magnetization in Fig. 9 obtained using the non-local flip-only update. It can be seen that whatever the initial value the running average of the (standard) magnetization becomes zero if one takes a long enough run so, as expected, the flip symmetry precludes a non-zero value. The fuki-nuke order parameters, on the other hand, should be invariant with respect to the plane-flip symmetry and this is clearly the case in Fig. 8(a), where the values of m x abs remain constant for the ordered, intermediate and disordered starting configurations.
The running average of the standard magnetic susceptibility χ m is plotted in a similar Running average of the standard magnetic susceptibility χ m plotted against the number t of plane-flips for the same three realizations. Note that the t-axis starts at 10 2 because χ m , being a variance, needs sufficiently many measurements to be meaningful. fashion in Fig. 8(b), where it can be seen that the ordered, intermediate and disordered configurations all converge after initial transients to values of χ m close to 1. The non-local plane-flips thus allow enough variability in the magnetization for the expected susceptibility value of χ m = 1 to be at least approximately attained even in the ordered phase, unlike the case of purely local spin flips. This suggests there might be some utility in incorporating such moves into a Metropolis simulation of the plaquette model to improve the ergodicity properties.
The histograms H(m) of the magnetization shown in Fig. 9 on a semi-logarithmic scale display interesting behaviour. They are symmetric around zero for all of the starting configurations because of the Z 2 symmetry of the Ising spins but they have rather different shapes in each case. The disordered starting configuration (e = −0.80) has a smooth maximum at m = 0, but both the intermediate (e = −1.29) and ordered (e = −1.46) starting configurations generate sharp peaks at m = 0. The pronounced peaks and valleys in the ordered histogram are presumably a consequence of the difficulty of reaching certain magnetization values (and the greater ease of reaching others) from an ordered starting configuration using only plane flips.
It would be interesting to construct and simulate a multimagnetic ensemble [17], where the weights give constant transition rates between configurations with different magnetizations to elucidate further on the magnetic and geometric barriers, as well as to confirm that m = 0 and χ = 1 for the low-temperature phase.

Conclusions
The multicanonical simulations presented here provide strong support for the idea that the plaquette gonihedric Ising model displays the same planar, fuki-nuke order seen in the strongly anisotropic limit of the model. In addition, the finite-size scaling analysis of the fuki-nuke order parameters gives scaling exponents in good agreement with the analysis of energetic quantities carried out in [2] and clearly displays the effect of the macroscopic lowtemperature phase degeneracy on the corrections to scaling. The transition temperature obtained here from the scaling of the fuki-nuke order parameters, and the amplitude for the leading correction to scaling term were found to be the same as those extracted from energetic observables. The analysis of the magnetic order parameters carried out here is thus complementary to the analysis of energetic observables in [2] and fully consistent with it.
The peculiar behaviour of the susceptibility of the standard magnetization χ m in the earlier Metropolis simulations of [3] was confirmed to be an artefact of the employed algorithm. However, even with multicanonical simulations, sampling the macroscopically degenerate low-temperature phase efficiently is difficult. Such problems could in principle be eased by introducing the plane-flips which provide a valid Monte Carlo update themselves and we investigated the behaviour of the standard magnetization and the fuki-nuke order parameters when such updates were applied.
With the work reported here on magnetic observables and the earlier multicanonical investigations of energetic quantities in [1,2] the equilibrium properties of the 3d plaquette gonihedric Ising model are now under good numerical control and the order parameter has been clearly identified. In the light of this clearer understanding, it would be worthwhile re-investigating non-equilibrium properties, in particular earlier suggestions [6,18,19] that the model might serve as a generic example of glassy behaviour, even in the absence of quenched disorder.