Anisotropy of the interface tension of the three-dimensional Ising model

We determine the interface tension for the 100, 110 and 111 interface of the simple cubic Ising model with nearest-neighbour interaction using novel simulation methods. To overcome the droplet/strip transition and the droplet nucleation barrier we use a newly developed combination of the multimagnetic algorithm with the parallel tempering method. We investigate a large range of inverse temperatures to study the anisotropy of the interface tension in detail.


Introduction
In many physical systems with discrete symmetry the anisotropy of the interface tension can play an important role for various phenomena, including equilibrium droplet shapes [1] and the interfacial roughening transition [2]. For sufficiently strong anisotropy, facets, edges, or even corners can be identified in the equilibrium droplet shape. Due to the anisotropy of the interface tension in the threedimensional (3D) Ising model, the shape of the equilibrium droplet at some finite temperature is not spherical and has, in principle, to be determined by the Wulff construction [3]. Since the 3D Ising model with nearest-neighbour interaction is not exactly solvable, no analytical results are available for the interfacial free energy and the Wulff construction can only be done using an effective model of the angle-depending interface tension. Only for temperatures not too far below the critical temperature one can use the spherical approximation and, therefore, it is important to know how large the anisotropy is for a given temperature. Whereas a lot of numerical results are available for the planar 100 interface tension of the simple cubic Ising model, see e.g. Refs. [4,5,6,7,8,9,10], there are only a few results in the literature for the 110 interface [11] and, to our knowledge, there are no results at all for the 111 interface.
The layout of the remainder of this paper is organized as follows. In Sec. 2 we discuss the results for the planar 100 interface and compare our results with previous estimates in the literature. Next we describe in Sec. 3 first the special boundary conditions employed for the simulation of the tilted interfaces with 110 and 111 orientation and then discuss the results of our finite-size scaling analysis. Finally, in Sec. 4 we conclude with a summary of our main findings.

Planar interface
We considered the Ising model on L × L × L simple cubic lattices with periodic boundary conditions in all three directions to simulate systems with planar 100 interfaces, for various temperatures below the Ising transition at β c ≡ 1/T c = 0.22165459 [12]. For a typical configuration see Fig. 1.
The interface tension can be measured using a multimagnetical (flat in the distribution of the magnetization m) simulation combined with parallel tempering [13], the result of which is after appropriate reweighting to the canonical ensemble a double-peaked magnetization density P (m). We simulated n = 26 replica of the system at different inverse temperatures β i , with β i = 0.195 + 0.005(i − 1) and i = 1, . . . , n. A planar interface of the 3D Ising model exhibits a transition at the roughening temperature [2] T R = 1/β R , with β R = 0.40758(1), above which the surface stiffness for the 100 interface is finite and below which it is infinite. Therefore, we restrict ourselves to the temperature range above this transition, i.e. T i = 1/β i > T R for all i . To construct the weight function for the multimagnetic part of the algorithm, we employed an accumulative recursion, described in detail in Refs. [14] and [15]. Statistical averages were taken over runs of 1 × 10 6 Monte Carlo (MC) steps, where one MC step consists of one full multimagnetical lattice sweep for all 26 replica and one attempted parallel tempering exchange of all adjacent replica.
The interface tension σ 100 can be estimated according to Refs. [16,4] by (d = 3) where P (L) min is the value of the magnetization density in the mixed phase region m ≈ 0 (strip phase) and P (L) max the value at its maxima located close to the equilibrium magnetization m = ±m 0 . Therefore, the finite-volume estimator is given by with c 1 = (ln A)/2β.
The power of L in the prefactor of Eq. (1) is a delicate problem and the knowledge of the pre-exponential behaviour fixes one free parameter of the fit. Using the capillary wave approximation [17,18,19] the exponent x = (d − 3)/2, i.e. x = 0 for d = 3. Therefore, we performed finite-size scaling fits according to and, to allow for higher-order corrections, also to We performed simulations for various lattice sizes ranging from L = 4 to L = 26. In Fig. 2 we show the magnetization density P (m) for β = 0.3, where the strip configurations, corresponding to the minimum between the two peaks, are suppressed by more than 175 orders of magnitude for the largest system and this suppression becomes even more pronounced for lower temperatures. Such unlikely configurations would not be a problem for a multimagnetical algorithm, but between the strip configuration and the droplet configuration there is an exponentially large barrier [20,21] that might not be overcome during the equilibration phase. Therefore, it is necessary to use the combined algorithm to overcome this barrier. A similar reasoning applies to the evaporation/condensation transition which is another hidden albeit weaker barrier in the multimagnetical system [21,22,23].
For every system the maximum and minimum probability P (L) max and P (L) min were read off, and by repeating the simulations 32 times the statistical error bars were obtained. For β = 0.3 and L ≥ 12 the resulting values for σ 100 (L) are plotted in Fig. 3. To check the stability of the fit results we performed fits with different lower bounds L min of the fit range. The upper bound of the fits was always the largest lattice L = 26. For the fits according to Eq. (4) we find due to the systematic variation of the lower bound a trend to larger values of σ 100 with increasing L min . This can be seen in the left panel of Fig. 4 where we also include the goodness-offit parameter Q into the figure to judge the quality of the fits. Above L min = 14 the goodness-of-fit parameter was well above 0.05, which we chose as cutoff value. Nevertheless, not yet reaching a constant value for σ 100 led us to include one more parameter in our fits, Performing fits with this ansatz and again varying the lower fit bound systematically, the resulting values for the planar interface tension σ 100 stay almost constant for reasonable fits (Q ≥ 0.05), as one can see in the right panel of Fig. 4. The infinite system size extrapolation in 1/L 2 according to Eq. (5) with L min = 12 yields for the particular inverse temperature β = 0.3 a value of σ 100 = 1.00933 (12) for the planar interface tension with goodness-of-fit parameter Q = 0.18, which is in good agreement with the result from Hasenbusch and Pinn [6] σ 100 = 1.009302(106). All results with Q ≥ 0.05 for the three different infinite system size extrapolations of the planar interface tension are collected in Tables 1 and 2.

