Phase transitions of antiferromagnetic Ising spins on the zigzag surface of an asymmetrical Husimi lattice

An asymmetrical two-dimensional Ising model with a zigzag surface, created by diagonally cutting a regular square lattice, has been developed to investigate the thermodynamics and phase transitions on surface by the methodology of recursive lattice, which we have previously applied to study polymers near a surface. The model retains the advantages of simple formulation and exact calculation of the conventional Bethe-like lattices. An antiferromagnetic Ising model is solved on the surface of this lattice to evaluate thermal properties such as free energy, energy density and entropy, from which we have successfully identified a first-order order–disorder transition other than the spontaneous magnetization, and a secondary transition on the supercooled state indicated by the Kauzmann paradox.


Introduction
The recursive lattice has been a classical methodology in statistical physics for several decades since its invention by Bethe [1] and Husimi [2]. With its impressive advantages of exact calculation and simple formulation due to the feature of recursive structure, many statistical or physical problems can seek an alternative but reliable solution from this modelling method, e.g. the travelling-salesman problem [3], the Ksatisfiability [4], the glass transition [5] and so on. For most cases, especially when involving complex structures, it is impossible to solve arbitrary models on a square lattice so one is forced to find approximated solutions. Usually, one attempts to solve the model in a mean-field approximation, while Gujrati [6] established that recursive lattice solutions are more reliable than the mean-field solutions, especially for the antiferromagnetic models.
On the other hand, the recursive feature of such lattices implies that the system is naturally homogeneous and suitable to describe bulky systems. Our group has already used the Bethe lattice to study inhomogeneous structures near a surface [7][8][9]. Here we extend to a one-/two-dimensional hybrid recursive lattice to study the phase transitions of Ising model on surface/interface, which is nonetheless still based on a symmetrical structure [10].
In this study, we have constructed an asymmetrical two-dimensional recursive lattice with a one-dimensional surface. By fractalizing the analogue of a diagonally cut regular square lattice, we can obtain a simple model with partial Husimi trees hung on a zigzag surface. Since this zigzag structure is infinite and homogeneous, it can be handled by recursive calculation [11,12], with approximating partial Husimi trees to be constant statistical contributions, which was derived from previous classical works, the exact calculation of this model appears to be feasible, like other Bethe-like models. An antiferromagnetic Ising model has been solved on the lattice, and we can locate the first-order phase transition and the Kauzmann paradox in the +1 spin system.

Lattice and model
Conventionally, the surface of a two-dimensional square lattice is just the one-dimensional line confining the bulk with a lower coordination number of 3. Imagine we cleave the lattice along square diagonals, a zigzag surface can be generated for the one-half bulk, and the basic unit of this lattice is square in the bulk and triangle on the surface. The triangle unit is surrounded by two identical triangles and a bulk square, thus a straight connection of triangle units can be adopted to represent the surface. This triangle chain is infinite and has the recursive feature required for calculation. Considering that bulk parts are merely structures hung on the zigzag line, on each triangle we may attach an infinite partial Husimi tree to represent the bulk part. The scheme of fractalizing the halved structure to be a recursive lattice is shown in figure 1.
Since this zigzag surface recursive lattice (ZSRL) has coordination of four inside the bulk and alternating two or four (averagely three) on the surface, we hope that it is a good approximation to the regular square lattice with a diagonal boundary. The site labelling is shown in figure 2: the vertices on surface triangle are labelled S and S 0 , where S 0 marks the site closer to an imaginary origin point on the surface, following the direction of recursive calculation. Note the S 0 in one triangle is S in the unit of lower level while the origin is marked as level 0. The base site that links the surface triangle and bulk tree is labelled S B .
Corresponding to the notation of neighbour interaction J, diagonal interaction J P and three spins interaction (triplet) J 0 , on the surface unit we label a bar on each term to differ them from bulk parameters. Then there are two neighbour interactions J, one diagonal interaction J P and one triplet ( J 0 ) in the Hamiltonian of one triangle unit a, which is where H is the magnetic field applied to each spin. Note the H of site S 0 is not included in the last term, because it will be counted in the unit of next level. In equation (2.1), a negative value of J will raise the lowest energy with different neighbouring spins states, the system is thereafter the antiferromagnetic case, or vice versa. The detailed effects of energy parameters will be discussed in §4.2.

