The model given by Eq. (1) is studied here on the square lattice and will be referred to as the pure BC model. However, our main focus, on the other hand, is the case with bond disorder given by the bimodal distri- bution so that r reflects the strength of the bond randomness. The resulting quenched disordered (random-bond) ver- sion of the Hamiltonian defined in Eq. (1) reads now as A first comparative study between the two versions (pure and random) of the 2d BC model has been already presented in our earlier paper [37] for two values of the crystal-field coupling corresponding to the second-order (Δ = 1) and first-order (Δ = 1 . 975) regimes of the pure model where for the random version the disorder strength r = 0 . 75 / 1 . 25 = 0 . 6 was chosen in both cases. In the next Sections our study on the 2d BC model is extended to several values of the crystal field and disorder (listed in the table of Sec. III). Details of our simulations are summarized in the next Section together with an up to date brief sketch of our entropic scheme. achieve an accurate estimation of all finite-size anoma- lies. Of course, taking the union of the corresponding subspaces, insures accuracy for the temperature region of all studied pseudocritical temperatures. The up to date version of our implementation uses a combination of several stages of the WL process. First, we carry out a starting (or preliminary) multi-range (multi-R) stage, in a very wide energy subspace. This preliminary stage is performed up to a certain level of the WL random walk. The WL refinement is G ( E ) → f ∗ G ( E ), where G ( E ) is the density of states (DOS) and we follow the usual modification factor adjustment f$_{j}$ +1 = √ f$_{j}$ and f 1 = e . The preliminary stage may con- sist of the levels : j = 1 , . . . , j = 18 and to improve accu- racy the process may be repeated several times. However, in repeating the preliminary process and in order to be efficient, we use only the levels j = 13 , . . . , 18 after the first attempt, using as starting DOS the one obtained in the first random walk at the level j = 12. From our expe- rience, this practice is almost equivalent of simulating the same number of independent WL random walks. Also in our recent studies we have found out that is much more efficient and accurate to loosen up the originally applied very strict flatness criteria [40, 41]. Thus, a variable flat- ness process starting at the first levels with a very loose flatness criteria and assuming at the level j = 18 the orig- inal strict flatness criteria is now days used. After the above described preliminary multi-R stage, in the wide energy subspace, one can proceed in a safe identification of the appropriate energy subspace using one or more al- ternatives outlined in Refs. [40, 41]. In random systems, where one needs to simulate many disorder realizations, it is also possible and advisable to avoid the identifica- tion of the appropriate energy subspace separately for each disorder realization by extrapolating from smaller lattices and/or by prediction from preliminary runs on small numbers of disorder realizations. In any case, the appropriate subspaces should be defined with sufficient tolerances. In our implementation we use such advance information to proceed in the next stages of the entropic sampling. The process continues in two further stages (two-stage process), using now mainly high iteration levels, where the modification factor is very close to unity and there is not any significant violation of the detailed balance condition during the WL process. These two stages are suitable for the accumulation of histogram data (for in- stance energy-magnetization histograms), which can be used for an accurate entropic calculation of non-thermal thermodynamic parameters, such us the order parame- ter and its susceptibility [41]. In the first (high-level) stage, we follow again a repeated several times (typically ∼ 5 − 10) multi-R WL approach, carried out now only in the restricted energy subspace. The WL levels may be now chosen as j = 18 , 19 , 20 and as an appropriate starting DOS for the corresponding starting level the av- erage DOS of the preliminary stage at the starting level may be used. Finally, the second (high-level) stage is ap- In the last few years we have used an entropic sam- pling implementation of the WL algorithm [38, 39] to study some simple [40, 41], but also some more complex systems [8, 42, 43]. One basic ingredient of this imple- mentation is a suitable restriction of the energy subspace for the implementation of the WL algorithm. This was originally termed as the critical minimum energy sub- space (CrMES) restriction [40, 41] and it can be carried out in many alternative ways, the simplest being that of observing the finite-size behavior of the tails of the energy probability density function (e-pdf) of the system [41]. Complications that may arise in random systems can be easily accounted for by various simple modifications that take into account possible oscillations in the e-pdf and expected sample-to-sample fluctuations of individ- ual disorder realizations. In our recent papers [8, 37, 43], we have presented details of various sophisticated routes for the identification of the appropriate energy subspace ( E 1 , E 2 ) for the entropic sampling of each random real- ization. In estimating the appropriate subspace from a chosen pseudocritical temperature one should be careful to account for the shift behavior of other important pseu- docritical temperatures and extend the subspace appro- priately from both low- and high-energy sides in order to B. An outline of our implementation of the Wang-Landau approach 3 P ( J$_{ij}$ ) = 1 2 [ δ ( J$_{ij}$ − J 1 ) + δ ( J$_{ij}$ − J 2 )] ; (2) J 1 + J 2 2 = 1 ; J 1 > J 2 > 0 ; r = J 2 J 1 , H = − ∑ <ij> J$_{ij}$s$_{i}$s$_{j}$ + Δ ∑ i s 2 $_{i}$. (3)