  Detailed knowledge of the different classes of stellar orbits that can be
accommodated in a given galactic potential is a prerequisite when building
self-consistent models using for instance the Schwazschild technique.
Furthermore, observational properties of galaxies depend on what these classes
of orbits are and on the presence of chaos in the systems. In the realistic
case in which the starting point for modeling is not a gravitational potential,
but an observed density distribution, we will require a gravitational theory to
make the connection between the stars that we see and the movement these stars
may be having. The argument can be turned upside down: understanding what
orbits may be allowed by each gravitational theory may give us a greater
insight on what these theories are and on how we can test them. Our aim is thus
to understand novel properties of orbits that are predicted by the latest
extension of the MOND phenomenology into the relativistic world. We integrated
orbits numerically in a fixed density distribution. The potential required for
such integration was obtained also numerically by assuming different
gravitational models. We find that thanks to the presence of a mass term in the
field equations, the theory can allocate new classes of orbits that do not
exist in Newtonian gravity nor standard MOND. We discuss consequences that
these new families of orbits can have in non-linear cosmological structure
formation as well as explore a possible alternative model for galactic
structure based on them.
