Neutron Star Cooling in Covariant f(Q) Gravity

We present a study of the thermal evolution of neutron stars in covariant f(Q) gravity, extending our previous analyses of neutron star structure within this framework (Alwan, 2025). Stellar structure profiles are obtained from a modified Tolman–Oppenheimer–Volkoff solver incorporating three representative f(Q) models (logarithmic, exponential, and quadratic) and subsequently used as input to the NSCool thermal evolution code. Since f(Q) gravity modifies the underlying spacetime geometry, its effects propagate into the neutron star cooling behavior through changes in stellar structure rather than through direct modifications to the microphysical processes governing heat transport and neutrino emission. We find that these structural modifications naturally give rise to distinct slow, threshold, and rapid cooling regimes, whose boundaries shift relative to the general relativistic case depending on the specific f(Q) model and its parameters. The logarithmic and exponential models allow for shifts in either direction depending on the choice of model parameters, whereas the quadratic model produces shifts in only one direction, indicating that the f(Q) gravity sector produces characteristic features in the thermal evolution of neutron stars that are absent in general relativity. We also compare the resulting cooling curves with available observational data from isolated neutron stars (Beznogov & Yakovlev, 2015).