Simple model for dice loading

Dice tossing is commonly believed to be random. However, throwing a fair cube is a dissipative process that is well described by deterministic classical mechanics. In Nagler and Richter (2008 Phys. Rev. E 78 036207; featured in 2008 Nature 455 434), we proposed a simplified model to analyze the origin of the pseudorandomness: a barbell with two masses at its tips with only two final outcomes. In order to keep things simple, we focused on the symmetrical case of equal masses. Here, we complete the picture by considering the general asymmetric case of unequal masses. We show how, depending on the initial conditions, dissipation during bounces, and mass asymmetry, the degree of unpredictability varies. Our analysis reveals, for the simplest possible non-trivial dice throwing model, the effect of dice loading. A surprising consequence of dynamical resonances is that an experienced player may benefit sometimes more from an unloaded than from a loaded barbell. In addition, we investigate the influence of loading on the symmetry breaking process causing one mass to come to rest earlier than the others.


Introduction
In an attempt to understand the nature of randomness in dice tossing, we recently analyzed a simple deterministic model [1] where apparent chaos is generated in the succession of free flight episodes and bounces off the ground. Our 'die' was only a caricature of a realistic cube [2]; in fact, it was a barbell consisting of a massless rod and point masses m 1 and m 2 attached to its tips. This barbell was allowed to perform translational and rotational motion in an (x, y)plane, with gravity acting in the negative y-direction and y = 0 being the ground. If released at some height above the ground, with a given initial orientation and angular velocity, the rod flies freely and without dissipation until one of the two masses hits the ground. At this point an inelastic reflection takes place, which throws back the barbell in the upward direction but with less energy than before. After a number of such bounces the barbell will come to rest in one of two orientations, mass 1 being either to the right or to the left of mass 2. These two final states are the system's two attractors. In addition, the system possesses two unstable, or hyperbolic, final states where the barbell stands upright on the ground, with either mass 1 or mass 2 on top.
We analyzed this model in terms of Poincaré sections and final state diagrams. The corresponding pictures contained structures reminiscent of classical chaos, with 'almost fractal' boundaries between the two basins of attraction [3]. The fractal nature is not complete as it would be in the complete absence of dissipation, because the loss of energy in each bounce prevents the stable and unstable manifolds of the two hyperbolic states from developing the infinitely fine entanglement of Hamiltonian horseshoe structures. However, with increasing initial height and decreasing inelasticity of the bounces, the phase space partitioning looks more and more chaotic. In this sense, randomness of the barbell's dynamics becomes a matter of degree and of the initial location in phase space. We believe the same holds true for realistic dice tossing.
A couple of reports on scientifically motivated experiments to test the fairness of dice and coins are worth mentioning [4]. Vulvović and Prange [5] compared the basins of attraction of two possible final configurations of a homogeneous rod with those calculated by Keller [6] for a simplified coin tossing model. Feldberg et al [7] modeled a rolling square and studied the corresponding final states for different initial conditions. Ford [8] and Zhang [9] studied the tosses of a homogeneous disc, and Murray and Teare [10] estimated the probability of a tossed coin landing on its edge.
In the present paper, we extend the analysis to include the effects of loading the barbell in terms of unequal masses. For that purpose, we found it appropriate to define two new kinds of diagrams: orientation flip diagrams (OFDs) and bounce diagrams (BDs). The OFDs are Example of a trajectory and its symbolic characterization. In contrast to [1], we define the symbols for bounces relative to the initial orientation of the barbell. Two possible orientations are distinguished, mass 1 being to the left or to the right of mass 2. A bounce is labeled by 0 when the orientation at time t = t bounce is the same as that initially at t = 0; a bounce where the barbell has flipped its orientation compared to the initial one is denoted by 1. Note that the symbolic code does not reveal which ball bounces. In addition, when during a flight the barbell passes the vertical axes clockwise, the passage is denoted by R. A counterclockwise passage is denoted by L. No further symbols are added to the symbol code when the barbell has too little energy for vertical passages, that is, when the motion is restricted to either infinitely many type-0 or type-1 bounces. In addition, we omit repeating symbols at the end of the code and do not consider the (singular) case of the barbell hopping vertically with vanishing angular momentum. The symbolic code for the example is R1L00R111L0.
closely related to the previous final state diagrams; the only difference is that the final state is not considered as such but in its relation to the initial state: has the orientation changed or not? The BDs, on the other hand, contain a new kind of information that was suggested to be relevant by observing the loaded barbell's behavior: which of the two masses bounces more often before the final state is reached? As will be seen, this feature has a stronger dependence on mass differences than the orientation flip.

