Geometry of Level sets of elliptic polynomials in two real variables Updated:

A polynomial $P(x,y)\in \mathbb{R}[x,y]$ is called an elliptic polynomial if its last homogeneous part does not vanish on $\mathbb{R}^2\setminus\{0\}$.The two answers to this post provide a proof for the following theorem:
Theorem: If $p$ is an elliptic polynomial whose last homogeneous part is positive definitive, then for $c$ sufficiently large
, $p^{-1}(c)$ is a  simple closed curve. Moreover if the centroid of  interior of $p^{-1}(c)$ is denoted by $e_c$ then $e_c$ is convergent as $c$ goes to $+\infty$. The limit $\lim_{c\to \infty} e_c$ can be written in terms of coefficients of $p$. If we drop the ellipticity condition then this convergence result is not necessarily true.

The previous version of the post:
Is there a polynomial function $P:\mathbb{R}^2 \to \mathbb{R}$ with the following property?

For  sufficiently large $c>0$, $P^{-1}(c)$ is a simple closed curve $\gamma_c$, homeomorphic to $S^1$, but as $c$ goes to $+\infty$. the centroid $e_c$ of the interior of $\gamma_c$
does not converge to any point of $\mathbb{R}^2$.

Motivation: The answer is negative if we consider this question for polynomials $p:\mathbb{R} \to \mathbb{R}$ whose eventual level sets are $2$-pointed set, i.e. $S^0$.(Namely a polynomial of even degree). The motivation comes from line -3, item III,  page 4 of Taghavi - On periodic solutions of Liénard equations, which can be generalized to every even degree polynomial with one variable.
 A: Concerning homogeneous polynomials: Let $P(x,y)=\sum_{j=0}^n a_j x^j y^{n-j}$ be such a polynomial, of degree $n$ such that $C:=P^{-1}(\{c\})$ is a simple closed curve for all large enough $c>0$. 
If $n$ is odd, then every line through the origin will have at most one point of intersection with $C$. So, then $C$ cannot be a simple closed curve for any real $c$ -- because every line through any point interior to a simple closed curve must intersect the curve in at least two points. 
It remains to consider the case when $n$ is even. Then $C$ is symmetric about the origin, and hence so is the interior of $C$. Then the centroid of the interior is the origin, and it does not depend on the level $c$. 

Consider now the case of an elliptic polynomial 
\begin{equation*}
 P(z)=P(x,y)=\sum_{j=0}^n a_j x^j y^{n-j}+\sum_{j=0}^{n-1}b_j x^j y^{n-1-j}+K|z|^{n-2} 
\end{equation*}
of (necessarily even) degree $n$, 
where $z:=(x,y)$ and $K=O(1)$ (as $|z|\to\infty$). The ellipticity here is understood as the following condition: 
\begin{equation*}
 \min_{|z|=1}\sum_{j=0}^n a_j x^j y^{n-j}>0. 
\end{equation*}
For any $d_*\in(0,1)$ and any real $D>0$, let $\mathcal P_{n,d_*,D}$ denote the set of all polynomials $p(x)=\sum_{j=0}^n d_j x^j$ such that $d_n\ge d_*$ and $\sum_{j=0}^n|d_j|\le D$. 
Then it is not hard to see that there is a real $c_*(n,d_*,D)>0$, depending only on $n,d_*,D$, such that 
for any polynomial $p(x)=\sum_{j=0}^n d_j x^j$ in $\mathcal P_{n,d_*,D}$ and for all real $c\ge c_*(n,d_*,D)$ the equation $p(x)=c$ has exactly two roots $x_\pm:=x_\pm(c)$ such that $x_-<0<x_+$ and, moreover, 
\begin{equation*}
 x_\pm=\pm\Big(\frac c{d_n}\Big)^{1/n}-(1+o(1))\frac{d_{n-1}}{nd_n} \tag{1}
