  In previous work, we have found new static, spherically symmetric boson star
solutions which generalize the standard boson stars by allowing a particular
superposition of scalar fields in which each of the fields is characterized by
a fixed value of its non-vanishing angular momentum number $\ell$. We call such
solutions "$\ell$-boson stars". Here, we perform a series of fully non-linear
dynamical simulations of perturbed $\ell$-boson stars in order to study their
stability, and the final fate of unstable configurations. We show that for each
value of $\ell$, the configuration of maximum mass separates the parameter
space into stable and unstable regions. Stable configurations, when perturbed,
oscillate around the unperturbed solution and very slowly return to a
stationary configuration. Unstable configurations, in contrast, can have three
different final states: collapse to a black hole, migration to the stable
branch, or explosion (dissipation) to infinity. Just as it happens with
$\ell=0$ boson stars, migration to the stable branch or dissipation to infinity
depends on the sign of the total binding energy of the star: bound unstable
stars collapse to black holes or migrate to the stable branch, whereas unbound
unstable stars either collapse to a black hole or explode to infinity. Thus,
the parameter $\ell$ allows us to construct a new set of stable configurations.
All our simulations are performed in spherical symmetry, leaving a more
detailed stability analysis including non-spherical perturbations for future
work.
