How to prove that $a^2b+b^2c+c^2a \leqslant 3$, where $a,b,c >0$, and $a^ab^bc^c=1$

$a,b,c >0$, and $a^ab^bc^c=1$, prove $$a^2b+b^2c+c^2a \leqslant 3$$

I don't even know what to do with the condition $a^ab^bc^c=1$. At first I think $x^x>1$, but I was wrong. This inequality is true, following by the verification from Mathematica

We employ of the rearrangement inequality.
First, since $x\mapsto x^2$ preserves the order for $x>0,$ we have $$a^2b+b^2c+c^2a\le a^3+b^3+c^3.$$
Next write the condition as $$a\ln a+b\ln b+c\ln c=0.$$ Again using the rearrangement inequality, as $x\mapsto\ln x$ preserves the order for $x>0,$ we have the inequalities: $$\begin{cases}a\ln b+b\ln c+c\ln a\le0\\ a\ln c+b\ln a+c\ln b\le0\\ a\ln a+b\ln b+c\ln c=0\end{cases}$$
Adding these together, we have $(a+b+c)\ln(abc)\le0,$ hence $abc\le1.$
Now let $\ln a+\ln b+\ln c=k\le0.$ Apply the Lagrange multiplier method with condition $g(a,b,c):=\ln a+\ln b+\ln c=k$ for a fixed $k\le0,$ to maximize $f(a,b,c):=a^3+b^3+c^3.$
Thus the Lagrange multiplier gives that the extreme of $f$ occurs when $$\begin{cases}3a^2-\lambda\frac{1}{a}=0\\ 3b^2-\lambda\frac{1}{b}=0\\ 3c^2-\lambda\frac{1}{c}=0\end{cases},$$ i.e. when $a^3=b^3=c^3=\lambda;$ then $a=b=c$ and hence $f(a,b,c)=3abc=3e^k\le3.$
$\square$

Hope this helps.

Edit:
As pointed out in the comment, the final part about Lagrange multipliers is incorrect; in fact, given the constraint $abc<1,$ it does not follow that $a^2b+b^2c+c^2a\le3.$
Also pointed out in the comment is to use Jensen's inequality to obtain $a+b+c\le3,$ but we are not getting answers yet.
Everything we tried so far to fix this fails. We shall update if we find a way to work around.

We shall prove the following inequalities:

If $x+y+z=3$ and $x,y,z>0$, then $$\tag{1}x^2y+y^2z+z^2x+xyz\le 4\label{1},$$ $$\tag{2}3x^xy^yz^z+xyz\ge 4\label{2},$$ i.e., $$\tag{*}x^2y+y^2z+z^2x\le 3x^xy^yz^z\label{*}.$$

Accept $\eqref{*}$ for a moment. Homogenizing the inequality by substitution

$$x=\frac{3a}{a+b+c},~y=\frac{3b}{a+b+c},~z=\frac{3c}{a+b+c},$$

we have $\forall a,b,c>0$,

$$\tag{**}a^2b+b^2c+c^2a\le 3(a^ab^bc^c)^{3/(a+b+c)}.\label{**}$$

If $a^ab^bc^c=1$, then we get the original inequality as desired.

Proof of $\eqref{1}$:

Note the LHS is cyclic, we cannot assume a specific order like $x\le y\le z$, but we can assume WLOG that $y$ is in the middle (neither the minimum nor the maximum), then

$$z(y-x)(y-z)\le 0\\ \Rightarrow x^2y+y^2z+z^2x+xyz\le y(x^2+2xz+z^2)=\frac{1}{2}(2y)(x+z)^2\le\frac{1}{2}\left(\frac{2(x+y+z)}{3}\right)^3=4,$$ by AM-GM. Equality holds when $(x,y,z)=(1,1,1),(2,1,0)$ along with the cyclic permutations.

Proof of $\eqref{2}$:

WLOG we assume $x=\min(x,y,z)$ and consider the following two cases:

• $x\ge 1/3$, then $y,z\ge 1/3$. We note that the function

$$f(t)=\left(t+\frac{1}{3}\right)\ln t$$

is convex for $t\ge 1/3$, as

$$f''(t)=\frac{3t-1}{3t^2}\ge 0.$$

By Jensen,

$$ \begin{align} &~~~~~~~x^{x+1/3}y^{y+1/3}z^{z+1/3}=\exp\left\{f(x)+f(y)+f(z)\right\} \ge\exp\left\{3f\left(\frac{x+y+z}{3}\right)\right\}=\exp\left\{3f(1)\right\}=1,\\ &\Rightarrow x^xy^yz^z\ge (xyz)^{-1/3},\\ &\Rightarrow 3x^xy^yz^z+xyz\ge 3(xyz)^{-1/3}+xyz\ge 4\sqrt[4]{(xyz)^{-1}(xyz)}=4, \end{align} $$

where we have applied AM-GM for the last line. Equality holds when $x=y=z=1$.

• $0<x\le1/3$. We note that the function

$$g(t)=t\ln t$$

is convex for $t>0$, as

