observations and vary ${M_\bullet}$ and $\Upsilon$ and use the relation $\rho_\mathrm{DM}(\Upsilon)$ determined from the first set of models [@Rusli13a]. In this way, we reduce the originally 4D parameter space to two smaller ones (3D and 2D, respectively).

In each trial potential, i.e., for each particular choice of ${M_\bullet}$, $\Upsilon$, and $\rho_\mathrm{DM}$, we compute an orbit library containing about $25\,000$ representative stellar orbits. Half of these orbits have positive angular momentum $L_z$ along the rotation axis; the rest have their $L_z$ flipped in the opposite direction. The contribution of each orbit to the best-fit model in the actual potential, the so-called orbital occupation number or weight, is determined by the maximization of $$\label{eq:maxs}
S-\alpha \chi^2,$$ subject to the constraints imposed by the deprojected galaxy luminosity density $\nu$. In Equation (\[eq:maxs\]), $S$ denotes the Boltzmann entropy and $\chi^2$ quantifies the deviations between the observed line-of-sight velocity distributions and the model. The entropy maximization and $\chi^2$ minimization is done iteratively for $\sim 30$ different values of the regularization parameter $\alpha$. We start with $\alpha=0$ and stop when an increase in $\alpha$ does not reduce the $\chi^2$ anymore. Models with different $\alpha$ reproduce the density constraints equally well but differ in their internal orbital 