DEVELOPING A TEST FOR RANDOM NUMBER GENERATORS USING A SIMULATION OF THE HIERARCHICAL POTTS DIAMOND MODEL AT THE CRITICAL POINT

This dissertation uses the hierarchical q-state Potts model at the critical point to develop a new random number generator test. We start with an exposition of renormalization group approach by means of which one can numerically exactly compute the free energy, speci c heat and susceptibility of large, but nite lattices. We then show that generalization of these standard techniques allows one to also compute probability distributions related to the energy and the order parameter. The various computed quantities can be compared with Monte Carlo estimates of the same quantities. We demonstrate that the structure of the hierarchical lattices used allows one to perform the Monte Carlo calculations by direct sampling. This avoids the usual critical slowing down that plagues Monte Carlo calculations at the critical point. As is well known, critical behavior is highly susceptible to perturbations. We expect that aws of the pseudo random number generator, such as correlations, will cause statistically signi cant discrepancies between the results of the simulations and the numerically exactly computed results. Details of the computer code generated for these tests are included.

The energy (pair-coupling) and heat capacity graphs per bond for hierarchal systems L 2 , L 3 , L 4 , L 5 , L 7 , and L 9 with κ c in- Carlo algorithms include optimization, integration, and making draws from a probability distribution. Results obtained by Monte Carlo methods can be compared with results obtained by numerically exactly calculating thermodynamic properties of statistical mechanical model systems dened on hierarchical lattices. We use this to design a new pseudo random number generator test. We shall henceforth use the acronym RNG in which pseudo is implicitly assumed.
In a classical study [1] titled Monte Carlo simulations: Hidden errors from 'good' random number generators, Ferrenburg demonstrated that random number generators deemed reliable produced systematic errors in Monte Carlo calculations of physical quantities when the generators were used to drive the fast Wol spin-cluster-ipping algorithm. That the results were biased was determined by comparison with the exact solution to the two-dimensional Ising model by Lars Onsager, which in later work by Ferdinand and Fisher [2] was applied to nite lattices. This dissertation proceeds in the same spirit, and adds some new features to the use of exactly solvable, statistical mechanical models as testing grounds for RNGs.
This dissertation has the following goals: 1. Application of a numerically exact method to a nite hierarchical q-state Potts diamond lattice to obtain thermodynamic properties such as specic heat and susceptibility as well as various probability distributions; 2. Design a direct sampling Monte Carlo method that produces statistical estimates of these quantities; 3. Combine the two approaches to verify the validity of calculated estimators of expectation values and histograms of probability distributions by comparison with dedecoration results. These comparisons are used to verify how well a random number generator performs.
We expect that the sensitivity to correlations of critical systems makes the proposed method a good tool for RNG testing and a useful addition to the standard test suites that are used for this purpose [3], [4].

The Layout of the Dissertation
Chapter 2 reviews the application of dedecoration to two specic lattice types.
These two lattices are the linear periodic chain lattice and the hierarchical diamond lattice. It is shown that thermodynamic quantities such as the free energy are constructed as series expressions that can be obtained by means of recursion in nite lattices or approximated to a desired degree of accuracy in innite lattices. In the thermodynamic limit, the process of dedecoration makes contact with renomalization group theory. The critical exponents can in principle be obtained by scale invariance used with the models discussed. This shows why various quantities diverge in the thermodynamic limit, but for our tests only nite quantities computed for nite lattices are used.

Review of Renormalization Group Theory
This dissertation restricts itself to two types of lattices; the linear chain lattice and the hierarchical diamond lattice. Both will be specied by notation L l where l = 1, 2, ... denotes the recursion level of lattice L. In the chain lattice conguration L 1 denotes a system with two sites connected by two bonds. The lattices L l+1 of the chain are recursively constructed by placing a new site between existing sites; see Figures 1,2,3. In the hierarchical lattice, L 1 denes a system of two sites connected by one bond. Hierarchical lattices L l+1 are recursively constructed by replacing this unit by a diamond, as illustrated Figures 4,5,6.
We dene n s l and n b l as the number of sites and bonds of lattice L l . The relationship of sites and bonds for both lattices based on level l are determined by recursion relations. In the chain lattice, Note, in the chain lattice, n s l = n b l , which corresponds to a chain with periodic boundary conditions. In the hierarchical lattice, n s l = (4 l + 8)/6, The microscopic variables located on the sites of the lattice L l are denoted in general by congurations S l = (s 1 , s 2 , . . . , s n s l ). The site variables take on discrete values called states: ±1 for the Ising model and 1, . . . , q for the q-state Potts model. The process of adding a site between two existing ones is called decoration.
The construction of the diamond lattice requires bond doubling combined with decoration.