Solutions on the surface
The general calculations of Ising model on Husimi square lattice can be found in previous works [6,13]. Here we follow the similar method of partial partition function (PPF) to achieve the solution x, which presents the probability of one site to be occupied by + spin on the surface site. The triangle unit has 2 3 ¼ 8 possible configurations, four of which are with the base site S 0 ¼ þ1 and the others are with royalsocietypublishing.org/journal/rsos R. Soc. open sci. 6: 181500 S 0 ¼ 21. We can derive two PPFs of the triangle on a level m from PPFs of the higher level m þ 1, the contribution from bulk Z B (S B ), and the local weight w(G), and where G is the index of configuration. Note that the total partition function (PF) of the entire system at the origin site S 0 is given by This relationship is fundamental when we derive thermal quantities, e.g. free energy from the PPF. We introduce ratios and as the solutions on surface at level m. By denoting royalsocietypublishing.org/journal/rsos R. Soc. open sci. 6: 181500 By introducing a set of polynomials and we can derive the solution x m at the mth level as a function of the ratio x mþ1 on higher level and the bulk solution x B , Taking x B as the solution on a joint site between surface triangle and bulk square, imagine we start the recursive calculation from a position deep in the bulk, then approach the thermodynamic contributions to the surface, thereby x B is simply the fix-point solution of Husimi lattice. In this way, we can count the contribution of the infinite bulk tree as a constant input and focus on the recursive approach of x along the surface. It should be addressed that this set-up of constant x B ignores the possible backward effect from surface to the bulk, which may bias the numerical value of x B . However, this approximation is acceptable to make the model simple and solvable.
With a constant x B , by equation (3.8) we can recursively calculate the solution on surface for a number of iterations until we reach a fix-point solution. The form of equation (3.8) implies that, regardless of the system being antiferro-or ferromagnetic, a uniform 1-cycle solution is expected on the surface, while antiferromagnetic Ising model presents an alternating 2-cycle solution as the ordered state [13]. A negative neighbour interaction J prefers to anti-align the S versus S B and S 0 versus S B pair. Unless we set the diagonal interaction J P also to be negative and large enough to outweigh J, the system will prefer the same spin states on S and S 0 .
Recall that bulk solutions can be either a 1-cycle solution to present the metastable state, or a set of 2cycle solutions as the ordered state [13]. Taking the 1-cycle x B , which is usually 0.5 with H ¼ 0, we will obtain a fixed x ¼ 0:5 also corresponding to a metastable surface; for the 2-cycle solutions, we can substitute either fixed solution as x B , and it will affect the calculation on x to bias the surface fix-point solution to 0 or 1. For example, due to antiferromagnetism, at T ¼ 0 with H ¼ 0, we will have 0 and 1 solutions in the bulk, then if we take x B ¼ 0, the surface solution will be captured as x ¼ 1, and vice versa. Our results confirmed that the thermodynamics calculated based on either selection are identical. royalsocietypublishing.org/journal/rsos R. Soc. open sci. 6: 181500 Figure 3 shows the ordered solution on the surface with comparison of 2-cycle bulk solutions, the energy terms are set as J ¼ 21, J ¼ À1 and others to be 0. The metastable 0.5 solution is not presented in figure 3. For convenience, we are still going to call the stable solution and corresponding thermodynamics '2-cycle' in the following discussion, although both stable and metastable x are actually in 1-cycle form for ZSRL.
It can be observed that at high T all the solutions are 0.5, that spins anywhere have an equal probability to be +1. Below the Curie point, the spins undergoes self-magnetization (even without an external field H), and the neighbouring spins prefer the þ/2 alternating arrangement, which is the lowest energy state. The value of 2-cycle x B of bulk Ising model is referred from previous reports [10,13]. Along with the bulk solutions, the spins on surface also present preference on unitary direction under the bulk Curie point; however, we will show that the actual phase transition occurs far below the T C (bulk).
One more thing should be addressed here. We know from exact solutions that for the square lattice the self-magnetization occurs at T C ¼ 2= log 1 þ ffiffi ffi 2 p 2:27. In the recursive lattice (RL), the bulk transition temperature is 2.8 as shown in figure 3, which is not all that close to the known T C , but closer than the results of mean-field theory. We believe this overestimate is acceptable. The aspect that the RL method provides results between the 'real' exact solution and the mean-field theory seems to be generic. Similar results had also been observed in other works, e.g. the three-dimensional cube lattice [13]. However, the nature of this overestimate has not been investigated yet.

