A numerical study of the nonlinear torsional solutions of convection in rotating, internally heated, self-gravitating fluid spheres is presented. Their dependence on the Rayleigh number has been found for two pairs of Ekman (E) and small Prandtl (Pr) numbers in the region of parameters where, according to Zhang et al. [J. Fluid Mech. 813, R2 (2017)], the linear stability of the conduction state predicts that they can be preferred at the onset of convection. The bifurcation to periodic torsional solutions is supercritical for sufficiently small Pr. They are not rotating waves, unlike the nonaxisymmetric case. Therefore they have been computed by using continuation methods for periodic orbits. Their stability with respect to axisymmetric perturbations and physical characteristics have been analyzed. It was found that the time-and space-averaged equatorially antisymmetric part of the kinetic energy of the stable orbits splits into equal poloidal and toroidal parts, while the symmetric part is much smaller. Direct numerical simulations for E = 10 −4 at higher Rayleigh numbers (Ra) show that this trend is also valid for the nonperiodic flows and that the mean values of the energies remain almost constant with Ra. However, the modulated oscillations bifurcated from the quasiperiodic torsional solutions reach a high amplitude, compared with that of the periodic, increasing slowly and decaying very fast. This repeated behavior is interpreted as trajectories near heteroclinic chains connecting unstable periodic solutions. The torsional flows give rise to a meridional propagation of the kinetic energy near the outer surface and an axial oscillation of the hot nucleus of the metallic fluid sphere.