A challenging inequality that involves the digamma function and polygamma functions Let $f(x)=x \psi(x+1)$, where $\psi$ is the digamma function. Define 
$$g(x)=(f(ax)+f((1-a)x)-f(x))-(f(ax+by)+f((1-a)x+(1-b)y)-f(x+y)),$$
where $0\le a,b\le 1$ and $x,y\ge 0$.
How to show that $g(x)$ is increasing in $x$ on $[0,\infty)$? Thank you!

Here is what I have tried. Let $h(x,y)=f(ax+by)+f((1-a)x+(1-b)y)-f(x+y)$. Then we have $g(x)=h(x,0)-h(x,y)$. To show that $g(x)$ is increasing, it suffices to show that $$\frac{\partial}{\partial x}(h(x,0)-h(x,y))\ge 0. $$
Then it is sufficient to show that $$
\frac{\partial^2 h(x,y)}{\partial y \partial x}\le 0,
$$
which is equivalent to $$
f''(x+y)\ge abf''(ax+by)+(1-a)(1-b)f''((1-a)x+(1-b)y),
$$
where $f''(x)=2\psi_1(1+x)+x\psi_2(1+x)$, where $\psi_1$ and $\psi_2$ are the first and second derivatives of the digamma function $\psi$. Then if I use the series representation of polygamma functions
$$\psi_m(z)=(-1)^{m+1}m!\sum_{k\ge 0}\frac{1}{(z+k)^{m+1}},$$
the above inequality is equivalent to 
$$
\sum_{k\ge 1}k\left[ \frac{1}{(x+y+k)^3}-\frac{ab}{(ax+by+k)^3}-\frac{(1-a)(1-b)}{((1-a)x+(1-b)y+k)^3} \right]\ge 0.
$$
I find that if $a,b,x,y$ are fixed, $\frac{1}{(x+y+k)^3}-\frac{ab}{(ax+by+k)^3}-\frac{(1-a)(1-b)}{((1-a)x+(1-b)y+k)^3} $ is negative for small $k$ and positive when $k$  is sufficiently large. So it seems that the inequality cannot be addressed term by term. I am thinking about if it can be proved by re-arranging the terms in the series and grouping them into positive groups. I also tried the similar strategy by using the integral representation of polygamma functions $$
\psi_m(z)=(-1)^{m+1}\int_0^\infty \frac{t^m e^{-zt}}{1-e^{-t}}dt,
$$
but it did not work either. 
 A: Max reduced the problem to proving that $(xf'')'\geqslant 0$. We have $f=x\psi(x+1)=1+x\psi(x)$, so $f''=(x\psi(x))''=x\psi''+2\psi'$, $f'''=x\psi'''+3\psi''$,$(xf'')'=xf'''+f''=x^2\psi'''+4x\psi''+2\psi'$. Exactly this guy is proved to be non-negative in the recent clever MO answer by Iosif Pinelis to a different question. I wonder whether this is coincidence.
A: Numerical evidence suggests that for any $a\in(0,1]$ and $x\geq 0$, we have
$$(\star)\qquad f''(ax)\leq \frac{f''(x)}{a}.$$
The inequality $(\star)$ allows to easily prove that
$$f''(x+y)\ge abf''(ax+by)+(1-a)(1-b)f''((1-a)x+(1-b)y).$$
Indeed, we have $f''(ax+by) = f''(\frac{ax+by}{x+y}(x+y)) \leq \frac{x+y}{ax+by} f''(x+y)$ and similarly $f''((1-a)x+(1-b)y) \leq \frac{x+y}{(1-a)x+(1-b)y}f''(x+y)$.
Since $f''(x+y)>0$, it remains to notice that
$$1 - ab\frac{x+y}{ax+by} - (1-a)(1-b)\frac{x+y}{(1-a)x+(1-b)y} = \frac{(a-b)^2xy}{(ax+by)((1-a)x+(1-b)y)} \geq 0.$$
I do not yet have a proof of $(\star)$.
A: Inspired by Max's answer, I finally have a proof for the inequality $$h(x)\triangleq f''(x)-a f''(ax)\ge 0.$$
Recall that $\psi_k(z)\triangleq \frac{d^k}{dz^k}\psi(z)$.
Using the two-sided approximation presented in https://arxiv.org/abs/0903.1984 $$\frac{(k-1)!}{\Bigl\{x+\Bigl[\frac{(k-1)!}{|\psi_k(1)|}\Bigr]^{1/k}\Bigr\}^k} +\frac{k!}{x^{k+1}}<\bigl|\psi_k(x)\bigr|<\frac{(k-1)!}{\bigl(x+\frac12\bigr)^k}+\frac{k!}{x^{k+1}},$$ we have $$
h(x)\le x \left(a^2 \left(\frac{4}{(2 a x+3)^2}+\frac{2}{(a x+1)^3}\right)-\frac{2}{(x+1)^3}-\frac{1}{\left(x+1+\sqrt{-\frac{1}{\psi ^{(2)}(1)}}\right)^2}\right)+2 \left(-\frac{a}{(a x+1)^2}-\frac{a}{a x+\frac{6}{\pi ^2}+1}+\frac{1}{(x+1)^2}+\frac{1}{x+\frac{3}{2}}\right)
$$ and $$
h(x) \ge x \left(a^2 \left(\frac{2}{(a x+1)^3}+\frac{1}{\left(a x+1+\sqrt{-\frac{1}{\psi ^{(2)}(1)}}\right)^2}\right)-\frac{2}{(x+1)^3}-\frac{4}{(2 x+3)^2}\right)+2 \left(-\frac{a}{(a x+1)^2}-\frac{2 a}{2 a x+3}+\frac{1}{(x+1)^2}+\frac{1}{x+\frac{6}{\pi ^2}+1}\right).
$$
By the squeeze theorem, we know that $$
\lim_{x\to \infty} h(x)=0.
$$
Now it remains to show that $h(x)$ is decreasing in $x$. Let us show an even stronger result that $h(x)$ is completely monotone. By Bernstein's theorem on monotone functions, it suffices to show that $$h(x)=\int_0^\infty e^{-tx} dg(t),$$
where $g(x)$ is the cdf of a non-negative finite Borel measure on $[0,\infty)$. 
Note that $$h(x)=\int_0^\infty\frac{\left(t^2 \left(a^2 \text{csch}^2\left(\frac{t}{2}\right)-\text{csch}^2\left(\frac{t}{2 a}\right)\right)\right) \exp (-t x)}{4 a^2}dt,$$ which completes the proof. This integral representation is obtained by using the integral representation of polygamma functions $$\psi_m(z)=(-1)^{m+1}\int_0^\infty \frac{t^m e^{-zt}}{1-e^{-t}}dt.$$ Plugging in the above integral representation and change the variable $at$ to $t$, we have $$h(x)= \int_0^{\infty } t (2-t x) \left(\frac{1}{\exp (t)-1}-\frac{1}{a \left(\exp \left(\frac{t}{a}\right)-1\right)}\right) \exp (-t x) \, dt. $$
Then we split it into two parts $h(x)=h_1(x)+h_2(x)$, where $$h_1(x)=\int_0^{\infty } 2t  \left(\frac{1}{\exp (t)-1}-\frac{1}{a \left(\exp \left(\frac{t}{a}\right)-1\right)}\right) \exp (-t x) \, dt$$ and $$h_2(x)=\int_0^{\infty } t (-t x) \left(\frac{1}{\exp (t)-1}-\frac{1}{a \left(\exp \left(\frac{t}{a}\right)-1\right)}\right) \exp (-t x) \, dt\\=\int_0^{\infty } t^2  \left(\frac{1}{\exp (t)-1}-\frac{1}{a \left(\exp \left(\frac{t}{a}\right)-1\right)}\right) d\exp (-t x)\\ =\int_0^\infty \exp (-t x) d\left(t^2  \left(\frac{1}{\exp (t)-1}-\frac{1}{a \left(\exp \left(\frac{t}{a}\right)-1\right)}\right)\right)\\=\int_0^\infty \exp(-tx) t \left(\frac{t}{a^2 \left(e^{t/a}-1\right)^2}+\frac{t-2 a}{a^2 \left(e^{t/a}-1\right)}-\frac{e^t (t-2)+2}{\left(e^t-1\right)^2}\right) dt\,.$$
Putting them together again and simplifying it yields the desired result.