Recursion by Dedecoration
We associate with each pair of sites (s i , s j ) the energy (s i , s j ; K). K is the set of all interaction parameters, which will be specied in the following chapters.
We also adopt the convention that denes the zero of the energy and absorb the Boltzmann factor, −1/k B T , into the energy. This reduced energy will simply be called the energy. We introduce a transfer matrix corresponding to the fact that we consider the energy of the system to be a sum of only single-site and nearestneighbor contributions: Dedecoration is the process of summing over all sites with two nearest neighbors in lattice L l . In the linear conguration, since connectivity does not vary, this process is nothing but summing over the microscopic variable at every other site in the chain. Here, and in following chapters, we deal mostly with nite systems; we will only consider the thermodynamic limit for theoretical purposes when needed.
The rst lattice we discuss is the linear chain. The process of dedecoration is illustrated in Figure 7. In terms of matrix multiplication, dedecoration takes the following form: We see that dedecoration in this system reduces the lattice from three to two sites and introduces a changed interaction parameter K → K (K) along with a shifted term g = g(K) to reimpose the zero energy convention used. Note, the set of parameters K must be suciently general enough to make sure that no new ones are generated in the process of dedecoration. For the Ising chain, nearest neighbor coupling satises this requirement both in the presence and absence of a magnetic eld. The case of the q-state Potts model is more complicated as we shall see in  The essential transformation required to generate the recursion relations of the hierarchical lattice is illustrated in Figure 8. All that needs to be added, compared to the case of the linear chain, is the bond doubling: T (s 1 , s 4 ; K)T (s 4 , s 2 ; K) = g T (s 1 , s 2 ; K ). In operations like Eq. (6) one can immediately see the relationship to single site dedecoration, g (K) = g(K) 2 and K (K) = 2K (K). The change in interactions is visually shown in Figure 8 where summing over s 3 and s 4 reduces the hierarchical L 2 lattice to the hierarchical L 1 lattice with two modied parameters of interaction K . For the rest of this chapter we consider hierarchical lattices only; the modications between the two lattices are straightforward and this will allow the notation to remain simple.
The general Hamiltonian that describes all systems specied by conguration S l and the set of all interaction parameters, K, is written as H(S l , K). We remind the reader that −1/k B T has been absorbed into the energy. In terms of the transfer matrix of Eq. (4), the Boltzmann weight of a conguration of a lattice of level l is given by where i, j ∈ L l denote all nearest neighbor pairs in L l . We now dene the partition function of lattice level l using Eq. (7), Consider the hierarchical lattice L 3 and apply dedecoration once. To accomplish this transformation from L 3 to L 2 in Figures 5 and 6 we open up S l and act out all sums on sites s 5 , s 6 , . . . , s 12 with operations like Eq. (6). After dedecoration, the remaining sums act on the lowered lattice of modied interactions. Therefore, dedecoration relates a starting partition function of level l to partition function of level l − 1 with modied interaction parameters, The factor G l−1 (K) is a by-product of normalizing the shift in energy that comes from operations Eq. (6). Each bond after dedecoration contributes to the overall shift in energy and so G l−1 (K) has a power of bond contributions of the dedecorated lattice, G l−1 (K) = (g (K)) n b l−1 . Using dedecoration, we generate a recursive method of evaluating the partition function. Note, dedecoration can be repeatedly applied until our system reaches the point at which one chooses to calculate the partition function by explicit summation.

