Let's examine a related problem first: finding a simple, closed curve $C$ from $(1,1)$ to $(-1,-1)$ and back again that maximizes $\displaystyle\oint_C\mathbf{F}\cdot d\mathbf{r}.$ Let $G$ be the region enclosed by $C.$ Using the notation $\mathbf{F}=(F_x, F_y),$ not to be confused with taking the partial derivatives, which I will make explicit, we have by Green's Theorem that
\begin{align*}
\oint_C \mathbf{F}\cdot d\mathbf{r}&=\iint_G\left(\frac{\partial F_y}{\partial x}-\frac{\partial F_x}{\partial y}\right)dA \\
&=\iint_G\left[\left(3-3x^2+2xy\right)-\left(3y^2-3+2xy\right)\right]dx\,dy \\
&=\iint_G\left[6-3x^2-3y^2\right]dx\,dy.
\end{align*}
To maximize this integral, we want to pick up all the region where the integrand is positive. But this is all of $D,$ since
\begin{align*}
6-3x^2-3y^2&\ge 0 \\
6 &\ge 3\left(x^2+y^2\right)\\
2&\ge x^2+y^2.
\end{align*}
This exactly describes $D.$ To pick up all of $D,$ we will just do the circle of radius $\sqrt{2},$ which is the boundary of $D.$ The integral will become
\begin{align*}
\iint_D\left[6-3x^2-3y^2\right]dx\,dy&=\int_0^{2\pi}\int_0^{\sqrt{2}}\left(6-3r^2\right)r\,dr\,d\theta \\
&=6\pi.
\end{align*}
But notice that the final integrand we had was $\theta$-independent. So from the area integral perspective, we wanted to pick up as much of $D$ as possible, and we needed to skirt around it on its boundary. Moreover, it didn't matter which direction we chose: we could have gone clockwise or counter-clockwise.
Now, to solve the original problem, we can see by angle symmetry that we just want to chop the region in half diagonally from $(1,1)$ to $(-1,-1),$ and go around the circle's edge. The maximum value will just be half of what we computed before, or $3\pi.$
This is not a rigorous proof, but I think it could be made into one fairly easily.
