The Behavior of the Three-Dimensional Askin–Teller Model at the Mixed Phase Region by a New Monte Carlo Approach

The new approach of performing Monte Carlo (MC) simulations, which eliminates large oscillations of the values of the thermodynamic quantities computed for a mixed phase region, is demonstrated. The results are presented on the example of the mixed phase region in the 3D Askin–Teller (AT) model, where within a certain range of parameters with equal probabilities there appear two different, but equivalent, ways of ordering two of the three order parameters showing independent behavior. This new approach allowed us to exploit magnetization and internal energy curves, Binder cumulant, Challa- and the Lee-Kosterlitz-like cumulants as well as the internal energy distribution histogram. According to the most effective strategy, in the critical region, we use our recently proposed cluster MC algorithm and the Metropolis algorithm beyond it wherever it is applicable. The existence of two tricritical points and the bifurcation point in this area of the phase diagram is confirmed, and their locations are determined. It is explained that although the system as a whole does not show the presence of latent heat at the boundary of the mixed phase region and the antiferromagnetic phase, it does occur for various order parameters. Specifically, the increase in the energy of the degrees of freedom of one kind is accompanied by an equal decrease in the energy of the degrees of freedom of the other kind.


Introduction
For many decades, the Ashkin-Teller (AT) model [1] is one of the important points of reference in statistical physics, as it is a nontrivial generalization of the intensively used Ising model. After Fan [2] we express the AT model in terms of two standard Ising models put on the same lattice with spins s i and σ i residing at each lattice site i, i.e., s i and σ i are the variables that can take values +1 or −1. This means that we take into account only two spin interactions of constant magnitude J 2 between the nearest neighbors. These two independent Ising models are coupled by the four-spin interaction of a constant magnitude J 4 , also only between couples of nearest neighboring spins. Thus, the effective Hamiltonian H is of the form Fig. 1 The investigated phase diagram area of the symmetric 3D AT model on a cubic lattice. The dotted curves denote the first-order phase transitions, the solid ones stand for continuous transitions, and the dottedbroken one indicates transitions of both types for different order parameters. In the phase labeled "Baxter" the system is ferromagnetically ordered with all order parameters s , σ and sσ nonzero, whereas in the phase labeled "para" they all are zero. In the phase " sσ AF " there is s = σ = 0 and only the product sσ is antiferromagnetically ordered. In the phase " σ " called the mixed phase region sσ = 0 and either s or σ is ferromagnetically ordered, but the other is not. The positions of labeled points inside the phase diagram are marked by + and A, H, H', K, K' are the tricritical points In the phase diagram in Fig. 1 there are also the Baxter and paramagnetic (labeled as "para") phases for which all order parameters, s , σ , and sσ , are ferromagnetically ordered and are zero, respectively. For the phase labeled as " sσ AF " s = σ = 0 and only the product sσ is antiferromagnetically ordered. First-order phase transitions are denoted by the dotted curves, whereas the continuous ones are represented by the solid curves. The labeled point positions are marked by + and A, H, H', K, and K' are the tricritical points. Arnold and Zhang [22] using MC simulations obtained the first more precise results along the line AP which occurs for the ferromagnetic four-spin interaction, i.e., for positive K 4 parameter values, and is only partially seen in Fig. 1 in the lower right corner. More accurate results for this area of the phase diagram are the subject of our recent paper [14]. Ising phase transitions occur along the continuous curve ending at the tricritical point K [20,21,32]. The 2D AT model shows the line of continuously varying phase transitions at K 4 ≤ K 2 first shown in the paper [33] which is still of current interest [34]. The MC simulation results suggest the possibility of occurrence of nonuniversal behavior also in the 3D AT model [17,19,21,[29][30][31] but our recent results indicate only a wide crossover along the AH line and the rare coexistence of continuous and first-order phase transitions along the HH' line [17] shown in Fig. 1. It is noteworthy that the character of continuous phase transitions along the HK' line is still an open question [31].
The first aim of this paper is to present our new approach of performing a MC computer experiment to study the order and phase transitions between the mixed phase region " σ " and the "para" phase on the example of the standard 3D AT model shown in Fig. 1. The " σ " phase occurs only in the symmetric AT model in 3D. Due to the difficulty in obtaining unambiguous results and their interpretation, the " σ " region is often omitted, as it is the most complex and the least recognized region of the 3D AT model phase diagram. To solve this problem, we have developed an appropriate strategy, precise tools such as three types of cumulants, and the energy distribution histogram, which enable a detailed analysis of this region. We have recently used these tools with success to analyze the first-order phase transitions to the right of point A shown in Fig. 1 [14]. The idea behind our method was recently explained in [35] on the example of a phase transition at K 4 = −0.3. Due to the recently announced presence of metastable and unstable states [4], we use our Wolff type cluster algorithm [36] in the critical region wherever it is applicable (i.e., for |K 4 | < K 2 ) and the Metropolis one beyond.
Using our method, we have examined the phase transitions along the HK'K b line between the mixed phase region " σ " and the paramagnetic phase labeled "para" in Fig. 1. Moreover, we have investigated the phase transitions along the KK b line between the antiferromagnetic phase " sσ AF " and the "para" one, as well as along the K b E line between the " sσ AF " phase and the " σ " region. This is the fragment of the phase diagram that is difficult to analyze and where there were only a few preliminary bibliography results [21,29]. These problems are important and topical, and have not yet been solved in the bibliography of the subject.

The Applied Method
To study the mixed phase region of the subject, we exploit the MC computer experiment with importance sampling of states and consider the finite-size cubic samples of the symmetric spin lattice AT model, the behavior of which is fully determined by the Hamiltonian (1). These samples of size L 3 with periodic boundary conditions are large enough to be able to compute the thermodynamic limit of our results. We perform our computer experiments to predict the equilibrium behavior of the 3D AT model according to the statistical mechanics methodology.
The detailed description of the MC computer experiment constructed by us for the 3D AT model based on the analysis of the dependencies of thermal averages of such thermodynamic variables as magnetization or internal energy, and three different cumulants on the coupling constant, and on the dependencies of the internal energy distribution histogram on energy can be found in our recent papers [14,37]. In this Section, apart from the computational aspects, we present only its key elements and we focus primarily on explaining how to realize it in the mixed phase region " σ ".
First, we bring our system to the state of thermodynamic equilibrium using the appropriate number of MC steps that we have analyzed in our paper [36]. Moreover, in our MC computer experiments, in contrast to simple MC simulations, we not only compute thermodynamic quantities but also carefully determine their error bars. For this purpose, one program run consists of the computation of 6 to 24 partial averages, each independently calculated from approximately 10 7 MC steps. However, only every kth step contributes to the thermodynamic calculations (6 < k < 10), which is enough to avoid correlations between sampled configurations of our system using the Metropolis algorithm [36]. The problem of these correlations is radically smaller in the case of the cluster algorithm [38], which is also the case in our version of this algorithm [36].
Of course, we get a true picture of the phase transition only in the thermodynamic limit. To obtain reliable extrapolations of our results to the thermodynamic limit, we perform computations in our system also with the largest possible size L, which take many weeks at sequential processing. In order to get the results in a reasonable time using the MPI library, we have parallelized the processing in our computer experiments, obtaining almost perfect speedup on symmetric multicomputers [39].
The presence of metastable and unstable states in the 3D AT model was recently signalized using the mean field method [4]. It is well known that the mean field theory is a solid tool, especially suitable for the first view of the problem, and it does not provide quantitative consistence with the precise results. Nevertheless, it gives good qualitative insight into the problem (see, e.g., [40]). So, we generate equilibrium configurations of finite-size cubic spin samples for fixed values of our model parameters described above in the Hamiltonian (1) using our recently constructed cluster algorithm of the Wolff type [36] in the critical region wherever it is applicable (i.e., for |K 4 | < K 2 , see also [41,42]) and the Metropolis one beyond. This is the best strategy, also to optimize the time to obtain results with comparable error bars.
To prelocate a temperature-driven phase transition point, we fix a particular value of K 4 coupling and analyze Binder cumulant Q α,L ( , [17,20,32,43]), where M n α L denotes the n-th power of the order parameter α component, with α = s, σ or their product sσ , which are averaged over an ensemble of independent samples of the size L 3 . The lack of characteristic minima in the course of the Q α,L (K 2 ) dependences indicates that the phase transition can be continuous [19,21,43].
To check if there occurs the latent heat during a phase transition, that is, to unambiguously determine the character of a phase transition and to more accurately determine the location of this transition point, we also compute the Challa [44] and the Lee-Kosterlitz [45] like cumulants. Here E n α L is the n-th moment of the entire Hamiltonian (α = H ) or the interaction energy of α-degrees of freedom (α = s, σ , or their product sσ ) separately, which is averaged over an ensemble of independent samples of the size L 3 . Thus, we can compute the latent heat l α coming also from each kind of degrees of freedom α with α = s, σ , or their product sσ separately [17,18,21].
We analyze the dependences V α,L (K 2 ) which show characteristic local minima V min α,L [44] and U α,L (K 2 ) with characteristic local maxima U max α,L [45] at a fixed value of K 4 coupling in the close critical region. When the thermodynamic limit V min α,L value with its error bar is different from 2/3 and the U max α,L value with its error bar is different from 1, we conclude that a phase transition is qualified to be of the first order; otherwise we assume that the phase transition is continuous [14,17,21,44,45]. The latent heat l α coming from the whole Hamiltonian (α = H ) or from the interaction energy E α of α degrees of freedom with α = s, σ , or their product sσ , in the thermodynamic limit where E α,± = E α (K 2 → K 2,c | ± ), are determined on the basis of the Lee-Kosterlitz formula [45,46] V min α,L→∞ = and using the method proposed in [21]. Equation (5) was also independently obtained by Borgs et al. [47] from a more rigorous point of view. K 2,c in Eq. (4) is the critical value of K 2 coupling with the fixed value of K 4 . The quantity A V in Eq. (5) stands for L independent expression of the complicated form [45]. Analogously, one can determine the latent heat l α on the basis of the U α,L cumulant maxima values scaled to the thermodynamic limit using the Lee-Kosterlitz formula [45] U max α,L→∞ = where α = s, σ , sσ , or H , and A U stands for L independent complicated expression. We conclude from Eqs. (5) and (6) that values and locations of cumulant V α,L minima and of cumulant U α,L maxima scale linearly versus L −3 . Using Eqs. (5) and (6) and the method proposed in [21] one can calculate the values of E α,+ and E α,− and estimate the latent heat from Eq. (4). Moreover, the thermodynamic limit K min 2,α values of minima and K max 2,α values of maxima are far better estimates of the critical K 2 values than the preliminary ones obtained on the basis of the Binder cumulant Q α,L (K 2 ) dependences mentioned above.
This method [14] gives good results when the system has unambiguously determined equilibrium configurations of finite-size cubic spin samples for fixed values of our model parameters described above in the Hamiltonian (1).
However, in the mixed phase region, marked as " σ " in Fig. 1, the situation is more complicated: sσ = 0 and there are two equally probable phases in which either s or σ is ferromagnetically ordered, but the other is not in the thermodynamic limit. This causes the values of the computed thermodynamic quantities to oscillate [21,30] as during sufficiently long simulations both phases will appear with approximately equal probability. Therefore, to obtain clear results, we have developed a new approach, the idea of which was presented in our recent paper [35] on the example of the phase transition at K 4 = −0.3. We propose a conventional division of our system into two sublattices: the ordered one which will be marked with a capital letter and the unordered one marked with a capital letter S. Importantly, the decision to allocate the real spins σ and s to these conventional sublattices is decided only after each MCS is performed and the results from the spins with greater magnetization are systematically added to the results of the conventional sublattice while the results from the second spins are added to the results of the conventional sublattice S. The number of MCS must be large enough to compensate for the separation of our system into these two artificial sublattices and S in the paramagnetic region. Computations for the product of spins sσ do not need to be changed.
Thus, bearing in mind that in the mixed phase region " σ ", two phases are equally probable in which either s or σ is ferromagnetically ordered but the other is not, thanks to the conventional division into two sublattices, the first one with greater magnetization of spins, and the second one with smaller magnetization of S spins, for the analysis we get smooth dependences of thermodynamic quantities and of cumulant values on the coupling K 2 with a fixed value of the coupling K 4 .
We have located phase transition points precisely enough to be able to use another independent method to check for the presence of the latent heat during a phase transition with greater accuracy. We compute the probability P α,L of the internal energy E α,L appearance in the samples of finite size L 3 . As in the case of cumulants, the P α,L (E α,L ) values are computed independently for each degree of freedom α = S, , or their product sσ , and also for the whole Hamiltonian (1) denoted by α = H , at a critical value K 2,c . A characteristic histogram of the distribution of this energy E α with two peaks in the close critical region for first-order phase transitions can be observed [14,45,48,49]. The maxima of these peaks appear at the energy value E α,−,L for the ordered state and at E α,+,L for the unordered one for the samples of finite-size L. In contrast, for continuous phase transitions, only a single peak of the probability dependence P α,L (E α,L ) appears in the thermodynamic limit. This is an important clue for determining the character of a phase transition.

Results
We demonstrate our new way of performing MC computer experiments at the mixed phase region " σ " boundaries HK'K b with paramagnetic and K b E with antiferromagnetic phases shown in Fig. 1. In the " σ " region with equal probabilities there appear two different, but equivalent, ways of ordering two of the three order parameters: sσ = 0 and either s or σ is ferromagnetically ordered but the other is not. Moreover, we have investigated the phase transitions along the adjacent KK b line between the antiferromagnetic phase " sσ AF " and the paramagnetic one.
For the point with K 4 = −0.358 on two boundaries between the mixed phase region " σ " and the paramagnetic phase "para" as well as between the phase "para" and the antiferromagnetic phase " sσ AF " shown in Fig. 1 we observe that the values of the cumulant V α,L (K 2 ) minima and the cumulant U α,L (K 2 ) maxima for α = S, , sσ , or H in the thermodynamic limit scale to 2/3 and 1, respectively, in terms of their error bars. This means that the phase transitions for all three order parameters , S , and sσ are continuous. To understand the situation, in Fig. 2 one can see the magnetization per site M 2 α L value changes for all three order parameters, , S and sσ . The latter two show no order in the paramagnetic phase, nor in the mixed phase region, as well as and S show no order in the " sσ AF " phase, in the thermodynamic limit. However, nonzero magnetization values result in the appearance of cumulant V minima and cumulant U maxima for all three order parameters. Therefore, in the lower graph of Fig. 2 for the order parameter S we observe a decrease in the magnetization value per node with an increase in the system size L, as opposed to the increase in magnetization for the order parameter labeled with α = in the " σ " region and for the order parameter sσ labeled with α = sσ AT computed from the one of two sublattices in the antiferromagnetic phase " sσ AF ". A similar decrease in the magnetization value per node for the order parameter sσ with an increase in L can be observed for the curves labeled with α = sσ in Fig. 2 computed for the entire lattice which indicates that there is the antiferromagnetic order in the phase " sσ AF ". Moreover, for large samples with L ≥ 72 the dependences V S,L (K 2 ) minima and the dependences U S,L (K 2 ) maxima disappear, as expected [14,21,35,44,45].
The abscissa K min 2,α,∞ of the minima positions of the dependences V α,L (K 2 ) and the abscissa K max 2,α,∞ of the maxima positions of the dependences U α,L (K 2 ) for α = , or sσ as well as the respective extremes computed using the energy of the whole Hamiltonian with α = H are consistent in the thermodynamic limit in terms of their error bars for each phase transition. The abscissa of these extremes for phase transition on the upper line HK' scale to K 2,c = 0.35983 (12) and on the lower line DK scale to K 2,c = 0.35716 (16) which are the average values calculated from the respective critical values of K min 2,α,∞ and K max 2,α,∞ at a fixed value K 4 = −0.358. The method of determining the error bars of the quantities extrapolated to the thermodynamic limit is explained below in Figs. 4 and 5. These are the more precise values than the preliminary estimates of K 2,c = 0.3605(8) and 0.3566(7) obtained from the analysis of the cumulant dependences Q α,L (K 2 ) (for the details, see [21,32,43]). We observe a similar behavior along the entire phase transition HK' line (see also [35]) and along the line DK. We have divided our further analyzes into three subsections. In the first we examine the behavior of our system in the vicinity of the tricritical points K and K', in the second subsection we focus on the bifurcation point K b region, and in the third we analyze the behavior of the system along the K b E line, which are shown in Fig. 1.

The Vicinity of Tricritical Points K and K'
With the value of the coupling constant K 4 = −0.366 (see Fig. 3) for L = 28, we observe two clear minima of the dependence of the cumulant V H ,L on the coupling constant K 2 . The index α = H represents the results obtained using the energy of the entire Hamiltonian and each minimum is the result of a different phase transition. The abscissa of the positions of these minima in the thermodynamic limit given in the third and fifth columns of Table 1 scale to different critical values of K 2 taking into account their error bars. The minimum on the left hand side for α = H in Fig. 3 corresponds to the transition between the antiferromagnetic and paramagnetic phases denoted respectively as " sσ AF " and "para" in Fig. 1. The abscissa of these minima in the thermodynamic limit scale to the critical value of K min,l 2,H ,∞ = 0.36592 (14) given in the fifth column of Table 1. The same phase transition corresponds to the upper left minimum in Fig. 3 for α = sσ and the abscissa of these minima in the thermodynamic limit scale to the critical value K min,l 2,sσ,∞ = 0.36582 (12) provided in the fourth column of Table 1. These two values agree within the limits of their error bars. On the other hand, the abscissa of the minima on the right-hand side in Fig. 3 scaled to the thermodynamic limit give respectively the critical values K min,u 2,H ,∞ = 0.36639 (14) and K min,u 2, ,∞ = 0.36645 (12) given in the third and second columns of Table 1. These K 2 values correspond to the phase transition between the mixed phase region marked as " σ " and the paramagnetic phase marked as " para" in Fig. 1. These two values also agree within the limits of their error bars.
We have determined the presence and amount of the latent heat during both phase transitions by rescaling according to Eq. (5) to the thermodynamic limit of the values of cumulant V α,L minima with α = sσ , and H for different system sizes L. The results of these ana-   Table 2. Similar conclusions can be drawn from the results of analyzes for K 4 = −0.362, one row higher.
The situation is different for phase transitions with K 4 = −0.358 and K 4 = −0.36, for which the results are presented in the first two rows of Table 2. Here, all the values of the cumulant V α,L minima rescaled to the thermodynamic limit, taking into account their error bars, are not different from the limit value 2/3, which have been marked in bold. Therefore, these are continuous phase transitions.
The results of these analyzes allowed us to determine the position of the tricritical point K at (−0.360(1), 0.35951(11)) on the phase diagram in Fig. 1, i.e., with greater precision than in the paper [21]. Moreover, thanks to the proposed method of separating the system into two sublattices, depending on the magnetization value, we have determined the position of the tricritical point K' at (−0.360(1), 0.36110(12)) on the phase diagram for the first time in the bibliography. In earlier papers, the possibility of the existence of the tricritical point K' was signaled, but its location was not determined due to large oscillations of the values of the thermodynamic quantities computed for a mixed phase region using MC experiments [21].
Analysis of the behavior of analogous cumulant U maxima with the use of the scaling relation Eq. (6) has led to the same conclusions. Figure 5 presents abscissas K min 2,α,L of minima of cumulants V α,L with α = sσ , and H explained in the legend which are extrapolated to the thermodynamic limit at the fixed value of the coupling K 4 = −0.37 (the lower graph) and at K 4 = −0.376 (the upper graph). The legend also indicates the areas of the phase diagram (antiferromagnetic phase " sσ AF " or the mixed phase region " σ ") on which these results were obtained. The dependences are fitted by straight solid lines using linear regression, according to the scaling relation (5). We have marked the scaled position values with the symbols ×, which constitute the critical value of the coupling constant K 2,c for individual phase transitions together with their error bars.

The Bifurcation Point K b Region
For K 4 = −0.37 (the lower graph in Fig. 5) we observe two critical values K min,l 2,sσ ∞ = 0.36992 (12) and K min,l 2,H ,∞ = 0.36982 (15) consistent within the error bars on the lower line KK b and two matching values K min,u 2, ∞ = 0.37031 (14) and K min,u 2,H ,∞ = 0.37043(18) on the upper line K'K b . However, the values for these two pairs are significantly different, taking into account their error bars. Thus, our results confirm that here the phase transitions for the order parameters sσ and still occur along different lines in the phase diagram in Fig. 1. One can observe a different situation in the upper graph in Fig. 5 for K 4 = −0.376. Here, the four critical values K min,l 2,sσ ∞ = 0.376006(9), K min,l 2,H ,∞ = 0.376006 (29), K min,u 2, ∞ = 0.375986 (26), and K min,u 2,H ,∞ = 0.375998(60) resulting from the extrapolations are consistent within their error bars. Thus, the phase transitions for the order parameters sσ and occur along the same line K b E in the phase diagram in Fig. 1.
The selected values of our analyzes for this part of the phase diagram are presented in rows 4 to 9 of Table 1 (12)). For the ordinate, we have assumed here the mean value of all four values of the K 2 coupling constant given in the seventh row of Table 1. Our result is in line with the preliminary results from the previous papers [21,29] determining the position of the K b E line satisfying the asymptotic relationship K 2 = −K 4 .
Analysis of the behavior of analogous U cumulant maxima with the use of the scaling relation Eq. (6) leads to the same conclusions.