The Free Energy and Renormalization
Taking the logarithm of the recursion relation for the partition function in Eq. (9) gives a recursive relationship of the reduced free energy where f l (K) = F l (K)/n b l , and f l−1 (K ) = F l−1 (K )/n b l−1 . The recursive form of the free energy can be used to compute numerically exact results for nite lattices and approximations correct to any desired accuracy for innite models; see Chapters 3, 4 for applications.
In the thermodynamic limit Eq. (11) shows that the hierarchical model is self-similar if the rst term on the right is interpreted as self-energy associated with each site. We write f for the free energy per site in the thermodynamic In Eq. (12) all double primes have been dropped; we settle on a unique notation of g, K for the remainder of the chapter.
Renormalization group (RG) theory was developed to obtain critical point exponents describing critical point singularities and explain the observed universality of critical behavior [2] of the free energy per bond f (K) from the regular functions g(K) and K (K). In general there is no guarantee that g(K) and K (K) are regular at the critical point. However, in the case of the hierarchical lattices we consider, this property is rigorously satised by the summation over sites with only two nearest neighbors at every level of dedecoration; the dedecoration operation is a modied RG transformation. It is well-known that the scaling transformation implies that the free energy has power-law singularities that are characterized by critical exponents. We therefore continue to review the standard approach.
The renormalization group equations are Fixed points are points in parameter space that are invariant under the RG transformation: There typically are several xed points for any transformation, but some, such as zero and innite temperature xed points, are only indirectly relevant for critical behavior. We are interested in non-trivial xed points for our analysis. The linearized form of the transformation of interaction to interaction K = K (K) is determined by the Jacobian matrix Here, α, β are matrix indexes of the Jacobian. The xed point of the linearized transformation is With the Jacobian, the linearized transformation reads To simplify the transformation in Eq. (17) one can introduce normal mode coor- Thus, we obtain In normal-mode coordinates the RG transformation takes the form The interaction parameters are expressed in terms of the scaling coordinates (or scaling elds) since the coordinates are functions of the interaction parameters, . Note that λ = 1 will change nothing; this value of λ is neither a relevant or irrelevant change. λ = 1 is the marginal value and it sets the bounds for the two cases λ > 1 (relevant) and λ < 1 (irrelevant). We write and we shall treat b as a continuous variable. The justication for this is that the RG transformation can be iterated, but this approach obscures some subtleties that are of no interest here. We refer to to the literature for details [2].

Scaling Theory
We next review the implications of the scaling equation for the free energy We restrict ourselves to the case of two scaling elds, u 1 , a thermal eld and u 2 a eld coupling to the order parameter.

The Regular and Singular Free Energy
We consider two parts to the free energy: the regular part and the singular part. To keep things simple we look at the case of a single, thermal scaling eld, only, As discussed above g = g(u) is a regular function even at the critical point. If the regular part of the free energy satises we nd that the singular part satises the homogeneous equation For details, see the work of Niemeijer and van Leeuwen [2].

