  Modeling the surface brightness distribution of stars is of prime importance
to interpret observations. Nevertheless, this remains quite challenging for
cool stars as it requires one to model the MHD turbulence that develops in
their convective envelope. In Paper I, the effect of the Coriolis acceleration
on the surface heat flux has been studied by means of hydrodynamic simulations.
In this paper, we aim to investigate the additional effect of dynamo magnetic
fields. We focus on an envelope thickness that is representative of either a
$\sim0.35~M_\odot$ M dwarf, a young red giant star or a pre-main sequence star.
We performed a parametric study using numerical MHD simulations of anelastic
convection in thick rotating spherical shells. For each model, we computed the
mean surface distribution of the heat flux, and examined the leading-order
effect of the magnetic field on the obtained latitudinal luminosity profile. We
identify three different regimes. Close to the onset of convection, while the
first unstable modes tend to convey heat more efficiently near the equator,
magnetic fields are shown to generally enhance the mean heat flux close to the
polar regions (and the tangent cylinder). By progressively increasing the
Rayleigh number, the development of a prograde equatorial jet was previously
shown to make the equator darker when no magnetic field is taken into account.
For moderate Rayleigh numbers, magnetic fields can instead inverse the mean
pole-equator brightness contrast (which means going from a darker to a brighter
equator when a dynamo sets in) and finally induce a similar regime to that
found close to the onset of convection. For more turbulent models with larger
Rayleigh numbers, magnetic fields alternatively tend to smooth out the
brightness contrast. This general behavior is shown to be related to the
quenching of the surface differential rotation by magnetic fields.