The model
Consider a barbell in the two-dimensional (x, y)-plane, see figure 1. Two point masses are connected through a massless rod of unit length. Their positions are denoted by (x 1 , y 1 ) and (x 2 , y 2 ), respectively. Furthermore, we assume gravity to pull in the negative y-direction, the floor being at y = 0. The barbell is released from an initial height with an initial orientation ϕ and given initial velocities; it then follows the laws of Newtonian mechanics until one of the masses hits the ground. The reflection process requires special consideration. The theory of this motion is given in [1] and need not be repeated here. We just mention the important facts.
To represent the three degrees of freedom (two for translation and one for rotation), we choose the center-of-mass coordinates (x, y) and the orientation angle ϕ (see figure 1). We assume that upon reflection no momentum is exchanged in the x-direction; hence,ẋ is a constant that we take to be 0, together with x. Thus, the system essentially has only two degrees of 4 freedom, y and ϕ. If we use appropriate units of length, time and energy, the total energy is given by where β 1,2 are the mass ratios β 1 = m 1 /(m 1 + m 2 ) and β 2 = m 2 /(m 1 + m 2 ). Since β 1 + β 2 = 1 there is only one independent mass parameter β = β 1 , but we find it convenient to also use β 2 for symmetry in the notation. It suffices to restrict the analysis to 1/2 β < 1. An important energy mark is E crit ≡ min(β 1 , β 2 ) = 1 − β 1 , the critical energy: at smaller values of E, no change of the barbell's orientation is possible.

