It is not clear to me why this optimization is expected to recover phi
. The discrepancy in your results is actually user error: in going from the solution of this
Solve[D[llog[\[Phi]], \[Phi]] == 0, \[Phi]]
to this
f2[{m0_, m1_, m2_}] := 2 ArcTan[Sqrt[m1 + 2 m2]/Sqrt[2 m0 + m1]]
you changed 4s to 2s. This is equivalent to redefining
L[\[Phi]_, m0_, m1_, m2_] := (Cos[\[Phi]/2])^(2 m0)*(Sin[\[Phi]]/2)^
m1*(Sin[\[Phi]/2])^(2 m2)
Note the first and third exponents have factors of 2 rather than 4. With this change, \[Phi] /. NMaximize[{d[[i]], 0 <= \[Phi] <= Pi}, \[Phi]][[2]]
will average pretty close to 1.5. (I still don't know why that would recover the 1.5, but will take on faith that it can do that. Must be a square root is useful somehow.)