Widom Scaling
Eq. (23) was obtained for one scaling parameter, but its form generalizes for multiple scaling parameters; see Section 2.2. The singular part that satises the well-known Widom homogeneity relation, implies scaling relations for critical exponents that describe the divergences of various quantities of the system [1], [3]. We briey review the main results here. The distance from the critical temperature, T c , can be dened as Near the critical temperature the singular behavior of the specic heat, spontaneous magnetization, magnetic susceptibility, and response to magnetic eld at T c are given by the exponents α, β, γ and δ: C |τ | −α , m (−τ β ), χ |τ | −γ , and h = |m| δ . The homogeneous singular piece of the free energy for two scaling elds, u 1 = u T , the thermal bonding scaling eld, and u 2 = u h as the external magnetic scaling eld gives for Eq. (25), where y T , y h are parameters that characterize the homogeneous function of degree d. Kadano [4] showed that d is the dimensionality of the lattice and in agreement with analysis of the renormalization approach to our models as discussed above.
First we can obtain the critical exponent of the magnetization, β, by dier- The thermal scaling eld u T is proportional to τ which yields Next, the degree of the critical isotherm, δ, is obtained by dierentiating Eq. (27) with respect to u h and setting u T = 0 and b = u Noting u h is proportional to the ordering eld h, we nd that Magnetic susceptibility is found through dierentiation of Eq. (27) twice with Hence, the critical exponent for susceptibility reads Finally, the critical exponent of the specic heat at constant eld is given by Again, setting u h = 0 and b = ( and, we conclude that The four critical exponents are obtained through the scaling parameters y T , y h . Combining the expressions for the critical exponents yields relationships between the four. Consider the combination of Eqs. (30), (32), and (35) which give From Eqs. (30), (32), and (38) we nd Hence, we obtain exact relationships for the critical exponents using scaling theory.
An identication used in getting the critical exponents came from operations . This choice of b is related to the fact that at critical points the correlation length of the system becomes innite; the uctuations in the system become correlated over all distances. In other words, the correlation length is observed to go as (τ ) −1/y T so we choose b on the order of the same length as the correlation length. We can think of repeatedly applying renormalization transformations to an innite system in the neighborhood of a critical point until the correlation length is of order unity.

Finite-Size Scaling
We can also obtain nite-size scaling relations [5]. We consider the inverse system size itself as a scaling eld. This analysis yields the specic heat and susceptibility as a function of system size. Consider Eq. (27), but with a third scaling eld, u 3 = 1/a, with a as the total linear system size measured in lattice units. If the system is scaled by b the inverse system size transforms as a −1 → ba −1 , Similar to the behavior of u T , u h at the critical point, the inverse system size also tends to zero through dedecoration. Following the same procedure of obtaining the specic heat and susceptibility previously will lead to results like Eqs. (33) and (36) The exponents on the right hand side of these equations can be positive or negative corresponding to a divergence in the innite system limit or approach to a nite limit.
List of References [1]

Classical Interpretation of the Ising Model
It is well-known that below the critical temperature microscopic spins can align over macroscopic distances in some systems. This is known as long-range order and in these systems it can produce macroscopic spontaneous magnetization.
Above the critical temperature the spins display only short-range order which produces no net macroscopic arrangement. The Ising model can be used to help describe this process seen in nature [2], [3], [4], [5]. Recall that Hamiltonian of L l is reduced to dimensionless form by absorbing β = −1/k B T in the coupling constants κ and h. The reduced energy of a spin conguration S l is In this sum i, j runs over all pairs of nearest-neighbors of the lattice L l and i runs over the sites. The coupling constant κ = −β determines the nearestneighbor pair interaction coupling for nearest neighbors and h is the coupling with the external magnetic eld. The partition function is given by where each s i in S l assumes the values ±1.

The Linear Ising Chain
In the linear Ising model the system is constructed as a periodic chain of dipoles; see Figures 1,2,3. The energy for a conguration S l is with the periodic boundary condition s n s l +1 = s 1 . The Hamiltonian in the absence of a eld simplies to pair-coupling energy only and the transfer matrix becomes Squaring the transfer matrix corresponds to dedecoration, as discussed in the previous chapter; see section 2.1, This matrix can be written as with g and κ determined by Solving Eqs. (51) and (52) gives The factor g introduced in Eq. (50) follows from the zero of energy convention implicit in the denition of the Hamiltonian.
If L l starts with κ everywhere then one operation of dedecoration, in terms of how coupling pair interactions change under subsequent steps of the dedecoration transformation, can be dened as Repeated applications of dedecoration continue to apply the same type of interaction transformations. Successive transformations can be obtained by recursion, Here, k is a positive integer that denotes the number of times the interactions transformed, i.e., a second dedecoration operation would give κ 2 = κ 1 [κ 1 (κ)] = κ 1 • κ 1 (κ). In the subscript notation applied to pair interactions it is understood that κ 0 = κ.
One application of dedecoration performs a summation over half of the sites.
The g term to the power of remaining bonds of the dedecorated lattice gives the overall energy shift, G, and the partition function according to Eqs. (53) and (54) can be constructed as To calculate the partition function by recursion from starting lattice L l we repeatedly square the transfer matrix. Any nite chain can be reduced by a nite number of dedecoration transformations to a system small enough to calculate the partition function directly. If the chain was innitely sized, then the process described in this section could be applied repeatedly until a quantity such as the free energy per site has converged to a desired accuracy. A property of the zero-eld linear Ising chain is that This property of the function κ shows that there is no phase transitions in the linear system. The only xed points in this system are the innite and zero temperature limit.

Linear Ising Lattice Free Energy, Total Energy, and Heat Capacity
As mentioned before, we use the term free energy as short for the reduced free energy, i.e., the free energy multiplied by −1/k B T ; see Section 2.1.1. The total energy and heat capacity that follow from our denition of the free energy are therefore not given by their usual expressions and are also in reduced form.
Using the reduced forms of these quantities does not alter the results presented in this dissertation other than by scaling in various powers of temperature, which is of no signicance for our purpose.
The rst derivative of free energy in the linear chain conguration for level l with respect to κ and no external eld gives the total energy U as We consider four systems, n b l = 4, 8, 16, and 32 bonds, representing L 2 , L 3 , L 4 , and L 5 ; the results of the energy are plotted in Figure 9. Another derivative with respect to κ gives the heat capacity of the system, C, which in the reduced quantities we use is nothing but the variance of the reduced energy, see Figure 10, Notice, the energy and heat capacity are constructed as a series that can be

The Hierarchical Ising Diamond
Here we consider the hierarchical diamond lattices in Figures 4, 5, and 6. The Hamiltonian with bond interaction κ and eld interaction h is given in Section 3.1; see Eq. (44). The transfer matrix method is applied as explained in Section 2.1, Eq. (6).

The Hierarchal Ising Lattice Phase Transition
In the absence of an external eld Eq. (44) reduces to Pair-site dedecoration for level l results in bond doubling in this case; see Figure 8.
Following Section 3.2 for how interactions change under subsequent steps of the dedecoration transformation, we similarly write how interactions change in the hierarchical diamond lattice as The Starting dedecoration above or below the critical point will tend the interactions towards the trivial xed points at zero or innite temperature. The problem has three xed points in total now, two are trivial, and the third, κ c , is the desired xed point associated with the phase transition.

Hierarchical Diamond Lattice Total Energy and Heat Capacity
The total energy and heat capacity of the hierarchical Ising diamond lattices are constructed by means of two coupled recursion relations In Figure 11 the energy and heat capacity per bond of the hierarchal Ising diamond lattice for various system sizes are shown. Notice that the specic heat has a rounded cusp at κ c that becomes sharper as system size increases. The linear chain which lacks a phase transition does not have this feature and is expected to be less sensitive to aws in the RNGs that we want to test by means of Monte Carlo simulation. Compare with Figure 10. Of course, in the hierarchical lattice, only the innite Ising system has a true cusp. Figure 11. The energy (pair-coupling) and heat capacity graphs per bond for hierarchal systems L 2 , L 3 , L 4 , L 5 , L 7 , and L 9 with κ c included. Notice as systems increase in size divergent phenomena appears about the critical point.

The Analysis of the Hierarchal Lattice System with Field
The previous sections excluded an external eld in the system. When the interaction K was just pair-coupling κ in the hierarchical models then dedecoration of L l to L l−1 would produce bond-doubling between all neighbors. To correctly account for bond-doubling in computing the partition function of the L l−1 hierarchical model, the general expressions for the shift of energy, g, and renormalized pair-coupling, κ , both needed to be squared; see Section 3.3.1. Upon dedecoration in systems with an external eld, sites with a eld will be summed over; these elds must be conserved and so they will be renormalized, split, and propagate to neighboring sites. The remaining sites with their elds will get contributions of the propagated elds depending on how many neighbors the remaining site originally had. This implies that if one starts with a system with uniform interactions that a single recursion step will destroy this property; this is something that does not happen in the absence of a magnetic eld. Repeated recursion generates a hierarchy of coupling constants, as will be shown in detail below.
We begin by dening the general results of dedecoration in the presence of an external eld; these results are the energy shift g(κ, h), the renormalized bonding κ (κ, h), and the renormalized eld h (κ, h). Until this point we have never considered a site interaction term in K using dedecoration and have only dealt with pair-coupling transformations and the shifts in energy that follow. Therefore, we step back from hierarchical models and look at a variation of our linear Applying dedecoration, summing over site s 2 , the eld transforms, splits, and propagates to the remaining sites. The pair-coupling also transforms and the energy shifts as well. The partition function after dedecoration is Z(κ, h, s 1 , s 3 ) = e χ+κ (s 1 s 3 )+h (s 1 +s 3 ) .
Here, χ is related to g mentioned previously by χ = log g. Now, any arbitrary spin function F can be expanded in terms of interactions between sites as where s a is the product of sites associated with interaction K a . The s a for paircoupling is nearest neighbor sites and for the eld it is single site interaction. The spin products s a form a complete, orthogonal basis in the space of spin functions.
This property can be used to invert Eq. (70) to yield, where n stands for the number of lattice sites. Using Eq. (71) the renormalized interactions are κ (κ, h) = 1 4 h (κ, h) = 1 4 with the explicit representation in relation with how we previously dened these parameters given as The results in Eqs. Application of dedecoration to the most general system of the hierarchical Ising model with an external eld requires the notation h k for the eld site interaction parameters. The subscript k denotes a eld at a site with 2 k connecting neighbors; using this notation, starting lattice L l will have l − 1 elds. Consider lattice L 4 with nearest neighbor interaction κ everywhere and elds h 1 , h 2 , and h 3 ; see Figure 13. The partition function for this system is Z 4 (κ, h 1 , h 2 , h 3 ).
The resulting renormalized elds after dedecoration at the remaining sites are One iteration applied to partition function Z 4 (κ, h 1 , h 2 , h 3 ) gives Recall that G comes from a shift in energy upon dedecoration and this time it is related to the new g(κ, h), i.e. G l−1 = g n b l−1 = g 2 n b l−1 . The general partition function Z l with l − 1 elds after one step of dedecoration is The result of Eq. (82) completes the general analysis of hierarchical models with external elds. In practice we start with uniform magnetic eld, but this property is not invariant under renormalization, which is rather unusual and more general than the method described in previous chapters.

The Complex Constants and Probabilities of Energy and Magnetization
The probability of nding the system in conguration S l is given by The probability of nding the system in a state with magnetization M 0 is Where M = i s i . The δ-function ensures that only those spin congurations that match the desired magnetization, M 0 , survive the sum giving the desired probability. The probability of nding the system in a state with pair-coupling Where I = i,j s i s j . In general Eqs. (84) and (85)  The Kronecker δ-function can be represented as a sum of exponentials. Consider two variables x, y that assume all integer values 0, 1, 2, . . . , N for some max value N ; these integers repeat in range n = N +1. For y xed, the general complex transformation of δ x,y is then given as If the variables δ x,y = δ 0,x−y do not initially assume values 0, 1, 2, . . . , N then they can be shifted and scaled to do so. In the following applications of this transformation all variables denoted with primes are variables that need to be shifted and scaled to satisfy the requirement of being subsequent integers.
We can write Numerically, Eqs. (88), (89) allow us to compute the full probability distributions recursively. See Figure 14 for the result of such a calculation.

List of References
where Θ(p) = 1(0) if p is true(false). Here κ and h are the previous pair and external eld interactions and λ is a second pair interaction, which is generated by the dedecoration operation as soon as the eld h diers from zero and q > 2.
Finally, the parameter w reects the freedom in choice of the zero of energy. It will be xed by convention and will generate the self-energy term in the scaling relation for the free energy. The q = 2 case can be related to the Ising model by means of a simple transcription. The transfer matrix for integral q is constructed as a q × q matrix that can be conveniently expressed in Boltzmann weight form [3] as Both Eqs. (90) and (91) have four parameters.

Potts Lattice Dedecoration in Zero Field
Dedecoration of the Potts diamond Hamiltonian in zero eld reduces Eq. (90) to only one coupling constant, κ. To obtain the conventional form of the Potts model for Eq. (91) in the absence of a eld we choose a = d = 1, and b = c = exp(κ): where g comes from the pair-site dedecoration energy shift and b is the renormalized Boltzmann weighted pair-coupling parameter. Similar to the Ising analysis, see Section 3.3.1, we can write b 1 (b, q) = b (b, q). The transformation has a xed with Note that we can generalize the model and consider q as a continuous variable.
We have not explored the simulation of lattices with continuous q and will refrain from doing so. The method we use, and as explained in detail below, does not lend itself for that purpose.

Hierarchical Diamond Lattice Total Energy and Heat Capacity
The partition function, the energy, and the heat capacity of the hierarchical q-state Potts model in zero eld are conveniently expressed in Boltzmann weight form. Let Z q l (b) be the partition of the hierarchical lattice of level l, state q. The partition function satises the following recursion relation The reduced free energy F q l = logZ q l satises the scaling relation From the free energy the thermodynamic quantities of interest follow by taking derivatives of F with respect to b, where all superscripts in parentheses denote rst and second derivatives with respect to the arguments b and b 1 . With Eqs. (100) and (101), the total energy and heat capacity of lattice l and state parameter q are To compare to Monte Carlo simulations in Chapter 5 we once again only use nite lattices. In calculations, the derivatives of the free energy can all be predened and calculated numerically; see Appendix A.
The q-state Potts model has specic heat that diverges at the critical point for q suciently large. We expect that this makes Monte Carlo computations for large values of q more sensitive to correlations in the RNGs. Figures 15 and 16 shows results for both the total energy and heat capacity per bond with q = 400 for various large nite systems. Figure 17 shows the q = 2 Ising results for the heat capacity per bond for the same nite systems for comparison. . This gure uses q = 400 and range the lattice level from l = 6 ,..., 12. The xed value b c for this selection of q is b c = 0.016615.

Hierarchical Potts Diamond Lattice in a Field
We now consider the RG transformation of the hierarchical Potts model in a eld. Using Eq. (90) and choosing w = 0, the following transfer matrix can be seen to have a structure that is invariant under renormalization, Recall that h 1 denotes the renormalized eld interaction at a remaining site after g(κ, λ, h, q) = 1 + (q − 2) e 2κ + e h+2κ+2λ , κ (κ, λ, h, q) = log e κ (2 + (q − 3) e κ + e h+κ+2λ ) 1 + (q − 2) e 2κ + e h+2κ+2λ , Combining Eqs. (106) and (110) gives the basic single site renormalization group equation of how sites transform in both the linear and hierarchical Potts lattices, h (κ, λ, h, q We remind the reader that pair-site dedecoration in hierarchical diamond lattices require the energy shift to be squared, and the pair-coupling to be doubled; this gives g = g 2 , κ = 2κ , and λ = 2λ . Additionally, we point out to the reader that the connectivity of neighbors in the hierarchical Potts diamond lattice is exactly the same as the hierarchical Ising diamond lattice. This implies that where G q l−1 (κ, λ, h) = g n b l (κ, λ, h, q).

Hierarchical q-state Potts Model Probability Distributions
Using the partition function Eq. (112), the probability distributions of the hierarchical Potts diamond lattices can be constructed closely following the procedure used for the Ising model. We restrict analysis to q ≥ 3 in what follows. The pair-coupling energy, I, with zero-eld, from the Hamiltonian of Eq. (90) is The probability of nding the system in a state with energy I 0 is For q ≥ 3 the energy distribution increases by steps of one. From Section 3.3.4 the complex transformation of the Kronecker delta allows Eq. (114) to be written as with the partition functions given by Eq. (112) for zero-eld. The order-like parameter, M , from Hamiltonian Eq. (90) is The probability of nding the system in a state with given value of M 0 is For q ≥ 3 this distribution also takes steps of one. By Section 3.3.4 the complex transformation then gives with the partition functions given by Eq. (112). Note, if the eld is taken as zero in Eq. (118), the shift acting on the eld interaction term still remains. In Figure   18 we show probabilities for the L 4 lattice with q = 3 evaluated at the critical point using Eqs. (115) and (118). In practice, we apply recursion to the partition functions until the L 2 lattice. The L 2 lattice is a sucient stopping point because it is the rst lattice in the dedecoration process where all site parameters will be equal to each other no matter what the starting conditions; brute force calculation of the partition function is easily applicable at the L 2 lattice. Note, computing the corresponding partition function of the L 2 lattice by brute force counting will still take time proportional to q 4 . This becomes long for large q and can be avoided by writing the partition function for the L 2 lattice as a polynomial in q; see Appendix B for how this is done for zero-eld. Generators (PRNGs). There are many dierent types of these generators [1].
A PRNG after being appropriately seeded, recursively generates a sequence of numbers that have the properties of random numbers. We shall not attempt to provide a denitive denition of what a random number is; let us just say that the sequence of numbers (r 1 , r 2 , ..., r n ) should contain no information about r n+1 for reasonable choice of n. In pseudo-random number generators, r n+1 is always determined by its predecessors and therefore cannot be a truly random number.
Although the numbers produced by PRNGs are pseudo-random numbers we will simply refer to both the generator and output as Random Number Generator (RNG) and random numbers" instead. If these random numbers are suciently correlated, the Monte Carlo will have systematic errors in addition to unavoidable statistical ones.
The bias introduced by correlated random numbers is not a problem when it is smaller than the statistical error of a particular computation. However, as more compute cycles become available, these statistical errors decrease so that the demands on the quality of the RNG increase. In practice, RNGs are subjected to various test suites [2] that have been developed in the past. We are interested in adding yet another random number generator test to the standard battery.
Thermal and magnetic quantities in the hierarchical q-state Potts diamond lattice can be expressed as sums of correlation functions that diverge on approach of the critical point in an innite system. In the nite system, which are the only ones where Monte Carlo is applicable, these quantities have no divergences.
Nonetheless, because of the incipient divergences, nite system near criticality are highly sensitive to correlations introduced by a bad random number generator. It is this property we intend to exploit.
In this chapter we discuss the details of a Monte Carlo algorithm with which one can compute estimates of observables like the heat capacity. These estimates are to be compared with the theoretical values computed by the methods discussed in the previous chapters. The probability distributions discussed previously are also estimated by constructing histograms.

Monte Carlo Simulation of The Hierarchical q-state Potts Diamond Lattice by Direct Sampling
To sample the Boltzmann distribution one can use the Metropolis-Hastings (MH) algorithm. However, there is an inherent problem in using such algorithms, a problem that is shared by other more sophisticated algorithms that have been designed with the problem of critical slowing down [3]. If any non-trivial Markov process is used to sample the Boltzmann distribution, subsequent congurations will be correlated. This typically reduces the eciency of a calculation, while also, burn-in is required because an initial conguration may be chosen with a probability that does not correspond to equilibrium.
Hierarchical lattices have the unique feature that the congurations they support can be sampled directly by the decoration procedure that generates the lattice geometry. At each step of decoration, at most, two frozen nearest neighbors determine the probability of the state of the site added next. In this way one can generate lattice congurations of arbitrary size.

The Direct Sampling Algorithm for the Zero Field Hierarchical q-State Potts Diamond lattice
The very rst site of our lattice of Potts variables, s 1 , has no neighbors, so that each state has equal probability P(s 1 ) = 1 q .
Once the rst site is generated according to Eq. (119) the second site is generated in the absence of a magnetic eld with probability Here, we remind the reader that b = e κ is the Boltzmann weight associated with The Implementation of the Direct Sampling Algorithm This subsection discusses how the sampling described in the previous section is implemented in detail; the reader may skip to section 5.2 to continue with the standard error analysis techniques we propose for the RNG test. Consider readily available the sequence of random numbers (r 1 , r 2 , ..., r n ) from the U(0,1) distribution.
1. Use r 1 to determine state of s 1 with Eq. (119).
2. Use r 2 against Eq. (120) to determine state of s 2 • If r 2 > P(s 2 = s 1 |s 1 ) then state of s 2 is same as s 1 .
If r 3 > P(s 3 = s 1 |s 1 , s 2 ) then state of s 3 is same as s 1 ,s 2 .
If r 4 > P(s 3 = s 1 |s 1 , s 2 ) then state of s 3 is same as s 1 or s 2 .
All newly included sites after s 2 will always only have two nearest neighbors like s 3 .

Standard Error Analysis
Using Monte Carlo we want to calculate the expectation value Q of some observable Q. The statistical ensemble average of a system of all states allowed, µ, denes Q , However, we assume that during a Monte Carlo run the system will not pass through every possible state. With MC direct sampling we generate congurations with a probability given by the Boltzmann distribution [4]. The ensemble average of the expectation value can then be written as a time average. The unbiased estimator, Q, of this expectation value, Q , is Direct sampling yields independent realizations of observable Q, and for large t, Q is normally distributed around Q according to the central limit theorem. The probability that Q exceeds deviation ∆ from Q is These standard techniques can be applied to any observable to determine whether observed estimates deviate signicantly from their expected values.

The Probability Distributions and χ 2
To estimate whether the estimates obtained for a whole probability distribution agree with what is theoretically expected, the usual χ 2 statistic is an obvious choice: If the x k are independent, standard normal stochastic variables, χ 2 is distributed according to the χ 2 distribution of w degrees of freedom. Consider during an MC run consisting of t observations we construct a histogram. The number of times that an observable assumes the value associated with the kth bin of this histogram is t k . The corresponding distribution is the binomial distribution with average p k t and variance p k t(1 − p k ). Here, p k is the probability of landing in bin k. The chi-square statistic is then The problem in applying this statistic to the bins is that many correspond to extremely improbable events. The condition of normality for applicability of the χ 2 -distribution is not satised for such bins. To deal with this problem we introduce super-bins formed by combining suciently many neighboring bins until the combined bin count exceeds the expected RMS error by a sizable factor for each bin. In practice, we use the following condition as it guarantees a normal distribution of the bin counts. We dene N [µ, σ 2 ] as the normal distribution with average µ and variance σ 2 . Each term in Eq. (127)   Consider the free energy and rst derivative of the recursive free energy of the q-state Potts model with respect to internal argument b. We drop the superscript dependence on q to simplify the notation, We write Eq. (A.2) in the following form,  In an innite system the nal term of the sum would be pushed innitely far away and one could approximate the energy by including more J terms in the sum.
However, in a nite system the dedecoration and therefore the sum of terms must end. The L 1 two site lattice is the nal stop point of nite dedecoration. Recursion with Eq. (A.7) in a nite system gives F (1) The nal term is the brute force derivative of the free energy of the L 1 system with respect to b k which is the nal dedecorated coupling parameter of the lattice, F The total energy of the system is now written in terms of derivatives of quantities that can be pre-dened. The partial derivatives, ∂b 1 /∂b, ∂b 2 /∂b 1 , etc, are all the same derivative of b 1 w.r.t. to b evaluated at a new Boltzmann pair-coupling point.
Recall, the renormalized expression for b 1 , (A.10) Taking the derivative of Eq. (A.10) w.r.t. b, which is evaluated at b to obtain the numeric derivative value. All partial derivatives of the new Boltzmann pair-coupling w.r.t. to the old Boltzmann pair-coupling will be of the same form. Next, the partial derivatives of the coecients ∂J(b)/∂b, ∂J(b 1 )/∂b 1 , etc. Recall, the coecient form, J l (b) = n b l−1 log g (b, q). Take the log, then dierentiate w.r.t. to b, , (A.14) and then evaluate at numeric point b. The quantity Eq. (A.14) will be the same form for all levels of dedecoration. Finally, the brute force partition function of the q-state Potts model at L 1 and arbitrary renormalized dedecorated coupling is and can be resolved for b k , q as, Z 1 (b k , q) = q 2 + q(b k − 1).