Orientation flip diagrams
We analyze the system in the following manner. At a given height h above the ground, the barbell is released with zero velocityẏ but arbitrary orientation ϕ and angular velocityφ. This produces a two-dimensional set of initial conditions (ϕ 0 ,φ 0 ), the energy E 0 = h + 1 2 β 1 β 2φ 2 0 increasing withφ 0 . In our OFD, see figure 2, a fine grid of these initial conditions is scanned, and for each initial point, the motion is computed until the final state is determined. Then, in contrast to the use of a color code for the final state as in [1], we choose the color according to whether the orientation in the final state is the same as the initial one (yellow color, state 0) or different (red color, state 1), see figure 1. If both masses were equal as in [1], this would introduce a symmetry under ϕ → ϕ + π into the final state diagrams there; the upper left panel of figure 3, where β = 1/2, exhibits this symmetry. The effect of loading the barbell can be clearly observed now: this symmetry gets broken.
The extent to which this happens is determined by the intersection of the plane of initial conditions with the stable and unstable manifolds of the hyperbolic planes (ϕ,φ) = (π/2, 0) with energy E β 1 , and (ϕ,φ) = (3π/2, 0) with energy E β 2 . These intersections are the white lines (stable manifolds) and black lines (unstable manifolds) in figure 2. The stable  manifolds are the boundaries between the basins of attraction, and the unstable manifolds are their mirror images under (ϕ,φ) → (ϕ, −φ). The linear approximation to these manifolds in the neighborhood of the hyperbolic planes was computed in [1]. The result for the eigenvectors 7 at A = (ϕ,φ) = (π/2, 0) (small mass at the bottom) is δϕ δφ and at B = (ϕ,φ) = (3π/2, 0) (large mass at the bottom) δϕ δφ These directions depend on the initial height h but not on the dissipation parameter f ; therefore the white and black lines are omitted in the lower two panels of figure 2. It is remarkable how far the linear approximation extends from the hyperbolic points. The dependence on β has the effect of turning the central rhomboid at β = 1/2 into a deltoid at β > 1/2. The implication for 'playing the loaded barbell' is evident: choosing the initial condition such that the angle ϕ 0 is near π/2 (the heavy mass on the top) makes it more probable to end up in the original orientation than near 3π/2. However, this statement only holds for relatively smallφ 0 . At larger initial velocities (and low friction), the mixing of yellow and red points rapidly becomes too intricate for the loading to be helpful. With larger friction, the patterns look more complicated but also more promising for a player who cares to practice a lot.
While figure 2 shows the dependence of OFDs on the dissipation parameter f at given mass β = 0.8, we show in the left parts of figures 3 and 4 how the pattern changes with β, at fixed f = 0.2. The asymmetry in the slope of the invariant manifolds becomes more and more pronounced, and the very chaotic mixing of colors is shifted to higher and higher values of the initial angular velocities.
It is interesting to ask how the probability P flip (E 0 ) of final state 1 (orientation has flipped) depends on the initial energy E 0 . To answer that question we sampled trajectories for random initial conditions at fixed initial energy E 0 . For each trajectory we chose the center-of-mass velocity to be zero,ẏ(t 0 ) = 0, initial orientation ϕ(t 0 ) to be random and uniformly distributed between 0 and 2π, the initial altitude h 0 also to be random and uniformly distributed between 1 and E 0 , whence the initial angular velocity must beφ(t 0 ) = + √ 2(E 0 − h 0 )/(β 1 β 2 ) (the other sign would give identical results). For friction strengths between f = 0.1 and f = 0.5, the results are shown in figure 5 (β = 0.8) and in figure 6 (β = 0.5).
As expected, all curves start from P flip = 0 at E 0 = 1 and approach the asymptotic value P flip → 1/2 as E 0 → ∞. More interestingly, at moderate initial energies they perform striking oscillations, see the insets of the figures. These oscillations depend both on friction f and mass ratio β; they are a consequence of dynamical resonances. Recall from (4) that the barbell can lose three quarters of the initial energy at the first bounce if f = 0.5. Hence, even for E 0 = 8, it may perform only two or three bounces before it can no longer flip over. Depending on the state before the last flip, for a given E 0 , slightly higher initial energies may or may not increase the flip probability. Since the final configuration is evaluated relative to the initial one, a larger E 0 can therefore both decrease or increase P flip (E 0 ). As a consequence, P flip (E 0 ) depends nonmonotonously on E 0 and so the 'phase' and 'frequency' of the oscillations depend on both β and f .