$$g''(t)=\frac{1}{t}> 0.$$

Again by Jensen, we have

$$ \begin{align} &~~~~~~~y^yz^z=\exp\left\{g(y)+g(z)\right\}\ge\exp\left\{2g\left(\frac{y+z}{2}\right)\right\} =\exp\left\{2g\left(\frac{3-x}{2}\right)\right\},\\ &\Rightarrow x^xy^yz^z\ge\exp\left\{g(x)+2g\left(\frac{3-x}{2}\right)\right\} \ge\exp\left\{g\left(\frac{1}{3}\right)+2g\left(\frac{4}{3}\right)\right\} =\frac{32}{27}\sqrt[3]{2}>\frac{4}{3},\\ &\Rightarrow 3x^xy^yz^z+xyz>4+xyz>4. \end{align} $$

The second line follows from that $h(t):=g(t)+2g((3-t)/2)$ is monotonically decreasing for $t\in(0,1)$, so $h(x)\ge h(1/3)$ as $x\le 1/3$. Proof of monotonicity:

$$h'(t)=\ln\left(\frac{2t}{3-t}\right)<0\iff 0<t<1.$$

Post-mortem:

1. Despite the flow of the proof, the starting point is the homogeneous version \eqref{**};

2. \eqref{1} is due to Vasile, which allows us to strengthen the inequality into a symmetric one as \eqref{2} (I am not able to find the originial post on artofproblemsolving, but see e.g., here). It seems that converting \eqref{*} into other symmetric form such as $x^3+y^3+z^3$ is too strong for it to hold;

3. For \eqref{2}, the case of $x\le 1/3$ is unfortunately necessary, as $x^xy^yz^z\ge(xyz)^{-1/3}$ does NOT always hold, e.g., check $x\to0,y-z\to 0$. Is there any better estimate of $x^xy^yz^z$ so we can avoid the derivatives? Another follow-up question is to determine the smallest $k$ such that $kx^xy^yz^z+xyz\ge k+1$ holds given the constraints. It definitely fails for $k\le 1/2$.

Brute Force ($200\times200\times200$ grid) Not a proof, but couldn't resist.
Of course, everybody knows where the maximum is, for reasons of symmetry:
with $\,(a,b,c) = (1,1,1)\,$ we have $\,a^ab^bc^c=1\,$ and $\,a^2b+b^2c+c^2a = 3$ .
The following (Delphi Pascal) program is supposed to be self documenting.

program HN_NH;
{
  Brute Force with Seven Point Stars
  ==================================
}
function pow(x,r : double) : double;
{
  x^r
}
begin
  pow := exp(r*ln(x));
end;

procedure test(veel : integer);
var
  i,j,k,ken : integer;
  a,b,c,d,f,min,max : double;

  procedure vertex(x,y,z : integer);
  var
    a,b,c,h : double;
  begin
    a := (2*i+x)*d; b := (2*j+y)*d; c := (2*k+z)*d;
    h := pow(a,a)*pow(b,b)*pow(c,c);
    if h < 1 then ken := ken*2 else ken := ken*2+1;
  end;

begin
{ Verify maximum (a,b,c)-value < 1.6 }
  Writeln(exp(2/exp(1)),' <',pow(1.6,1.6));
  d := 1.6/veel/2; { half voxel size }
  min := 3; max := 0; { initialize }
  for i := 1 to veel-1 do
  begin
    for j := 1 to veel-1 do
    begin
      for k := 1 to veel-1 do
      begin
        ken := 0; { Binary number for collecting <> }
      { Each vertex of a 7-point star }
        vertex(-1,0,0); vertex(+1,0,0);
        vertex(0,-1,0); vertex(0,+1,0);
        vertex(0,0,-1); vertex(0,0,+1);
        if (ken = 0) or (ken = 63) then Continue;
      { Midpoint of star is near a^a*b^b*c^c = 1 }
        a := 2*i*d; b := 2*j*d ; c := 2*k*d;
        f := sqr(a)*b + sqr(b)*c + sqr(c)*a;
        if f < min then min := f;
      { Determine maximum of f(a,b,c) }
        if f > max then max := f;
      end;
    end;
  end;
  Writeln(min,' < f(a,b,c) <',max);
end;

begin
  test(200);
end.

 2.08706522863453E+0000 < 2.12125057109759E+0000
 9.12537600000000E-0003 < f(a,b,c) < 3.00026265600000E+0000

Note. Explaining the estimate $\{a,b,c\} < \{1.6\}$ in the program: $$ f(x) = x^x = e^{x\ln(x)} \quad \Longrightarrow \quad f'(x) = [1+\ln(x)]e^{x\ln(x)} = 0 \quad \Longrightarrow \quad x=1/e \\ \Longrightarrow \quad f(1/e) = e^{-1/e} $$ This means that the maximum $x$ of each one of the coordinates in the product $a^ab^bc^c=1$ is: $$x^x = e^{2/e} < (1.6)^{1.6}$$ Therefore each of the coordinates $\;x < 1.6\,$ (or $\,1.58892154635044$ , to be double precise).