  A set of unified relativistic mean-field equations of state for hyperonic
compact stars recently built in [M. Fortin, Ad. R. Raduta, S. Avancini, and C.
Providencia, Phys. Rev. D {\bf 101}, 034017 (2020)] is used to study the
thermal evolution of non-magnetized and non-rotating spherically-symmetric
isolated and accreting neutron stars under different hypothesis concerning
proton $S$-wave superfluidity.
  These equations of state have been obtained in the following way: the slope
of the symmetry energy is in agreement with experimental data; the coupling
constants of $\Lambda$ and $\Xi$-hyperons are determined from experimental
hypernuclear data; uncertainties in the nucleon-$\Sigma$ interaction potential
are accounted for; current constraints on the lower bound of the maximum
neutron star mass are satisfied.
  Within the considered set of equations of state, the presence of hyperons is
essential for the description of the cooling/heating curves.
  One of the conclusions we reach is that the criterion of best agreement with
observational data leads to different equations of states and proton $S$-wave
superfluidity gaps when applied separately for isolated neutron stars and
accreting neutron stars in quiescence.
  This means that at least in one situation the traditional simulation
framework that we employ is not complete and/or the equations of state are
inappropriate.
  Another result is that, considering equations of state which do not allow for
nucleonic dUrca or allow for it only in very massive NS, the low luminosity of
SAX J1808 requires a repulsive $\Sigma$-hyperon potential in symmetric nuclear
matter in the range $U_\Sigma^{(N)}\approx 10-30$ MeV. This range of values for
$U_\Sigma^{(N)} $ is also supported by the criterion of best agreement with all
available data from INS and XRT.
