  The two-dimensional $q$-state clock model for $q \geq 5$ undergoes two
Berezinskii-Kosterlitz-Thouless (BKT) phase transitions as temperature
decreases. Here we report an extensive worm-type simulation of the
square-lattice clock model for $q=$5--9 in a pair of flow representations, from
the high- and low-temperature expansions, respectively. By finite-size scaling
analysis of susceptibility-like quantities, we determine the critical points
with a precision improving over the existing results. Due to the dual flow
representations, each point in the critical region is observed to
simultaneously exhibit a pair of anomalous dimensions, which are $\eta_1=1/4$
and $\eta_2 = 4/q^2$ at the two BKT transitions. Further, the approximate
self-dual points $\beta_{\rm sd}(L)$, defined by the stringent condition that
the susceptibility like quantities in both flow representations are identical,
are found to be nearly independent of system size $L$ and behave as $\beta_{\rm
sd} \simeq q/2\pi$ asymptotically at the large-$q$ limit. The exponent $\eta$
at $\beta_{\rm sd}$ is consistent with $1/q$ within statistical error as long
as $q \geq 5$. Based on this, we further conjecture that $\eta(\beta_{\rm sd})
= 1/q$ holds exactly and is universal for systems in the $q$-state clock
universality class. Our work provides a vivid demonstration of rich phenomena
associated with the duality and self-duality of the clock model in two
dimensions.
