What values of $x,y,z$ maximize $f(x,y,z,p) = x^4+y^4+z^4 + p(xy+xz+yz)$ with constant $p \geq 0$ with the constraint $x^2+y^2+z^2=1$?
Some preliminary studies in Mathematica showed that the behavior of the solution changes around $p=1$. For $p>1$, the problem seems solved by $x=y=z=\frac{1}{\sqrt 3}$. For $0 \leq p \leq 1$ it appears the solution (up to permutations of x,y,z) looks like x > y = z. These solutions appear to favor one or the other of the two terms in $f$, and I find the $p>1$ behavior striking and would appreciate a proof that $x=y=z=\frac{1}{\sqrt 3}$ indeed maximizes.
UPDATE
For at least $p>4/3$, one can show that $x=y=z=\frac{1}{\sqrt 3}$ maximizes f.
The proof for $1<p<4/3$ evades me currently (the left bound, as mentioned in comments, may be slightly higher than 1).
I'll assume without loss of generality that $x,y,z$ are all positive. Note that $f(\frac{1}{\sqrt 3},\frac{1}{\sqrt 3},\frac{1}{\sqrt 3}, p) = p+\frac{1}{3} $, so if we show that for some $p$ we enjoy $f(x,y,z,p) \leq p+\frac{1}{3}$, we've shown $x=y=z=\frac{1}{\sqrt 3}$ maximizes $f$.
Now consider that via the constraint, $$f(x,y,z,p) = 1-2(x^2y^2+x^2z^2+y^2z^2) + p(xy+xz+yz)$$ The root mean square compared to arithmetic mean inequality tells us that for nonnegative $a,b,c$ we enjoy $\sqrt{\frac{a^2+b^2+c^2}{3}} \geq \frac{a+b+c}{3},$ so $-(a^2+b^2+c^2) \leq -\frac{(a+b+c)^2}{3}$. Then $$f(x,y,z,p) \leq 1 + p(xy+xz+yz) -\frac{2}{3}(xy+xz+yz)^2.$$ Now, let's call $(xy+xz+yz)=q$. Note that it is easy to prove with our constraints that $0\leq q \leq 1$. Then we can bound the maximum of $f(x,y,z,p)$ by finding the maximum of $g(q) = 1+pq-\frac{2}{3}q^2$ on $0\leq q \leq 1$. The derivative $g'(q) = p-\frac{4}{3}q$, so $g(q)$ increases up to $q=\frac{3}{4}p$. Then on $0\leq q \leq 1$, we have $q_{max} = min(\frac{3}{4}p, 1)$. Thus for $p \geq \frac{4}{3}$ we have $q_{max}=1$, so $$ f(x,y,z,p) \leq 1+p*1-\frac{2}{3}*1^2 = p + 1/3$$ Thus since $f(\frac{1}{\sqrt 3},\frac{1}{\sqrt 3},\frac{1}{\sqrt 3}, p) = p+\frac{1}{3}$, we have for $p \geq \frac{4}{3}$ that $x=y=z=\frac{1}{\sqrt 3}$ maximizes $f$.
I have attempted the method of Lagrange multipliers, which led me to three cubic equations:
I set out to maximize $4(x^3+y^3+z^3)+p(xy+xz+yz)+\lambda(x^2+y^2+z^2-1)$.
This yielded $$4x^3+2x\lambda + p(y+z)=0$$ $$4y^3+2y\lambda + p(x+z)=0$$ $$4z^3+2z\lambda + p(y+x)=0$$ But now I feel a little at a loss on how to solve this system of equations. One can verify that $x=y=z=\frac{1}{\sqrt{3}}$ solves these equations at any p with the right choice of lamda. I'm unsure about pairing Hessian techniques to verify the local nature of the critical point when using lagrange multipliers.
I have also tried incorporating the constraint by writing $z=\cos(\theta)$, $x=\sin(\theta)\cos(\phi)$ and $y=\sin(\theta)\sin(\phi)$.
This yields $$f(\theta, \phi, p) = \cos(\theta)^4+(\sin(\theta)\cos(\phi))^4+(\sin(\theta)\sin(\phi))^4+p(\cos(\theta)\sin(\theta)\sin(\phi)+\cos(\theta)\sin(\theta)\cos(\phi)+(\sin(\theta))^2\sin(\phi)\cos(\phi) ) $$ which I am unsure how to maximize. However, with motivation from Michael Behrend, I've found that studying $\phi=\frac{\pi}{4}$ shows interesting critical point behavior slightly above p=1 where three critical points undergo some sort of bifurcation leaving just a single critical point at the eternal $x=y=z=\frac{1}{3}$ .
Lagrange function is $$L=x^4+y^4+z^4+p(xy+xz+yz)+\lambda(x^2+y^2+z^2-1).$$ Solving system $$L'_x=0,\;L'_y=0,\;L'_z=0,\;L'_x=0,\;L'_\lambda=0$$ we get that if $p\ge1$ global maximum is $$f_{max}=p+\frac13$$ Maximum points is:
if $p>1$ then $[x=\frac{1}{\sqrt{3}},y=\frac{1}{\sqrt{3}},z=\frac{1}{\sqrt{3}}]$ or $[x=-\frac{1}{\sqrt{3}},y=-\frac{1}{\sqrt{3}},z=-\frac{1}{\sqrt{3}}]$
if $p=1$ then $[x=\frac{\sqrt{6}}{3},y=\frac{1}{\sqrt{6}},z=\frac{1}{\sqrt{6}}],[x=-\frac{\sqrt{6}}{3},y=-\frac{1}{\sqrt{6}},z=-\frac{1}{\sqrt{6}}],[x=\frac{1}{\sqrt{6}},y=\frac{\sqrt{6}}{3},z=\frac{1}{\sqrt{6}}],[x=\frac{1}{\sqrt{6}},y=\frac{1}{\sqrt{6}},z=\frac{\sqrt{6}}{3}],[x=-\frac{1}{\sqrt{6}},y=-\frac{\sqrt{6}}{3},z=-\frac{1}{\sqrt{6}}],[x=-\frac{1}{\sqrt{6}},y=-\frac{1}{\sqrt{6}},z=-\frac{\sqrt{6}}{3}],[x=\frac{1}{\sqrt{3}},y=\frac{1}{\sqrt{3}},z=\frac{1}{\sqrt{3}}],[x=-\frac{1}{\sqrt{3}},y=-\frac{1}{\sqrt{3}},z=-\frac{1}{\sqrt{3}}]$
If $p=0$ then $$f_{max}=1,$$ $[x=\pm1,y=0,z=0],[x=0,y=\pm1,z=0],[x=0,y=0,z=\pm1]$.
If $0<p<1$ we can't exact solve system.