The K b E Line
Along the line K b E, shown in the phase diagram in Fig. 1, i.e., for the coupling constant value K 4 ≤ −0.371, the parameter describing four-spin interactions, rescaled cumulant V min ,∞ and V min sσ,∞ minima in the thermodynamic limit, taking into account their error bars, are significantly different from the limit value 2/3. This means that the phase transitions for the order parameters sσ and along the K b E line are of first order. Important information is provided here by the values of the cumulant minima V min,σ H ,∞ and V min,sσ H ,∞ , which in the thermodynamic limit within the error bars are in line with 2/3. These values are shown in bold on the last three lines of Table 2. Thus, the system as a whole shows no latent heat in the transition from the antiferromagnetic phase " sσ AF " to the mixed phase region " σ ". This is due to the fact that when going from the " sσ AF " phase to the " σ " region, the product sσ undergoes the phase transition from order to disorder, increasing the associated internal energy, while for the same transition, the internal energy related to the degrees of freedom is reduced by the same amount, because it is a transition from disorder to order. Thus, the internal energy of the system as a whole does not change within its error bar.
These energy changes are plotted on the left-hand side of Fig. 6 at K 4 = −0.376 and for L = 40. As we have explained in the last paragraph of Sect. 2, at the critical point, one can compute the latent heat value with good precision using the energy distribution histogram method. The energy distribution histograms for the system size L = 40 for the critical point K 2,c = 0.375999(31) (the average value of the critical values K 2,c provided in the last row of Table 1) and for the fixed value K 4 = −0.376 are shown in the graphs on the right hand side of Fig. 6.  Fig. 6 The energy dependences of the E α,L distribution histogram P α,L (E α,L ) shown on the graphs on the right-hand side for the 3D AT model for degrees of freedom specified in the legend at K 4 = −0.376 and for L = 40. The left-hand side graphs present internal energy dependences E α,L (K 2 ) with α = H , S, sσ , also at K 4 = −0.376 and for L = 40. The critical coupling value K 2,c = 0.375999(31) is marked with a vertical dashed line. Labels " sσ AF " and " σ " inform about the type of order in the system in a given interval of K 2 coupling values. For clarity only selected points are shown, and error bars are omitted, as they are less than magnitudes of symbols Using this method, one can more accurately determine the latent heat l α as the change of internal energy l α = E α,+,L − E α,−,L for individual degrees of freedom α = , sσ , as well as for the entire system α = H . From our analyzes, we have l sσ = 0.090(10), l = 0.086 (5) and both of these values agree within their error bars. Due to the opposite direction of change of the internal energy of the product sσ and of degrees of freedom during this phase transition, a decrease in the value of the former is accompanied by an increase in the value of the latter, and the total energy change of the system as a whole is 0.004 (15), which within the error bar is consistent with 0. Therefore, we have observed no latent heat not only for the S degrees of freedom but also for the system as a whole l H , as one can see in Fig. 6, where for α = S and H the single probability peaks P α are observed. The presence of single probability peaks of the energy occurrence indicates the absence of latent heat during a phase transition as mentioned in the last paragraph of Sect. 2.
The results obtained by the method of the internal energy distribution histogram are consistent within their error bars with the results determined with the use of the cumulant V and U properties.

