We present an extended version of the Constant Time Lag analytical approach for the tidal evolution of circumbinary planets introduced in our previous work. The model is self-consistent, in the sense that all tidal interactions between pairs are computed, regardless of their size. We derive analytical expressions for the variational equations governing the spin and orbital evolution, which are expressed as high-order elliptical expansions in the semimajor axis ratio but retain closed form in terms of the binary and planetary eccentricities. These are found to reproduce the results of the numerical simulations with arbitrary eccentricities very well, as well as reducing to our previous results in the low-eccentric case. Our model is then applied to the well-characterised Kepler circumbinary systems by analysing the tidal timescales and unveiling the tidal flow around each different system. In all cases we find that the spins reach stationary values much faster than the characteristic timescale of the orbital evolution, indicating that all Kepler circumbinary planets are expected to be in a sub-synchronous state. On the other hand, all systems are located in a tidal flow leading to outward migration; thus the proximity of the planets to the orbital instability limit may have been even greater in the past. Additionally, Kepler systems may have suffered a significant tidally induced eccentricity damping, which may be related to their proximity to the capture eccentricity. To help understand the predictions of our model, we also offer a simple geometrical interpretation of our results.

We present a self-consistent model for the tidal evolution of circumbinary planets. Based on the weak-friction model, we derive expressions of the resulting forces and torques considering complete tidal interactions between all the bodies of the system. Although the tidal deformation suffered by each extended mass must take into account the combined gravitational effects of the other two bodies, the only tidal forces that have a net effect on the dynamic are those that are applied on the same body that exerts the deformation, as long as no mean-motion resonance exists between the masses. We apply the model to the Kepler-38 binary system. The evolution of the spin equations shows that the planet reaches a stationary solution much faster than the stars, and the equilibrium spin frequency is sub-synchronous. The binary components evolve on a longer timescale, reaching a super-synchronous solution very close to that derived for the 2-body problem. After reaching spin stationarity, the eccentricity is damped in all bodies and for all the parameters analyzed here. A similar effect is noted for the binary separation. The semimajor axis of the planet, on the other hand, may migrate inwards or outwards, depending on the masses and orbital parameters. In some cases the secular evolution of the system may also exhibit an alignment of the pericenters, requiring to include additional terms in the tidal model. Finally, we derived analytical expressions for the variational equations of the orbital evolution and spin rates based on low-order elliptical expansions in the semimajor axis ratio and the eccentricities. These are found to reduce to the 2-body case when one of the masses is taken equal to zero. This model allow us to find a close and simple analytical expression for the stationary spin rates of all the bodies, as well as predicting the direction and magnitude of the orbital migration.