Bounce diagrams
In the course of our computer experimentation, we noticed that with unequal masses, the light and heavy barbell tips tend to have different numbers of bounces before the final state is determined, with an increasing preference for the heavy mass as β grows. This may not be relevant for the ordinary tossing game because this property is independent of which final state is reached, but it does reflect an interesting aspect of the loading and might be used for another kind of game.
As it may not be obvious that the two masses can have different numbers of bounces on average, let us demonstrate this fact with a simple calculation for the case where both of them stay close to the floor, i.e. for small energy E 1 and |ϕ| 1 (the case where ϕ is close to π would give the same result). Using the approximation sin ϕ ≈ ϕ, we obtain Together with the energy formula, equation (1), this yields where E 1,2 (y 1,2 ,ẏ 1,2 ) = 1 2 β 1,2ẏ 2 1,2 + β 1,2 y 1,2 .
This shows that the Hamiltonian separates into one of two independent freely falling mass points. Ifẏ 0 is the initial velocity of any of them immediately after a bounce, then T = 2ẏ 0 is the time at which the next bounce occurs and A =ẏ 2 0 /2 is the amplitude reached in between. Combining this with the friction modelẏ →ẏ = (1 − f )ẏ, we see that each mass point will perform oscillations with geometrically decreasing amplitude and period, If the initial T is T 1,2 for the two masses, then they come to complete rest, after an infinite number of bounces, at times T 1 / f and T 2 / f , respectively. If we stop counting at some earlier time, then we will find that the mass with smaller T i and amplitude A i has undergone more bounces than the others. It seems intuitively clear that the heavier mass dominates the relaxation process: the reflection has a bigger effect on the lighter mass than vice versa. Therefore, one expects the larger mass to come to rest earlier than the smaller and, hence, to exhibit more bounces before the final state is determined. This is indeed what the BDs in the right panels of figures 3 and 4 show. The same (ϕ,φ)-plane as in the OFDs is scanned, but now we monitor the number of bounces of the two masses before the orientation can no longer change. If mass 1 bounces more often, the color is white to gray; if mass 2 bounces more often the color is red to gray; the brightness reflects the number of bounces as indicated in the color code bars on the upper right of each panel. At β = 0.5, the numbers on the whole are equal, but as β grows, bounces of the heavy mass become more and more dominant.
In figure 7, we plot the fraction F h of initial configurations in the BDs for which the heavy mass bounces more often, for three different friction strengths. The data points are determined by pixel counting in diagrams like those of figures 3 and 4, for the respective values of β. The result of this numerical evaluation is that F h is a little lower than β. It appears that the deviations from F h = β originate from the low-energy part of the BDs.  In order to gain an understanding of the linear law, we assume a sort of Boltzmann statistics in the spirit of Levin [11]. Consider the two unstable upright configurations of the barbell, with energies E A = E crit = 1 − β (at ϕ = ϕ A = 3π/2, heavy mass at the bottom) and E B = β (at ϕ = ϕ B = π/2, light mass at the bottom), respectively. These states are the centers of channels on the barbell's route from E = E 0 to the final state at E = 0. When the barbell is in the upright configuration A with the heavy mass touching the floor, the heavy mass obviously bounces more frequently than the light mass. Hence, we estimate the fraction of white points in the BDs, figures 3 and 4, as F h (β) ≈ P A (β), where P A (β) is the probability that the final state is reached through channel A. A rough approximation to P A (β) may be derived from the assumption of a heat bath in which the barbell's energy fluctuates on a 'thermal' energy scale kT ≈ (E A + E B )/2 = 1/2. This implies Boltzmann factors Together with this yields 12 An expansion of equation (13) around β = 1/2 yields P A (β) = β. In figure 7, the gray line represents this simple behavior, whereas the black line follows equation (13). The agreement with the data is satisfactory, given the roughness of the argument. A slightly more involved argument would take into account a finite size of the two channels and integrate over certain ranges around ϕ A and ϕ B , say of size π/2. This would lead to assuming The numerical evaluation with kT = 1/2 gives a curve that is very similar to the black curve of figure 7.

Summary and conclusions
We investigated a simple model for throwing loaded dice in terms of a barbell with unequal masses at its tips. In comparison to our previous study with equal masses, the emphasis here is on the effects of loading. This required the introduction and statistical evaluation of two new kinds of diagram, orientation flip diagrams (OFDs) and bounce diagrams (BDs). OFDs coding for the relative orientation of final and initial states show that the flip probability is an intricate feature. Its overall behavior, as shown in figures 5 and 6, is to increase from 0 to 0.5 as the initial energy increases, but with complicated oscillations superimposed, which depend on both the friction and the mass ratio. Surprisingly, these oscillations can be bigger for equal masses than for unequal masses, so an experienced player may benefit more from an unloaded barbell than from a slightly loaded barbell, by choosing an optimal initial energy. The bounce diagrams of figures 3 and 4 reveal a strong tendency, increasing with mass difference, for the heavier mass to bounce more frequently than the lighter one, see figure 7. This feature is unrelated to the final state, but it suggests that a player may benefit from loaded dice by choosing an optimal initial angle.
The experience gathered with this simple model may help us to analyze more realistic three-dimensional bodies such as cubes.