Tilted interfaces
To generate a tilted interface we used a 2L × L × L simple cubic lattice under two sets of boundary conditions. We chose a rather unusual combination of boundary conditions for this study.    Table 1. Results for the planar interface tension σ 100 using fits according to Eqs. (3), (4), and (5), respectively. In the last column we include for comparison the results of Hasenbusch and Pinn [6]. fit ansatz (3) fit ansatz (4) fit ansatz (5) Ref. [6] β fit range Q σ 100 fit range Q σ 100 fit range Q σ  we imposed periodic boundary conditions in the x and y direction and shifted boundary condition in the z direction. To be more precise, the neighbour in negative z direction of a spin in the first xy layer of the system s(x, y, 1) is s(x+L, y, L) for 1 ≤ x ≤ L and s(x, y, L) for L < x ≤ 2L. A typical configuration for such a system below the Ising transition is depicted in Fig. 5.
For the lattice with an interface along the 111 direction we imposed periodic boundary conditions in the x direction and shifted boundary condition in the y and z direction. Therefore, the neighbour in negative y direction of a spin in the first xz layer of the system s(x, 1, z) is s(x + L, L, z) for 1 ≤ x ≤ L and s(x, L, z) for L < x ≤ 2L. A typical configuration for such a system below the Ising transition is depicted in Fig. 6.
Using the same setup as described above we measured the probability density of the magnetization for various lattice sizes ranging from L = 4 to L = 20. From these distributions, we determined the interface tension in the 110 and 111 direction by means of infinite-system size extrapolations via Eqs.  Tables 3, 4, 5, and 6, respectively.
For σ 111 the difference between the final estimates for the interface tension using the different infinite-volume extrapolations is larger than for the other two interface directions. This is an indication that the finite-size effects are more pronounced in this setup. One reason for this is that the distance between the two interfaces is too small and therefore the fluctuations of the interfaces are correlated and the effect becomes more pronounced as the critical temperature is approached. In spite of the computational effort, the system sizes are still too small to give equally accurate values for σ 111 as for σ 100 and σ 110 , respectively. Nevertheless, the statistical error of the interface tension of the 111 interface is only roughly 0.01% − 0.04%, which is more than one order of magnitude smaller than the effect which we are interested in, namely the anisotropy of the interface tension.
Using the results from the fits according to Eq. (3) we determined the anisotropy of the interface tension for different directions, namely the 110 and 111 direction. Due to the large finite-size effects near the critical point, the accessible temperature range is limited to a window at low temperatures. Therefore, we cannot give reliable results for the region near the transition. In Fig. 9 we show the anisotropy as a function of the reduced temperature and include for comparison results for the twodimensional Ising model. In the two-dimensional case we calculated the anisotropy using the exact expressions for the interface tension of the 10 "surface" [24], and for the 11 "surface" [25], Our results show that the anisotropy of the interface tension as a function of the reduced temperature grows faster in three dimensions than in two. Nevertheless, the absolute value is still very small (< 3%) for the temperature range investigated in this work.

Summary
We have presented a careful analysis of the interface tension in 100, 110 and 111 direction of the simple cubic Ising model with nearest-neighbour interaction. Using a newly developed combination of the multimagnetic algorithm with the parallel tempering method, we were able to measure the highly suppressed configurations of the strip phase for systems up to 26 3 at β = 0.32, i.e. T ∼ 0.7 T c . We show that at given T /T c the anisotropy of the interface in three dimensions is larger than in two dimensions. However, down to 0.7 T c it never exceeds 3%, so that in most cases the isotropic approximation for droplet condensation phenomena should be sufficiently accurate.

Acknowledgements
We wish to thank Martin Hasenbusch for helpful discussions.   Fig. 9. The anisotropy of the interface tension for the 110 and 111 direction as a function of the temperature T = 1/β (upper scale), respectively. When approaching the 3D Ising transition temperature indicated by the arrow, the interface tension becomes isotropic and finally vanishes. To allow a comparison with the 2D Ising model, the lower scale shows the reduced temperature T /T c = β c /β. The inset shows another comparison, where the abscissa is proportional to the asymptotic correlation length ξ ∝ |1 − T /T c | −ν (with ν = 1 in 2D and ν = 0.63 in 3D). Table 3. Results for the 110 interface tension using of fits according to Eqs.