$$
\newcommand{\d}{\partial}
\newcommand{\r}{\mathbf{r}}
$$
Yes, there are many other terms involving higher powers of $\omega$, though
the powers of $\omega$ are more often encountered as $\cos 2 \omega$, $\sin
3 \omega$, etc., and the size of most of them is usually ignorable,
because they are multiplied by things which are small, like $J_2$ cubed or
$J_5$. There are recipes for generating however many terms you want to use,
with arbitrarily high powers of all the orbital elements, but they are
complicated to explain. For the first part of the derivation, you might
read Robert Battin's An
Introduction to the Mathematics and Methods of Astrodynamics (1999), where
this material is most of Chapter 10. It's also in pretty much every other
solid astrodynamics text you might like, though all with slightly different
notation. The second part is not covered well in Battin, Vallado, or other
such books I've read, with the exception of William Kaula's Theory of
Satellite Geodesy, so I recommend you try to find a copy of that.
I tried to boil down the derivation to something I could post here, but I
couldn't manage it. The algebra is just too dense. The main idea is to
write the full dynamics as the simple Keplerian two-body solution plus
something else, where all of the extra accelerations are known and the task
is to find the change to the two-body motion which matches those
perturbations. This is called variation of
parameters, which
is a general method for solving certain types of differential equations,
invented by Euler and extended by Lagrange and Gauss, specifically to treat
this problem in determining the orbits of Jupiter and Saturn. They figured
out the equations describing how the orbital elements change with time,
based on all of the things other than the simple Keplerian two-body
point-mass that you want to consider. That results in Lagrange's planetary
equations.
Of the six, the one you asked about is
$$
\frac{d \omega}{d t}=\frac{-\cos i}{\sqrt{\mu a (1-e^2)} \sin i} \frac{\d R}{\d i} + \frac{\sqrt{1-e^2}}{e \sqrt{\mu a}} \frac{\d R}{\d e}
$$
in Lagrange's form, where $R$ is the ''disturbing function'', whose
derivatives give the accelerations we want to study. In the case where the disturbance is caused by a third gravitating body, such as the sun or moon, $R$ equals
$$
\frac{Gm}{\rho}\left(\frac{\rho}{d} - \frac{r \cos \alpha}{\rho}\right)
$$
where $m$ is the third body's mass, $\rho$ is its distance from the central body, $d$ is its distance from the orbiting body, and $\alpha$ is the angle between the
vectors $\r$ and $\rho$. Finding an $R$ that gives rise to the acceleration you wish to study can be a
difficulty, so Gauss's alternate form uses the accelerations directly, but
I'm not going to introduce it here. Instead, let's look at where we get $R$
and how we take its derivatives. If we're looking at the effect of the
central body not being a point mass, then it is traditional to write $R$ as
$$
\frac{GM}{r} \sum_{n=1}^\infty\left(\frac{a}{r}\right)^n \sum_{m=0}^n \bar{P}_{nm}(\sin\phi)[\bar{C}_{nm}\cos(m\lambda)+\bar{S}_{nm}\sin(m\lambda)]
$$
as described in this
answer.
Keep in mind that we didn't have to do it this way. We chose to do it
this way, because this is actually a lot easier than many other ways we
could have chosen! The coefficients $\bar{C}_{nm}$ and $\bar{S}_{nm}$ are
defined as the result of integrating appropriate weighting functions times
the gravitational potential, with domain the mass distribution of the
central body, but they are actually measured by comparing this formula to
the observed motion of satellites. This process is the "satellite geodesy"
of Kaula's title.
What I called "the second part" above is the next step, which is actually
taking the derivatives of the disturbing function, and then fiddling around
until we come up with some approximation we can solve. Much of the work at
this point is classical mechanics at the level of Herbert Goldstein's
book,
making a bunch of "canonical"
transformations,
which really means changing variables until you transform the
Hamiltonian into
something with only trivial solutions. This is the origin of mean element
theories like Brouwer's
work that
began with Delaunay
variables
and went on from there. The other half,
which I've only seen worked out to this level of detail in Kaula, is to
rewrite various pieces of the disturbing function as things whose
derivatives can more easily be taken. It so happens that the sine of the
latitude $\phi$ equals $\sin i \sin(\omega + f)$, where $f$ is the true
anomaly, and $\sin \phi$ is the argument of the associated Legendre
polynomials! Furthermore, one can rewrite the terms in sin and $\cos m
\lambda$ (longitude) as polynomials in sin and cos of ($m(\Omega-\theta)+c(\omega+f)$),
where c is an integer depending on the power of cos or sin, and $\theta$ is
Greenwich hour angle. This means every term in the expansion contains
$\omega$ explicitly, except for the few that are symmetric about the
rotational axis of the spherical coordinate system.
================================================================
As an aside, something I find even more interesting happens at the so-called critical inclination: if $\cos^2 i = 1/5$, which happens at approximately 63.435 degrees, then by that formula $\dot{\omega}$ is exactly zero, no matter what value $J_2$ takes, and thus the argument of perigee does not change, no matter what body is being orbited.
This is used practically for the so-called Molniya orbit, a highly-eccentric (typically $e$ is 0.6 to 0.75) orbit with argument of perigee at 270 degrees (near the south pole), so that the apogee is near the north pole, and could thus give much better geographic coverage of the Soviet Union than could a (necessarily near-equatorial) geostationary orbit. High eccentricity exploits Kepler's equal area equal time law to arrange it so that the satellite spends most of its 12-hour orbit moving very slowly across the northern sky, and then zooms very quickly around the southern half of the globe before popping back up over the north where it wants to be and slowing down again. In such an orbit, you need $\dot{\omega}$ to vanish for the whole concept to work, or else your satellite that starts out dwelling over the northern hemisphere will start rotating around so that apogee becomes first equatorial and then over the southern hemisphere, where you don't need the communication coverage.
Now, of course $\dot{\omega}$ is never exactly zero, so these vehicles do have to do their own kind of stationkeeping, but the forces on them that they maneuver to cancel are very different than those found in other kinds of named orbits. The exact nature of this oddity has attracted a great deal of research interest over the years, some with a rather frightening degree of mathematical sophistication. To give a flavor of the sort of thing one may find, here are links to three very different papers that touch on interesting parts of the issue.
da Costa, de Moraes, Carvalho, and Prado, "Artificial satellites
orbiting planetary satellites: critical inclination and
sun-synchronous
orbits",
Journal of Physics: Conf. Series 911 (2017) 012018 is modern,
accessible, and discusses Earth's moon, Io, and Europa.
Jupp, "The critical inclination problem - 30 years of
progress",
Celestial Mechanics 43 (1989) 127-138, is an engagingly chatty history of the work in the 60s and 70s which provides a
somewhat gentle introduction to some of the advanced mathematics
necessary to understand my third example,
Coffey, Deprit, and Miller, "The Critical Inclination in Artificial
Satellite
Theory",
Celestial Mechanics 39 (1986) 365-406, which starts with the words
"the critical inclination in the main problem of artificial satellite
theory is an intrinsic singularity. Its significance stems from two
geometric events in the reduced phase space on the manifolds of
constant polar angular momentum and constant Delaunay action. In the
neighborhood of the critical inclination, along the family of circular
orbits, there appear two Hopf bifurcations, to each of which there
converge two families of orbits..."
and gets harder from there.
