  We investigate the dynamical phases and phase transitions arising in a
classical two-dimensional anisotropic $XY$ model under the influence of a
periodically driven temporal external magnetic field in the form of a symmetric
square wave. We use a combination of finite temperature classical Monte Carlo
simulation, implemented within a CPU + GPU paradigm, utilizing local dynamics
provided by the Glauber algorithm and a phenomenological equation-of-motion
approach based on relaxational dynamics governed by the time-dependent free
energy within a mean-field approximation to study the model. We investigate
several parameter regimes of the variables (magnetic field, anisotropy, and the
external drive frequency) that influence the anisotropic $XY$ system. We
identify four possible dynamical phases -- Ising-SBO, Ising-SRO, $XY$-SBO and
$XY$-SRO. Both techniques indicate that only three of them (Ising-SRO,
Ising-SBO, and $XY$-SRO) are stable dynamical phases in the thermodynamic
sense. Within the Monte Carlo framework, a finite size scaling analysis shows
that $XY$-SBO does not survive in the thermodynamic limit giving way to either
an Ising-SBO or a $XY$-SRO regime. The finite size scaling analysis further
shows that the transitions between the three remaining dynamical phases either
belong to the two-dimensional Ising universality class or are first-order in
nature. The mean-field calculations yield three stable dynamical phases, i.e.,
Ising-SRO, Ising-SBO and $XY$-SRO, where the final steady state is independent
of the initial condition chosen to evolve the equations of motion, as well as a
region of bistability where the system either flows to Ising-SBO or $XY$-SRO
(Ising-SRO) depending on the initial condition. Unlike the stable dynamical
phases, the $XY$-SBO represents a transient feature that is eventually lost to
either Ising-SBO or $XY$-SRO.