Conclusions
We conclude that our MC computer experiment and the method explained in our recent paper [14] have been successfully extended to investigate the mixed phase region. The results of our computer experiments clearly confirm that in the mixed phase region " σ " of the 3D AT model there are two equally probable phases in which either s or σ is ferromagnetically ordered but the other is not in the thermodynamic limit. As a consequence, the values of the computed thermodynamic quantities oscillate when we apply our method presented in [14] in its original form, since during sufficiently long simulations both phases appear with approximately equal probability [21,30]. These oscillations of the values of the computed thermodynamic quantities make their interpretation very difficult. Therefore, to obtain clear results, we have demonstrated and carefully checked our new approach, the idea of which was presented in our recent paper [35] on the example of the phase transition at K 4 = −0.3. We have used a conventional division of our system into two sublattices: the ordered one and the unordered one and the decision to allocate the real spins σ and s to these conventional sublattices and S is decided only after each MCS is performed and the results from the spins with greater magnetization are systematically added to the results of the conventional sublattice while the results from the second spins are added to the results of conventional sublattice S.
In contrast to the results published so far in the bibliography, the implementation of our new approach in the MC experiment allowed us to use Binder cumulant, the Challa-like and the Lee-Kosterlitz-like cumulants, as well as the internal energy distribution histogram to obtain clear results for a mixed phase region. Along the line HK'K b E (see Fig. 1), according to the most effective strategy, in the critical region we have used our recently proposed cluster MC algorithm [36] and the Metropolis algorithm beyond the critical region, which are suitable for both the first-order and continuous phase transitions in the 3D AT model. The KK b line is outside the applicability range of our cluster algorithm of the Wolff-type |K 4 | < K 2 [36,41,42].
We have demonstrated our method on the mixed phase region " σ " boundary where there are only a few preliminary bibliography results [21,29], i.e., along the line HK'K b E. Moreover, to complete our study of this area of the 3D AT model phase diagram, we have examined the behavior of the system along the adjacent line KK b . It is worth mentioning here that the AT model is an important reference point in statistical physics. However, this method can be successfully applied to other spin lattice models whose phase diagram contains a mixed phase region.
Since for the conventional degrees of freedom and S we have gathered contributions from both equally probable phases: this with σ nonzero and s = 0 as well as with σ = 0 and s nonzero, thus we have shown that the phase transitions for K 4 ≤ −0.36 are continuous. These continuous phase transitions end at the tricritical points K with the (K 4 , K 2 ) coordinates (−0.360(1), 0.35951 (11)) and K' with (−0.360(1), 0.36110 (12)). These values are in line with those obtained earlier [21] and are more accurate. We have unequivocally shown that for K 4 < 0.36 the first-order phase transitions take place.
The preliminary results of Ditzian et al. [29] indicate the existence of the single K point from which the boundary between the antiferromagnetic phase sσ AF and the mixed phase region σ begins. In turn, the results of Musiał et al. [21] already have indicated the existence of two different tricritical points K and K' with the lines of phase transitions KE and K'E asymptotically approaching each other on both sides of the relation K 2 = K 4 . However, the results of this paper clearly indicate the existence of the bifurcation point K b , and its location is (−0.371(1), 0.37103 (12)). For K 4 < −0.371, the " sσ AF " phase and the " σ " region are separated by the line K b E.
The results presented show that to explain the complex behavior of the 3D AT model, it is not enough to compute the latent heat of the entire system, but to fully understand the behavior of the system, it is necessary to consider individual contributions to the latent heat from different degrees of freedom, which was originally proposed in the paper [21]. We have discovered and explained that although the system as a whole does not show the presence of the latent heat at the boundary F b E between the mixed phase region " σ " and the antiferromagnetic phase " sσ AF ", this heat does occur for various order parameters. Specifically, the increase in the energy of the degrees of freedom is accompanied by an equal decrease in energy for the product sσ .