Milky Way dwarf spheroidal satellites are a prime target for Dark Matter (DM) indirect searches. Recently the importance of possible long-range interactions has been recognized, as they can boost the expected DM gamma ray signal by orders of magnitude through an effect commonly known as the Sommerfeld enhancement. However, for such analyses precise modelling of DM phase-space distribution becomes crucial and can introduce large uncertainties in the final result. We provide a pioneering attempt towards a comprehensive investigation of these systematics. First, the DM halo profiles are constrained using Bayesian inference on the available stellar kinematic datasets with a careful treatment of observational and theoretical uncertainties. We consider both cuspy and cored parametric DM density profiles, together with the case of a non-parametric halo modelling directly connected to observable quantities along the line-of-sight. After reconsidering the study case of ergodic systems, the basic ingredient of all previous analyses, we investigate for the first time scenarios where DM particles are allowed to have anisotropic velocity distributions. Referring to a generalized J-factor, sensitive to velocity-dependent effects, an enhancement (suppression) with respect to the 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. This work presents a state-of-the-art analysis of the aforementioned effects and falls within the interest of current and future experimental collaborations involved in DM indirect detection programs
