This is a question of the symmetric type, such as listed in:
With a constraint $\;x+y+z=1\;$ and $\;x,y,z > 0$ . Sort of a general method to transform such a constraint into the inside of a triangle in 2-D has been explained at length in:
Our function $f$ in this case is:
$$ f(x,y,z) = \dfrac{x+y}{\sqrt{x^2+xy+y^2+yz}}+\dfrac{y+z}{\sqrt{y^2+yz+z^2+zx}}+\dfrac{z+x}{\sqrt{z^2+zx+x^2+xy}} - \left[2+\sqrt{\dfrac{xy+yz+xz}{x^2+y^2+z^2}}\right]$$

Places where $\;|f(x,y,z)|\;$ is less than, say, $\;\epsilon\;$ ("eps", see below), are colored $\color{blue}{blue}$.
The blue spots in the pictures are a proof without words
that the only minimum is indeed close to zero and near the midpoint $M$ of the equilateral triangle: $(x_M,y_M,z_M) = (1/3,1/3,1/3)$ .
- Picture at top left: geometry of the conditions $\;x,y,z > 0\;$ and $\;x+y+z=1$ .
- Picture at top right: contour lines of $f(x,y,z)$ , as seen in the plane of the $\color{red}{red}$ triangle,
are at $25$ equidistant levels, between the minimum and the maximum as found within the viewport
- Picture at bottom left: enlargement of the $\color{blue}{blue\, spot}$ in the picture at top right
- Picture at bottom right: same enlargement but with the value of $\epsilon$ reduced to a much smaller value
A few further details:
- Function values close to viewport minimum are white
- Function values close to viewport maximum are black
- Contour grey values are just the other way around
A = (1,0,0); O = (0,0,0); eps := 0.0003;
B = (0,1,0); M = (1/3,1/3,1/3); xmin := -0.7; xmax := 0.7;
C = (0,0,1); ymin := -0.6; ymax := 0.8;
xmin := -0.02; xmax := 0.02; eps := eps/4096;
ymin := -0.02; ymax := 0.02;
Two or three iterations result in pictures that, quite remarkably,
can hardly be distinguished from the two pictures at the bottom:
- Picture at bottom left: enlargement $\times 64$ of the $\color{blue}{blue\,spot}$ in the picture at bottom right
- Picture at bottom right: same enlargement but with the value of $\;\epsilon\;$ reduced to $\;\epsilon/4096\;$
xmin := xmin/64; xmax := xmax/64; eps := eps/4096;
ymin := ymin/64; ymax := ymax/64;
With calculations ending in a "range check error" because $\epsilon$ ( = eps) has become zero (too) quickly. Which - alright? - is a proof of the statement up to no less than machine precision.
A local and planar orthogonal coordinate system for the equilateral triangle can be obtained as:
$$
\vec{e}_\xi = \vec{OB}-\vec{OA} = \left[ \begin{array}{c} -1 \\ +1 \\ 0 \end{array} \right] / \sqrt{2}\\
\vec{e}_\eta = \vec{OC}-\left(\vec{OA}+\vec{OB}\right)/2 = \left[ \begin{array}{c} -1/2 \\ -1/2 \\ +1 \end{array} \right] / \sqrt{3/2}
$$
Giving for the relationship between global coordinates $(x,y,z)$ and local triangle coordinates $(\xi,\eta)$:
$$
x = -\frac{1}{\sqrt{2}} \xi -\frac{1}{\sqrt{6}} \eta + \frac{1}{3}\\
y = +\frac{1}{\sqrt{2}} \xi -\frac{1}{\sqrt{6}} \eta + \frac{1}{3}\\
z = \sqrt{\frac{2}{3}} \eta + \frac{1}{3}
$$
It's easily verified that $x+y+z=1$ ; the origin of the local coordinate system is at $(\xi,\eta) = (0,0)$ .
In this way, the function can be expressed in local triangle coordinates, but (at least for me) that didn't simplify the puzzle much.
A common method to study the function $f(x,y,z)$ in the neighborhood of $M$ is to take derivatives at that point (without a computer algebra system it's a tedious job). To begin with, the zero'th order derivative is zero:
$$
f(M) = 0
$$
The first order derivatives are zero as well, meaning that there is an extremum to be expected there:
$$
\left[ \begin{array}{c}
\frac{\partial f}{\partial x} \\ \frac{\partial f}{\partial y} \\ \frac{\partial f}{\partial z}
\end{array} \right] (M) =
\left[ \begin{array}{c} 0 \\ 0 \\ 0 \end{array} \right]
$$
The second order derivatives are summarized as a Hessian matrix:
$$
\left[ \begin{array}{ccc}
\frac{\partial^2 f}{\partial x^2} & \frac{\partial^ f}{\partial x \partial y} & \frac{\partial^2 f}{\partial x \partial z} \\
\frac{\partial^2 f}{\partial y \partial x} & \frac{\partial^ f}{\partial y^2} & \frac{\partial^2 f}{\partial y \partial z} \\
\frac{\partial^2 f}{\partial z \partial x} & \frac{\partial^ f}{\partial z \partial y} & \frac{\partial^2 f}{\partial z^2}
\end{array} \right] (M) =
\frac{51}{64} \left[ \begin{array}{ccc} 2 & -1 & -1 \\ -1 & 2 & -1 \\ -1 & -1 & 2 \end{array} \right]
$$
Therefore $f(x,y,z)(M)$ can be approximated by the following Taylor series expansion:
$$
f(x,y,z)(M) \approx \frac{51}{32} \left( x^2 + y^2 + z^2 - x y - x z - y z \right)
$$
To investigate the nature of this expression,
eigenvalues and eigenvectors of the Hessian matrix are determined:
$$
\left| \begin{array}{ccc} 2-\lambda & -1 & -1 \\ -1 & 2-\lambda & -1 \\ -1 & -1 & 2-\lambda \end{array} \right| =
-\lambda^3 + 6 \lambda^2 - 9\lambda = 0 \qquad \Longrightarrow \qquad \lambda \in \left\{ 0,3 \right\}
$$
The normed eigenvector corresponding with $\;\lambda=0\;$ is $\;(1,1,1)/\sqrt{3}\;$ and the normed eigenvectors corresponding
with $\;\lambda=3\;$ is anything in the plane $\;x+y+z=0\;$ , which is parallel to the plane of our equilateral triangle. So we can choose the local triangle coordinate system for the latter:
$$
\vec{e}_\xi = \left[ \begin{array}{c} -1 \\ +1 \\ 0 \end{array} \right] / \sqrt{2} \qquad
\vec{e}_\eta = \left[ \begin{array}{c} -1/2 \\ -1/2 \\ +1 \end{array} \right] / \sqrt{3/2} \qquad
\vec{e}_\zeta = \left[ \begin{array}{c} +1 \\ +1 \\ +1 \end{array} \right] / \sqrt{3}
$$
If we express the "conic section" into those local coordinates, then what we get is the following:
$$
x^2 + y^2 + z^2 - x y - x z - y z = 3 \xi^2 + 3 \eta^2 = \frac{\xi^2}{1/3} + \frac{\eta^2}{1/3} + \frac{\zeta^2}{\infty}
$$
The iso-surfaces are a circular cylinder with axis in the direction of $\;\vec{OM} = (1,1,1)/3 $ . Thus, as a first approximation,
the isolines in the neighborhood of $M$ must be circles, which is exactly what's observed in the above pictures:
$$
f(x,y,z)(M) \approx \frac{153}{32} \left( \xi^2 + \eta^2 \right) \qquad \mbox{for} \qquad x+y+z=1
$$
Behaviour of the abovementioned iterations is explained quite well with this theory.
It is also trivial that the only minimum found with this approximation is $\,0\,$ at the midpoint $M$ with $(\xi,\eta) = (0,0)$ .
Guess something more is needed, though, to convert the last argument into a die hard proof.