Solve a PDE defining a limit cycle of a nonlinear DE
Let me summarize here some ideas from comments and one weird finding of mine. This is not a complete answer, but I think I have found some explanation for what's happening.
As @Dmitry told, numerics show that system has a rotational limit cycle. Standard ways to prove existence of limit cycle for some parameter values include understanding what bifurcation could give rise to this limit cycle. To get a grasp on which bifurcation exactly happens it is reasonable to vary parameters and see what happens. Changing $k$ is a good idea: when $k=0$ the system is particularly simple and some analytics could be done. Naturally I was trying to understand what happens with limit cycle when you decrease $k$ to $0$. This is important because when $k = 0$ no such limit cycle is possible, I'll explain why further. The weird finding is that the more I decrease $k$ to $0$, the higher and higher this rotational limit cycle goes. That's weird and slightly suggests that bifurcation happens at infinity when $k = 0$. I am not good at dealing with bifurcations at infinity, for me it looks like some sort of weird "Andronov-Hopf at infinity", but I've never met myself such bifurcation before.
I think the picture is a bit more tractable when you start increasing $k$: at some value close to $k = 0.215$ it seems that rotational limit cycle collides with a heteroclinic trajectory connecting two saddles. It would be a good idea to figure out how stable separatrices of saddles behave before and after this bifurcation: they are always part of boundaries between different attraction basins and can help figure out multistability. You know about stable focus in this system (which is unique on the cylinder), thus the presence of another attractor might mean that a limit cycle exists.
Here I'll try to explain why rotational closed trajectories exist when $k \approx 0$. I'll start with "unperturbed" version of OP's equations, i.e. when $k = 0$. In that case system takes a form $\dot{x} = y, \, \dot{y} = -f(x)$. If you consider an equation $$\ddot{x} + f(x) = 0,$$ where $f(x)$ is at least continuous, it is well know that all such equations are integrable: the first integral is simply $\frac{\dot{x}^2}{2} + F(x)$, where $F(x)$ is such that $F'(x) \equiv f(x)$. If $f(x)$ is periodic, then we can consider an equivalent system $$ \dot{x} = y, $$ $$ \dot{y} = -f(x),$$ which is naturally a system on a cylinder. As a system on plane it also has a first integral which is $\frac{y^2}{2} + F(x)$. Note that $F(x)$ might be not periodic and hence the system on cylinder wouldn't have a first integral. Using the system on the plane and the first integral we can compute Poincaré map for $x = 0$ and $x = 2 \pi$ as cross-sections. Namely, $$ \frac{\lbrack y(0) \rbrack^2}{2} + F(0) = \frac{\lbrack \overline{y} \rbrack^2}{2} + F(2\pi),$$ or $$ \overline{y} = \sqrt{y^2 - 2(F(2\pi)-F(0))}$$ if we start with $y > 0$.
For the parameter values that @Dmitry mentioned in comments $F(2\pi) - F(0)$ is negative and doesn't depend on $k$. So I'll just write map as $\overline{y} = \sqrt{y^2 + \alpha}$, where $\alpha > 0$. Any rotational closed trajectory would correspond to a fixed point of this Poincaré map, i.e. to a solution of $y = \sqrt{y^2 + \alpha}$. It is quite obvious that there is no solution to this equation because $\sqrt{y^2 + \alpha} > \sqrt{y^2} = y$. However what happens, if we perturb this mapping a bit? For example, does $\beta + \sqrt{y^2 + \alpha} = y$ has solutions for $\beta \approx 0$? The answer is "no" when $\beta > 0$, but when $\beta < 0$ the answer is "yes". The function $\beta + \sqrt{y^2 + \alpha}$ has line $y + \beta$ as an asymptote, thus $(\beta + \sqrt{y^2 + \alpha}) - (\beta + y)$ tends to $0$ as $y \rightarrow +\infty$. From this follows that $\beta+\sqrt{y^2 + \alpha} - y \rightarrow \beta$, thus it is negative at some values of $y$, but positive when $y = 0$ it is positive. An existence of fixed point follows from continuity. Note that this fixed point ceases to exist when $\beta = 0$, but it exists for small $\beta < 0$: smaller $\beta$ corresponds to bigger coordinate of fixed point, "it comes from infinity".
My idea is that although the real perturbation of Poincaré map would be much different than my model example, something quite similar happens when $k$ is non-zero in your system. Probably it could be made more rigorous, but I didn't delve much into it.
This is an illustration of level sets for $\frac{y^2}{2} + F(x)$. Black level sets contain saddle equilibria. $F(x)$ was evaluated numerically using quad function from SciPy.
Note that the point at which level set intersects $x = 0$ is lower than the point of intersection with $x = 2\pi$. It happens due to aperiodicity of $F(x)$; in particular, $F(2\pi) - F(0) < 0$.
There is no limit cycle. For $k>0$ the system either moves continuously in one direction towards infinity or spirals down toward a stationary point.
Your system can be put back together as a second order scalar equation of the type $$ \ddot y+k\dot y +P'(y)=0 $$ for some potential function $P$. That is, loosely speaking, the system describes a mechanical motion of some object in the landscape of height $P(y)$ with a linear friction with friction coefficient $k$. The intuitive behavior is that this system will continuously lose energy to friction. If it gets trapped in a valley of the potential function it will oscillate inside that valley with continuously descending maxima until it settles at one of the local minima of $P$. Or the initial energy level, that is, the initial velocity, is so high that the object never settles and eventually moves towards infinity.
To derive that result formally, multiply the equation with $\dot y$ and integrate to find that the energy function $$ E(t)=\frac12\dot y(t)^2+P(y(t))=E(0)-k\int_0^t\dot y^2(s)\,ds $$ is constantly declining (as long as it remains in motion, which is indefinitely) towards one of the points with minimal value of $P(y)$ and $\dot y=0$.