  Fully non-linear numerical evolutions of the Einstein-(multi)--Klein-Gordon
equations are performed to study head-on collisions of $\ell$-boson stars.
Despite being spherically symmetric, $\ell$-boson stars have a (hidden) frame
of reference, used in defining their individual multipolar fields. To assess
the impact of their relative orientation, we perform simulations with different
angles between the axes of the two colliding stars. Additionally, two scenarios
are considered for the colliding stars: that they are composites of either the
same or different scalar fields. Despite some model-specific behaviours, the
simulations generically indicate that: 1) the collision of two sufficiently
(and equally) massive stars leads to black hole formation; 2) below a certain
mass threshold the end result of the evolution is a bound state of the
composite scalar fields, that neither disperses nor collapses into a black hole
within the simulation time; 3) this end product (generically) deviates from
spherical symmetry and the equipartition of the number of bosonic particles
between the different scalar fields composing the initial boson stars is lost,
albeit not dramatically. This last observation indicates, albeit without being
conclusive, that the end result of these collisions belongs to the previously
reported larger family of equilibrium multi-field boson stars, generically
non-spherical, and of which $\ell$-boson stars are a symmetry enhanced point.
We also extract and discuss the waveforms from the collisions studied.
