Tidal torque drives the rotational and orbital evolution of planet–satellite and star–exoplanet systems. This paper presents one analytical tidal theory for a viscoelastic multi-layered body with an arbitrary number of homogeneous layers. Starting with the static equilibrium figure, modified to include tide and differential rotation, and using the Newtonian creep approach, we find the dynamical equilibrium figure of the deformed body, which allows us to calculate the tidal potential and the forces acting on the tide generating body, as well as the rotation and orbital elements variations. In the particular case of the two-layer model, we study the tidal synchronization when the gravitational coupling and the friction in the interface between the layers is added. For high relaxation factors (low viscosity), the stationary solution of each layer is synchronous with the orbital mean motion (n) when the orbit is circular, but the rotational frequencies increase if the orbital eccentricity increases. This behavior is characteristic in the classical Darwinian theories and in the homogeneous case of the creep tide theory. For low relaxation factors (high viscosity), as in planetary satellites, if friction remains low, each layer can be trapped in different spin-orbit resonances with frequencies $$n/2,n,3n/2,2n,\ldots $$ . When the friction increases, attractors with differential rotations are destroyed, surviving only commensurabilities in which core and shell have the same velocity of rotation. We apply the theory to Titan. The main results are: (i) the rotational constraint does not allow us to confirm or reject the existence of a subsurface ocean in Titan; and (ii) the crust-atmosphere exchange of angular momentum can be neglected. Using the rotation estimate based on Cassini’s observation (Meriggiola et al. in Icarus 275:183–192, 2016), we limit the possible value of the shell relaxation factor, when a deep subsurface ocean is assumed, to $$\gamma _s\lesssim 10^{-9}\,\hbox {s}^{-1}$$ , which corresponds to a shell’s viscosity $$\eta _s\gtrsim 10^{18}\,\hbox {Pa}\,\hbox {s}$$ , depending on the ocean’s thickness and viscosity values. In the case in which a subsurface ocean does not exist, the maximum shell relaxation factor is one order of magnitude smaller and the corresponding minimum shell’s viscosity is one order higher.