The alternative forms of the given inequality are
$$x\ln x + y \ln y \ge (x+y)\ln(x+y)+(x+y)\ln\left(1-\dfrac{\sqrt{xy}}{x+y}\right),$$
$$x\ln\,\dfrac x{x+y} + y \ln\, \dfrac y{x+y} \ge (x+y)\ln\left(1-\dfrac{\sqrt{xy}}{x+y}\right),$$
$$\dfrac x{x+y}\,\ln\,\dfrac x{x+y} + \dfrac y{x+y}\, \ln\, \dfrac y{x+y}
\ge \ln\left(1-\sqrt{\dfrac x{x+y}}\;\sqrt{\dfrac y{x+y}}\right).\tag1$$
Let WLOG $\;x\le y,\;$
$$\dfrac yx = \dfrac{1+z}{1-z},\quad \dfrac x{x+y} =\dfrac{1-z}{2}
,\quad \dfrac y{x+y} =\dfrac{1+z}{2},\quad z\in[0,1],\tag2$$
then $(1)$ can be presented in the forms of
$$\dfrac{1-z}{2}\,\ln\dfrac{1-z}{2}+\dfrac{1+z}{2}\,\ln\dfrac{1+z}{2} \ge \ln(1-\,^1\!/_2\sqrt{1-z^2}\,),$$
$$\dfrac{1-z}{2}\,\ln(1-z)+\dfrac{1+z}{2}\,\ln(1+z) \ge \ln(2-\sqrt{1-z^2}),$$
$$-\dfrac{z}{2}\,\ln(1-z)+\dfrac{z}{2}\,\ln(1+z) \ge \ln(2-\sqrt{1-z^2})
-\ln\sqrt{1-z^2},$$
$$f(z) = (1-z)^{-\large\,^z\!/_2}(1+z)^{\large\,^z\!/_2}+1- \dfrac 2{\sqrt{1-z^2}}\ge 0. \tag3$$
Using binomial series in the forms of
$$(1-z)^{-\large\,^z\!/_2} = 1 + \dfrac12 z^2 + \dfrac1{2!}\,\dfrac z2\,\dfrac{2+z}2 z^2
+ \dfrac1{3!}\,\dfrac z2\,\dfrac{2+z}2\dfrac{4+z}2 z^3 +\dots$$
$$= 1+\dfrac12z^2+\dfrac14z^3+\dfrac7{24}z^4+\dfrac14z^5+\dfrac{113}{480}z^6+\dots,$$
$$(1+z)^{\large\,^z\!/_2} = 1 + \dfrac12 z^2 + \dfrac1{2!}\,\dfrac z2\,\dfrac{-2+z}2 z^2
+ \dfrac1{3!}\,\dfrac z2\,\dfrac{-2+z}2\dfrac{-4+z}2 z^3 +\dots$$
$$= 1+\dfrac12z^2-\dfrac14z^3+\dfrac7{24}z^4-\dfrac14z^5+\dfrac{113}{480}z^6+\dots,$$
$$(1-z^2)^{-\large\,^1\!/_2} = 1+\dfrac12z^2+\dfrac38z^4+\dfrac5{16}z^6+\dots,$$
one can get the Maclaurin series in the form of
$$f(z) = \left(1+\dfrac12z^2+\dfrac7{24}z^4+\dfrac{113}{480}z^6+\dots\right)^2
-\left(\dfrac14z^3+\dfrac14z^5+\dots\right)^2$$
$$+1-2\left(1+\dfrac12z^2+\dfrac38z^4+\dfrac5{16}z^6+\dots\right),$$
$$f(z) = \dfrac1{12}z^4+\dfrac3{40}z^6+\dots$$
(see also WA series),
with the all positive coefficients.
Proved!
$\color{brown}{\textbf{About binomial series.}}$
Newton's binomial formulas exist independently of the derivatives associated theory.
If the exponent is the rational, then the coeffitients of the binomial series $(1+x)^{\large\frac pq} = 1+a_1x+a_2x^2+\dots$ are defined via the identity
$$(1+a_1x+a_2x^2+\dots)^q = (1+x)^q,$$
by the comparation of the polynomials in the LHS and RHS.
If the exponent is real, then the binomial series are defined via the limiting transition.
These elementary knowledges point that the binomial series exist independently of the derivatives associated theory.