The question is "how are singularities avoided for $e→0$ by defining $(e_x,e_y)$, which is still a function of $ω$, thus still affected by the singularity?".
(Worth noting first is that here $ω$ means argument of perigee, i.e. distance from node to perigee, and $e_x, e_y$ are rectangular components of the eccentricity vector. The formulation given in the question
$$(e_x, e_y) = (e.\cos\omega, e.\sin\omega)$$
also presupposes that the xy plane is the orbit plane, also that the x-axis itself is in the direction of the node (vector) -- but these are not conditions often used in practice.)
While $e_x, e_y$ can be defined through a relation similar to that given in the question, they are not in practice calculated that way.
It would perhaps be too limited to leave the matter here, without saying how in practice such quantities have more suitably been calculated. For that, there are numerous ways, involving various different sets of orbital elements.
Relations such as given in the question can arise e.g. in problems of perturbed elliptic motion, where it is necessary to calculate relationships between many instantaneous positions and velocities of an orbiting body, and corresponding instantaneous values of some chosen set of elliptical orbital elements. It was found long ago that older sets of Keplerian elements including mean anomaly, argument of perigee and longitude of node could give rise to numerical difficulties where $e\rightarrow 0$ or $i\rightarrow 0$.
The problems, and their avoidance, can be made clearer when discussed not so much in terms of singularities, but rather of conditions under which some relations between orbital elements and position and velocity of the orbiting body become poorly determined. The 'offending' elements where the inclination and/or eccentricity are small are the locations of the perigee and node.
Thus, the relation given in the question would not be used in practice to calculate $e_x, e_y$, rectangular components of the eccentricity vector, partly because the location $ω$ of the perigee becomes poorly determined in terms of position and velocity when the magnitude $e$ of the eccentricity is close to zero. Also, the location of the node $Ω$ becomes poorly determined in terms of the orbiting body's position and velocity when the inclination is close to zero. Both those occurrences leave the relation given in the question with poor numerical definition.
For such reasons, other sets of orbital elements not involving the node, the argument of perigee/perihelion, and mean anomaly have often been used in place of the older orbital elements, so that the singularities or poor determinacies of the locations of node and perigee no longer arise.
The topic receives good explanation in A E Roy 'Orbital motion', (2005, 4th ed., ch.8, p.239 onwards), in terms that the "usual transformation to avoid this disadvantage" [i.e. the singularities] "consists of adopting the variables $h, k, p, q$ instead of the offending elements $e, ϖ, i$ and $Ω$" (where variable $ϖ$ is $(ω + Ω)$), and
$$h = e.\sin ϖ$$ $$k = e.\cos ϖ$$ $$p = \sin i.\sin Ω$$ $$q = \sin i.\cos Ω .$$
Also, in place of mean anomaly $M$, mean longitude $\lambda = (M+ω+Ω)$ is used. So too are vectors $H$ and $e$, defined in terms of positions and velocities of the orbiting body in rectangular coordinates, where the instantaneous angular momentum or areal velocity, and eccentricity vector, are respectively
$H$ = ($r$ x $v$) ,
$e$ = ($v$ x $H$)/$\mu$ -$r$/r ,
with components and magnitudes respectively $H_x, H_y, H_z$ and $H_m$, and $e_x, e_y, e_z$ and $e_m$ .
As an example of principle to compare with the relation given in the question, the quantities $e_m.\cosω$ and $e_m.sinω$ would not be used, but instead the following related elements of the alternative set might be calculated when required from components of $e$ and $H$, in turn evaluated from the body's position and velocity:
$$k (= e_m.\cosϖ) = e_x - e_z.H_x/(H_m + H_z)$$
$$h (= e_m.\sinϖ) = e_y - e_z.H_y/(H_m + H_z)$$
(cf. Roy (2005) 8.30) (noting that if inclination i = 0 then e_z, H_x, and H_y are all zero as well), and these relations involve $ϖ$ i.e. $(ω + Ω)$ and do not depend on the individual parameters $ω, Ω,$ or $M$ which may become poorly-determined.
[Postscript/edit:]
A visual form of explanation how the elements $ ϖ = (\omega + \Omega)$ and $\lambda = (M+\omega+\Omega)$ can remain well-determined even when the location of the node $\Omega$ and/or argument $\omega$ of perihelion become nearly indeterminate is presented by the diagram below.
Angles $ϖ$ and $\lambda$ are often portrayed as if they were (artificial) sums of angles in two different planes. In reality, they are also sums of angles in the same plane, i.e. the orbit plane, taken to define auxiliary coordinate axes, here $(F G W)$, inclined to the $(x y z)$ coordinate axes in main use (e.g. ecliptic for planetary orbits, equatorial for earth-satellites). The $F G$ plane is the orbit plane, and its intersection/node $N$ with the $x y$ plane is equally distant, by angle $\Omega = xN = FN$, from the initial points for reckoning longitudes in both planes, i.e. the $x$-vector of the $(x y z)$ coordinates and the $F$-vector of the $(F G W)$ coordinates. Argument of perigee $\omega$ is angle $NP$, and the mean position of the orbiting body is at $B$ with mean anomaly angle $PB$ = $M$.
It can be seen that if the inclination angle $i$ at node $N$ becomes very small, and the location of $N$ poorly determined, any error in the location of $N$ leaves the angle $FP$ practically unaffected, $XN$ and $FN$ remaining equal by construction. So the errors in angles $\omega$ and $\Omega$ are practically equal and opposite, and their sum, the distance $FP$ = $ϖ = (\omega + \Omega)$ practically unaffected. On similar considerations, $\lambda = (M + \omega + \Omega)$, angle $FB$, can remain well-determined even with poorly-determined locations of intermediate points $N$ and $P$.
The $(F G W)$ coordinate system is described for example in R M L Baker (1967) "Astrodynamics: Applications and Advanced Topics", especially Appendix C (without diagram), where $x- y- z-$ coordinates of unit vectors $W$ and $F$ are derived as
$ W_x = + \sin \Omega . \sin i $ ( = $H_x/Hm$),
$ W_y = - \cos \Omega . \sin i $ ( = $H_y/Hm$),
$ W_z = \cos i $ ( = $H_z/Hm$);
$ F_x = 1 - W_x^2 / (1 + W_z) $ ,
$ F_y = - W_x W_y / (1 + W_z) $ ,
$ F_z = - W_x $ ;
and $G$ = ($W$ x $F$).
The eccentricity vector $e$ can also be expressed by elements $k$ and $h$, seen to be its components resolved in the directions of the orthogonal $F-$ and $G-$ axes:
$(k, h) = (e_F, e_G) = (e_m.cosϖ, e_m.sinϖ)$ .
This form of relation is clearly analogous to the formulation given in the question.
It has in addition the desirable property that the quantities involved do not become poorly-determined at low inclinations and/or eccentricities, as they can be computed as shown above from the elements of vectors $H$ and $e$ derived from the orbital positions and velocities.
[end edit:]
A number of related alternative sets of orbital elements have also been proposed and used for similar reasons -- see the same chapter in Roy (2005) already cited; and a previous related set proposed in Broucke & Cefola (1972), 'Equinoctial orbital elements', Cel.mech. 5:303-310, which adopts
$$p = \tan i/2.\sin Ω$$ $$q = \tan i/2.\cos Ω$$
and a yet further variant set which has been used and which adopts
$$p = \sin i/2.\sin Ω$$ $$q = \sin i/2.\cos Ω .$$
[Postscript/edit:]
The last-mentioned variant of elements $p$ and $q$, based on the factor $\sin i/2$, has been effectively used in celestial-mechanics studies at the Observatoire de Paris/Bureau des Longitudes (see "Introduction aux éphémérides astronomiques: supplément explicatif à la Connaissance des Temps", 1997, esp. ch.8). This set of elements has an advantage in generality for inclinations that may exceed 90° (in which case the ($p$,$q$) elements based on $\sin i$ can result in ambiguity for the value of $i$), or may reach 180° (where the elements based on the factor $\tan i/2$ encounter a singularity.
[end edit:]
Brouwer & Clemence (1961) 'Methods of Celestial mechanics', ch.XI on 'Variation of Arbitrary Constants' (a/k/a variation of parameters) (sec.7 at p.287 on) has further discussion of situations where alternative forms of the orbital elements are desirable in cases of small eccentricity and/or small inclination, and how they are used.
