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.)