Let the ellipse be in standard position, with (fraction-free) equation
$$b^2 x^2 + a^2 y^2 = a^2 b^2$$
Let our equilateral triangle have circumcenter $(p,q)$ and circumradius $r$. Note that maximizing the area of the triangle is equivalent to maximizing $r$.
For some angle $\theta$ ---actually, for three choices of $\theta$--- the vertices of the triangle have coordinates
$$
(p,q) + r \; \mathrm{cis}\theta \qquad (p,q)+r\;\mathrm{cis}\left(\theta+120^{\circ}\right) \qquad (p,q) + r\;\mathrm{cis}\left(\theta-120^{\circ}\right)
$$
where I abuse the notation "$\mathrm{cis}\theta$" to indicate the vector $(\cos\theta,\sin\theta)$.
Substituting these coordinates into the ellipse equation gives a system of three equations in four parameters $p$, $q$, $r$, $\theta$. I used Mathematica's Resultant[] function to help me eliminate $r$ and $\theta$, arriving at a huge polynomial equation in $p$ and $q$. One factor of the polynomial gives rise to this equation:
$$p^2 b^2\left( a^2+3b^2 \right)^2 + q^2 a^2\left(3a^2+b^2\right)^2= a^2b^2\left(a^2-b^2\right)^2$$
This says that the family of circumcenters $(p,q)$ lie on their own ellipse! We can therefore write
$$p = \frac{a\left(a^2-b^2\right)}{a^2+3b^2}\cos\phi \qquad q = \frac{b\left(a^2-b^2\right)}{3a^2+b^2}\sin\phi$$
for some $\phi$. Back-substituting into the system of equations gives this formula for $r$:
$$r = \frac{4 a b \sqrt{a^2\left(a^2+3b^2\right)^2-\left(a-b\right)^3\left(a+b\right)^3 \cos^2\phi}}{\left(a^2 + 3 b^2\right)\left(3a^2+b^2\right)}$$
The maximum value, $R$, is attained when $\cos\phi = 0$, so
$$R := \frac{4a^2b}{3a^2+b^2}$$
The area of the maximal triangle is
$$\frac{3\sqrt{3}}{4}R^2 = \frac{12a^4 b^2\sqrt{3}}{\left(3a^2+b^2\right)^2}$$
Perhaps-unsurprisingly, the corresponding triangles are centered horizontally within the ellipse, with a vertex at either the top or bottom of the minor axis.
Now, I should point out that my big $pq$ polynomial has other factors, namely, $p$ itself, $q$ itself, and a giant I'll call $f$.
One can verify that the cases $p=0$ and $q=0$ lead to the same results as above. (Specifically, they correspond to the respective cases $\cos\phi=0$ and $\sin\phi=0$.) Intuitively, if a circle's center lies on an axis of the ellipse, then the points of intersection with the ellipse have reflective symmetry over that axis. If there are four distinct points (or two, or none), then we cannot choose three to be the vertices of our equilateral triangle; consequently, there must be only three points of intersection, with one of them on the axis, serving as the point of tangency for the circle and ellipse.
As for the case $f=0$ ... I'll just irresponsibly call it extraneous. (The method of resultants tends to spawn such things.)
In the standard notation we need to prove that
$$\frac{a^2+b^2-c^2}{2ab}<\cos60^{\circ}$$ or
$$c^2>a^2-ab+b^2$$ or
$$\sqrt[3]{(a^3+b^3)^2}>a^2-ab+b^2$$ or
$$(a+b)^2(a^2-ab+b^2)^2>(a^2-ab+b^2)^3$$ or
$$ab>0,$$ which is true.
Id est, $$\measuredangle ACB>60^{\circ}$$ and we are done!
Best Answer
Reflect $M$ with respect to the sides of $ABC$. You get an hexagon whose area is $2S$:
The hexagon can be decomposed as the union between $A'B'C'$ (whose side lengths are $\sqrt{3}MA,\sqrt{3}MB,\sqrt{3}MC$) and the isosceles triangles $CA'B',BC'A',AB'C'$. It follows that
$$ 2S = \frac{\sqrt{3}}{4}(AM^2+BM^2+CM^2)+ 3 S'$$
where $S'$ is the area of a triangle with side lengths $MA,MB,MC$. By Weitzenbock's inequality $AM^2+BM^2+CM^2 \geq 4\sqrt{3}S'$, hence $S'\leq \frac{S}{3}$ as wanted.