He told me later that he looked online for loaded dice, that is, dice that would favor some of the faces and not have a uniform distribution. After a while, we thought it would be a nice problem to design a loaded die and to somehow calculate the theoretical probabilities for each face. This is a rather complex problem, so I started to think on reasonable simplifications to get something sound.
The main idea is to load a die by changing the center of mass of it. This would definitely impact the chances of getting a particular face. Usually we get that a regular die has a probability of $1/6$ of displaying a given face based on symmetry arguments. But when talking about a loaded die things get loaded.
Instead of considering a 3D die, I tried to analyze the 2D case, by considering a unit square with uniform mass density $\sigma$ and with a point mass $\mu$ at a point $(a,b)$. This will make the square to have its center of mass at
$$\left(\frac{a \mu +\frac{\sigma }{2}}{\mu +\sigma },\frac{b \mu +\frac{\sigma }{2}}{\mu +\sigma }\right)\,,$$
and not at the center of the square. With this, the symmetries of the square are broken and a new description of its behavior has to be made.
In order to figure out the probabilities, first we have to define what do we mean by the probability of getting a specific face. One way of doing this is by considering the frequentists point of view of the system: lets roll the die many times and calculate the probabilities as relative frequencies of the outcomes.
Being this a a posteriori point of view, analyzing the physics of such situation becomes the best approach.
We want to have the square freely rotating, landing on a flat line, and then recording which face comes up. In order to do this, we have to prescribe an angular velocity of the rotating square, the landing vertex, and the incidence angle. We can simplify the model by just considering the picture when the die lands on the line, forgetting about its spinning past.
The argument is that by ignoring air resistance, friction, and everything else a good physicist will ignore, the landing angular velocity the square will land with is just a scalar multiple of its original velocity in the air (using conservation of momentum). So a good toy model is starting with a vertex fixed on the line, with an initial angle $\theta_0$ and an initial angular velocity $\omega_0$. Using Newton's second law, we have that the sum of all torques is equal to the moment of inertia times the angular acceleration of the square. That is
$$I\frac{d^2\theta(t)}{dt^2}=\sum \tau\,,$$
were the moment of inertia for the square in this configuration is given by
$$I=\frac{\sigma}{3}+\mu(a^2+b^2)\,,$$
and where we only have one torque acting on the square due to gravity. Putting these together, we have that the system is described by the ODE
$$\frac{d^2\theta(t)}{dt^2}=-\frac{g || r ||}{I}\cos\left(\theta(t)+\delta\right)\,,$$
with initial conditions $\theta(0)=\theta_0$ and $\theta'(0)=\omega_0$, and where $g$ is gravity, $r$ the position of the center of mass of the square, and $\delta$ its argument.
This is a nonlinear differential equation which can be numerically solved in mathematica. Solving this equation gives us the rotation angle of the square in time. If the angle reaches $\pi/2$, it means that it completely tilted over the next face, and we would have to make the same analysis again but now with new initial conditions $\Theta(0)=0$ and $\omega(0)=\kappa \theta'(T)$, where $T$ is the time that took the previous face to completely reach the next one and $\kappa$ is a constant describing the energy loss the square experiences when landing on the next face.
We can keep this process as long as the equation
$$\theta(T)=\frac{\pi}{2}$$
has positive real solutions and the number of times this procedure is carries over will tell us on which face the square is going to stop rolling.
This approach is completely deterministic, as given any pair of initial conditions $(\theta_0,\omega_0)$, we can calculate on which face the square is going to stop rolling. In order to introduce probabilities, we can asume a probability distribution on the initial angles and the initial velocities.
It is reasonable to pick a uniform distribution of the incidence angles $\theta_0$, so I choose a uniform probability distribution on $(0,\pi/2)$. As for the initial velocities, the argument relies on the system being cyclic. It is true that all angular velocities are not equally probable, as when someone grabs an actual die the angular velocities would have a distribution centered towards a preferred value, but when considering the problem only from the landing perspective, what really matters is that the system presents an almost periodic behavior on the initial velocity. So we can take a uniform distribution on $(0,\omega_{\text{max}})\,.$ We also have to take into account that the die could rol clockwise or anti-clockwise with uniform probability.
Taking these assumptions into account would make the final probability not to depend on how we roll the die, that is, we are making it up for the lack of the symmetry the system has.
Finally we have to account for the probability of the die landing on a particular vertex. In order to calculate this, we can think that the probability of it landing on a given vertex is proportional to the time this vertex is the closest one to the landing line, while the die is rotating in the air. We can use the moments of mass about the face, and then calculate the probability of a vertex to be the average of its adjacent faces.
Performing some numerical simulations, we get the following probabilities ($\sigma=1, \mu=100$)
| (a,b) | 0 | 1 | 2 | 3 |
| (0.5,0) | 0.268748 | 0.34257 | 0.0461125 | 0.34257 |
| (0.1,0.1) | 0.468528 | 0.468551 | 0.031472 | 0.0314486 |
| (0.5,0.5) | 0.25 | 0.25 | 0.25 | 0.25 |
and here is the code I used in Mathematica (the last function gives a the probability vector {0,1,2,3}):
T[{a_, b_}] = {1 - b, a};
INC[a_, b_, \[Mu]_, \[Sigma]_, c_] := \[Sigma]/
3 + \[Mu] Norm[Nest[T, {a, b}, c]]^2;
Sol[a_, b_, \[Mu]_, \[Sigma]_, \[Omega]0_, \[Theta]0_, c_] :=
NDSolve[{\[Theta]''[
t] == -((9.8 Norm[CM[Nest[T, {a, b}, c], \[Mu], \[Sigma]]])/
INC[a, b, \[Mu], \[Sigma], c])
Cos[\[Theta][t] +
ArcTan[Nest[T, {a, b}, c][[2]]/
Nest[T, {a, b}, c][[1]]]], \[Theta][
0] == \[Theta]0, \[Theta]'[0] == \[Omega]0}, \[Theta][t], {t, 0,
100}];
Rol[a_, b_, \[Mu]_, \[Sigma]_, W_, A_, c_] := Module[{w, Flag, n, X},
X = Sol[a, b, \[Mu], \[Sigma], W, A, c];
w = (0.9) ((D[(\[Theta][t] /. X),
t]) /. (NSolve[(\[Theta][t] /. X) == \[Pi]/2, t]))[[1]][[1]];
Flag = NumericQ[w] && Positive[w];
n = c;
While[Flag,
n++;
X = Sol[a, b, \[Mu], \[Sigma], w, 0, n];
w = (((D[(\[Theta][t] /. X),
t]) /. (NSolve[(\[Theta][t] /. X) == \[Pi]/2, t]))[[1]][[
1]]) (0.9);
Flag = NumericQ[w] && Positive[w];
];
n
]
Fr[a_, b_, \[Mu]_, \[Sigma]_, \[Theta]0_, c_] :=
Table[Mod[Rol[a, b, \[Mu], \[Sigma], w, \[Theta]0, c], 4], {w, 0, 1,
0.01}]
Pr[a_, b_, \[Mu]_, \[Sigma]_, \[Theta]0_, c_] := Module[{temp},
temp = Fr[a, b, \[Mu], \[Sigma], \[Theta]0, c];
{Count[temp, 0]/Length[temp], Count[temp, 1]/Length[temp],
Count[temp, 2]/Length[temp], Count[temp, 3]/Length[temp]}
]
Pro[a_, b_, \[Mu]_, \[Sigma]_, c_] := Module[{temp},
temp = Table[
Pr[a, b, \[Mu], \[Sigma], \[Pi]/2 (i), c], {i, 0, 1, 0.1}];
Total[temp]/Length[temp]
]
M0[a_, b_, \[Mu]_, \[Sigma]_] = (\[Mu] b + \[Sigma]/2)/(
2 \[Mu] + 2 \[Sigma]);
M1[a_, b_, \[Mu]_, \[Sigma]_] = (\[Mu] a + \[Sigma]/2)/(
2 \[Mu] + 2 \[Sigma]);
M2[a_, b_, \[Mu]_, \[Sigma]_] = (\[Mu] (1 - b) + \[Sigma]/2)/(
2 \[Mu] + 2 \[Sigma]);
M3[a_, b_, \[Mu]_, \[Sigma]_] = (\[Mu] (1 - a) + \[Sigma]/2)/(
2 \[Mu] + 2 \[Sigma]);
C0[a_, b_, \[Mu]_, \[Sigma]_] = (
M0[a, b, \[Mu], \[Sigma]] + M1[a, b, \[Mu], \[Sigma]])/2;
C1[a_, b_, \[Mu]_, \[Sigma]_] = (
M1[a, b, \[Mu], \[Sigma]] + M2[a, b, \[Mu], \[Sigma]])/2;
C2[a_, b_, \[Mu]_, \[Sigma]_] = (
M2[a, b, \[Mu], \[Sigma]] + M3[a, b, \[Mu], \[Sigma]])/2;
C3[a_, b_, \[Mu]_, \[Sigma]_] = (
M3[a, b, \[Mu], \[Sigma]] + M0[a, b, \[Mu], \[Sigma]])/2;
Prob[a_, b_, \[Mu]_, \[Sigma]_] := (1/
2) (C0[a, b, \[Mu], \[Sigma]] Pro[a, b, \[Mu], \[Sigma], 0] +
C1[a, b, \[Mu], \[Sigma]] Pro[a, b, \[Mu], \[Sigma], 1] +
C2[a, b, \[Mu], \[Sigma]] Pro[a, b, \[Mu], \[Sigma], 2] +
C3[a, b, \[Mu], \[Sigma]] Pro[a, b, \[Mu], \[Sigma], 3]) + (1/
2) (C0[1 - a, b, \[Mu], \[Sigma]] Pro[1 - a, b, \[Mu], \[Sigma],
0] + C1[1 - a, b, \[Mu], \[Sigma]] Pro[1 - a,
b, \[Mu], \[Sigma], 1] +
C2[1 - a, b, \[Mu], \[Sigma]] Pro[1 - a, b, \[Mu], \[Sigma],
2] + C3[1 - a, b, \[Mu], \[Sigma]] Pro[1 - a,
b, \[Mu], \[Sigma], 3])[[{1, 4, 3, 2}]]
Prob[a, b, \[Mu], \[Sigma]]