Free energy calculation
We follow the Gujrati trick to calculate the free energy of a local area by recursive approach [6,13]. The scheme will be briefly described here, as shown in figure 4. Owing to the uniform structure and solution on the surface, we can randomly select a site as the origin point O. The local area is chosen to be two triangle units joint on the origin. Imagine we cut off two sub-trees contributing to the point A and A 0 then rejoin them together to make an identical but smaller ZSRL, and hook up two partial bulk trees hung on B and B 0 together to make a full Husimi square lattice; therefore we have F total ¼ F local þ  royalsocietypublishing.org/journal/rsos R. Soc. open sci. 6: 181500 With the definition of Helmholtz free energy F ¼ 2k B T log Z, the F per site in the local area is given as since there are three full sites in the local area (on whole site S 0 and four half-shared sites S 1 , S B ) and k B is normalized to be 1. Recalling equations (3.1)-(3.3), it is easy to break down the local free energy into a function of PPFs. Then from the relations between Z(+), Q and x derived in the previous section, we can calculate the above function as then the entropy and energy ( per site) can be easily achieved by S ¼ 2dF/dT and E ¼ F þ TS.

The thermodynamics and transitions on the surface
The thermal behaviours of the reference set-up with J ¼ 21, J ¼ À1 and other parameters to be 0 is shown in figure 5. The free energy of two solutions differ at T ¼ 2.8. Usually this bifurcation indicates the spontaneous magnetization (Curie point) in normal spin models. However, here on the surface the magnetization does not bend the free energy of alternating spins arrangement below the disordered 0.5 solution but upward, which is unphysical with a sharply increased entropy. Therefore, the 2-cycle solution at this point is only a numerical existence and the system must follow the curve of 1-cycle solution, until a cross point is reached at T ¼ 1.33, where the system makes a transition from the amorphous to the crystalline ordering, i.e. the order-disorder transition at critical temperature T C . Unlike the conventional self-magnetization where entropy is continuous, the entropy here must undergo a discontinuous jump like a first-order transition.
In this way, we may conclude that the critical order-disorder transition on the surface is not the Curie point, but much lower than that of the bulk. This can be understood as that, due to the asymmetry and  Below the real T C (bulk) and the apparent Curie point on the surface, the spins on surface can maintain their 'melted' state in a deep temperature region, and this does not fall into the concept of supercooled liquid, since the corresponding crystal state in particular temperature region is unfeasible. Intuitively, this phenomenon may easily refer to the analogue of ice premelting [14]; nonetheless, this work only focuses on the preliminary modelling and we are not going to further expand this point.
For the 1-cycle solution below T C , the system can undergo cooling without any phase transition and shares features of a supercooled state. With further cooling process, the entropy of 2-cycle solution approaches zero, while the entropy of 1-cycle becomes negative at T ¼ 0.69, which is the Kauzmann paradox at T K [13].
The result indicates that a free surface dramatically decreases the transition temperatures. Figure 6 shows the free energy comparison of Husimi bulk system and ZSRL. This observation agrees with others' work on phase transitions on the surface or thin film, for example, the glass transition of polymer system in confined geometry [7][8][9]15]. The fact that similar reduction can be observed in our monoatomic model implies that the lower transition temperature on surface/thin film basically originates from the dimension reduction and less interaction constraints.

The effects of secondary energy parameters
Besides mathematical curiosity, there is always a primary expectation to describe and study real systems for establishing a theoretical model. In this way, other than the interaction J and J between the nearest neighbours, further interactions are included to make the model more versatile to describe various systems, as shown in equation (2.1). With adjustable combinations of energy parameters set-up, we can manipulate the thermal behaviour of system and the transition temperatures to better match the reality in particular situations. In other words, more adjustable parameters may serve as useful tools for the theoretical modelling to be correlated with experimental parameters. By the control of J ¼ 21 and other parameters to be 0 inside the bulk, we explored the effects of J, J P and J 0 . The secondary parameters are expected to either comply or compete with J, and some interesting phase behaviours are found with particular set-ups.

The surface nearest-neighbour interaction J
Considering the feature of asymmetry on the boundary, we set the nearest-neighbour interaction on the surface, denoted as J, to be differentiated and adjustable from the J, to make the model capable to describe some particular situations, e.g. the surface tension. The transition temperatures with J ¼ 21, other parameters to be 0, and J ¼ À0:5, 20.7, 20.9, 21.1, 21.3 and 21.5 are shown in table 1. As negative J complies with antiferromagnetic set-up J ¼ 21, the larger absolute value of J makes the system more stable and increases both critical and ideal glass transition temperature; and the relative length of supercooled region, indicated by the ratio T C /T K , gradually increases. However, for the J ¼ À0:5 decreased from 20.7, both T C and T K fall while the latter changes more dramatically,  royalsocietypublishing.org/journal/rsos R. Soc. open sci. 6: 181500 raising a much larger supercooled region. This critical phenomenon observes that a loosened surface with too weak interactions is easier to be supercooled.