  Milky Way dwarf spheroidal satellites are a prime target for Dark Matter (DM)
indirect searches. There have been recent reassessments of the expected DM
gamma-ray signals in case of long-range interactions, commonly known as
Sommerfeld enhancement. Since details of the underlying DM phase-space
distribution function become critical, there are potentially large
uncertainties in the final result. We provide here a first attempt towards a
comprehensive investigation of these systematics, addressing the impact on the
expected DM flux from Milky Way dwarfs via Bayesian inference on the available
stellar kinematic datasets. After reconsidering the study case of ergodic
systems, we investigate for the first time scenarios where DM particle orbits
may have a radial or tangential bias. We consider both cuspy and cored
parametric DM density profiles, together with the case of a non-parametric halo
modelling directly connected to line-of-sight observable quantities. The main
findings of our work highlight the relevance of the assumed phase-space
distribution: Referring to a generalized J-factor, namely the line-of-sight
convolution of the spatial part in case of velocity-dependent annihilation
rate, an enhancement (suppression) with respect to the limit of isotropic
phase-space distributions is obtained for the case of tangentially (radially)
biased DM particle orbits. We provide new estimates for J-factors for the eight
brightest Milky Way dwarfs also in the limit of velocity-independent DM
annihilation, in good agreement with previous results in literature, and derive
data-driven lower-bounds based on the non-parametric modelling of the halo
density. The outcome of our broad study stands out as a representative of the
state-of-the-art in the field, and falls within the interest of current and
future experimental collaborations involved in DM indirect detection programs.