\end{equation*}
uniformly over all polynomials $p(x)=\sum_{j=0}^n d_j x^j$ in $\mathcal P_{n,d_*,D}$; here and in the sequel the asymptotic relations are for 
$$c\to\infty,$$ 
unless otherwise specified. 
This uniformity can be obtained by refining this reasoning.
Moreover, without loss of generality (wlog), 
\begin{equation*}
 \text{for all $p\in\mathcal P_{n,d_*,D}$ and all real $c\ge c_*(n,d_*,D)$ we have $p'(x_\pm)\ne0$.} \tag{1.5}
\end{equation*}
Indeed, because (1) holds uniformly over all $p\in\mathcal P_{n,d_*,D}$, wlog 
\begin{equation*}
 |x_\pm|\ge\Big(\frac cD\Big)^{1/n}-2\frac D{nd_*}\to\infty, \tag{1.6}
\end{equation*}
so that $|x_\pm|\to\infty$ uniformly over all $p\in\mathcal P_{n,d_*,D}$. 
Also, taking any polynomial $p(x)=\sum_{j=0}^n d_j x^j$ in $\mathcal P_{n,d_*,D}$ and writing $p'(x)=\sum_{j=1}^n d_j jx^{j-1}$, we see that for $|x|\ge1$
\begin{equation*}
 \frac{|p'(x)|}{|x|^{n-1}}\ge nd_n-\sum_{j=1}^{n-1} |d_j| j|x|^{j-n}
 \ge nd_*-n D |x|^{-1}\underset{x\to\infty}\longrightarrow nd_*>0. 
\end{equation*}
So, by (1.6), wlog (1.5) holds indeed.  
Let us now turn back to the elliptic polynomial $P(x,y)$. For each real $t$ consider the polynomial
\begin{equation*}
 p_t(r):=P(r\cos t,r\sin t). 
\end{equation*}
By the ellipticity of the polynomial $P(x,y)$, there exist $d_*\in(0,1)$ and a real $D>0$ such that $p_t\in\mathcal P_{n,d_*,D}$ for all real $t$. Take now any real $c\ge c_*(n,d_*,D)$. Then, by the paragraph right above, for each real $t$ the equation $p_t(r)=c$ has exactly two roots $r_\pm(t):=r_\pm(c;t)$ such that 
$r_-(t)<0<r_+(t)$ and, moreover,
\begin{equation*}
 r_\pm(t)=\pm\Big(\frac c{a(t)}\Big)^{1/n}-(1+o(1))\frac{b(t)}{na(t)} 
\end{equation*}
uniformly in real $t$, 
where 
\begin{equation*}
 a(t):=\sum_{j=0}^n a_j \cos^jt\, \sin^{n-j}t,\quad 
 b(t):=\sum_{j=0}^{n-1} b_j \cos^jt\, \sin^{n-1-j}t.   
\end{equation*}
Moreover, $\frac{dp_t(r)}{dr}|_{r=r_\pm(t)}\ne0$. So, by the implicit function theorem, the functions $r_\pm$ are continuous (in fact, infinitely smooth). 
Also, the functions $r_\pm$ are periodic with period $2\pi$, since for each real $t$ we have $p_{t+2\pi}=p_t$ and the values $r_\pm(t)$ of the the functions $r_\pm$ at $t$ are uniquely determined by the polynomial $p_t$. Furthermore, for all real $r$ and $t$ we have $p_{t+\pi}(r)=p_t(-r)$, which implies $r_+(t+\pi)=-r_-(t)$. So, letting 
$$z_\pm(t):=r_\pm(t)(\cos t,\sin t),$$ 
we see that $z_\pm(t+2\pi)=z_\pm(t)$ and 
$z_+(t)=z_-(t-\pi)$ for all real $t$. So, 
the $c$-level curve of $P(x,y)$ is 
\begin{align*}
C=P^{-1}(\{c\})&=\{z_+(t)\colon t\in\mathbb R\}\cup\{z_-(t)\colon t\in\mathbb R\} \\ 
 &=\{z_+(t)\colon t\in[0,2\pi)\}\cup\{z_-(t)\colon t\in[0,2\pi)\} \\ 
 &=\{z_+(t)\colon t\in[0,2\pi)\} \\ 
&=\{z_+(t)\colon t\in[0,\pi)\}\cup\{z_-(t-\pi)\colon t\in[\pi,2\pi)\} \\ 
&=\{z(t)\colon t\in[0,2\pi)\}, 
\end{align*}
where 
\begin{equation*}
 z(t):=R(t)(\cos t,\sin t), \quad
 R(t):=
 \begin{cases}
 r_+(t)>0&\text{ for }t\in[0,\pi],\\
 |r_-(t-\pi)|>0&\text{ for }t\in[\pi,2\pi]. 
 \end{cases}
\end{equation*}
So, the level curve $C$ is closed and simple, and its interior is 
\begin{equation*}
 I(c):=\{r\,(\cos t,\sin t)\colon0\le r<R(t)\}. 
\end{equation*}
The main idea for the elliptic polynomial case is to consider, for all real $c\ge c_*(n,d_*,D)$, the two opposite infinitesimal sectors of the interior $I(c)$ of the simple closed curve $C=P^{-1}(\{c\})$ between the rays $t$ and $t+dt$ and between the rays $t+\pi$ and $t+\pi+dt$, where $t$ is the polar angle in the interval $[0,\pi)$. The centroid of the union of these two sectors of $I(c)$ is at (signed) distance 
\begin{equation*}
 d(t)\sim \frac23\,\Big(r_+(t)\frac{|r_+(t)|^2}{|r_+(t)|^2+|r_-(t)|^2}
 +r_-(t)\frac{|r_-(t)|^2}{|r_+(t)|^2+|r_-(t)|^2}\Big) \tag{2}
\end{equation*}
from the origin. 
Formula (2) follows because (i) the centroid of an infinitesimal sector of radius $r>0$ between the rays $t$ and $t+dt$ is at distance $\frac23\,r$ from the origin, (ii) the area of this sector is $\frac12\,r^2\,dt$, and (iii) the centroid of the union of the two sectors is the weighted average of the centroids of the two sectors, with weights adding to $1$ and proportional to the areas of the sectors, and thus proportional to the squared radii of the sectors. 
Simplifying (2), we get
\begin{equation*}
 d(t)\sim-\frac{2b(t)}{na(t)}. 
\end{equation*}
Averaging now over all the pairs of opposite infinitesimal sectors, we see that the centroid converges to 
\begin{align*}
 &-\int_0^\pi dt\,\frac{2b(t)}{na(t)}(\cos t,\sin t)\frac12\,\Big(\frac c{a(t)}\Big)^{2/n}
 \Big/\int_0^\pi dt\,\frac12\,\Big(\frac c{a(t)}\Big)^{2/n} \\ 
 &=-\int_0^\pi dt\,\frac{2b(t)}{na(t)}(\cos t,\sin t)\Big(\frac1{a(t)}\Big)^{2/n}
 \Big/\int_0^\pi dt\,\Big(\frac1{a(t)}\Big)^{2/n}. \tag{3}
\end{align*}

I have checked this result numerically for $P(x,y)=x^4 + y^4 + 3 (x - y)^4 + y^3 + x y^2 + 10 x^2$, getting the centroid $\approx(-0.182846, -0.245149)$ for $c=10^4$ and $\approx(-0.189242,-0.25)$ for the limit (as $c\to\infty$) given by (3). From the above reasoning, one can see that the distance of the centroid from its limit is $O(1/c^{1/n})$; so, the agreement in this numerical example should be considered good, better than expected. 

One may also note that in general the level sets $P^{−1}([0,c])$ will not be convex, even if $P$ is a positive elliptic homogeneous polynomial. E.g., take $P(x,y)=(x−y)^2(x+y)^2+h(x^4+y^4)$ for a small enough $h>0$. Here is the picture of this level set for $c=1$ and $h=1/10$:

Clearly, the shape of this level set does not depend on $c>0$. 
This non-convexity idea can be generalized, with 
$$P(x,y)=P_{k,h}(x,y)
:=\prod_{j=0}^{2k-1}\Big(x\cos\frac{\pi j}k-y\sin\frac{\pi j}k\Big)^2+h(x^{4k}+y^{4k})$$
for natural $k$ and real $h>0$. Here is the picture of the curve $P_{k,h}^{-1}(\{1\})$ for $k=5$ and $h=(3/10)^{4k}$:

A: $(y-x^2)^2+x^2\phantom{aaaaaaaaaaaaaaaaaaaaa}$
