In their famous paper in 1985 (link), Damour&Deruelle describe the orbital motion for a binary system taking into account first-order post-Newtonian corrections (1PN). The solution is given in their eqs. (7.1-7.2), according to which the binary separation $a_R$ (just to make an example) is:
\begin{equation}\tag{1} a_R= - \dfrac{GM}{2E}\left[ 1- \dfrac{1}{2}(\nu-7)\dfrac{E}{c^2} \right] \end{equation}
where $\nu \doteq (m_1 m_2) / (m_1+m_2)^2$. However, in their Appendix B, the authors show the equation of motion by using the Lagrange's method of variation of constants. In this case, all the orbital parameters are functions of time: for example, the binary separation is now (see Eqs. B11):
\begin{equation}\tag{2} a(t)= \bar{a}\left[ 1+ \dfrac{GM}{c^2 p (1-e^2)}*f(e,\phi) \right] \end{equation}
where $\bar{a}$ is a constant of integration, $e$ is the eccentricity (Keplerian or PN?), $p=a(1-e^2)$, $\phi$ is the true anomaly, and $f$ is a given function of $e$ and $\phi$ that I omitted for brevity.
My question is:
(a) Which is the physical difference between $a_R$ and $a(t)$?
(b) I know that relativistic corrections make the orbital parameters time-dependent; so why (1) does not depend on time while (2) does?
(c) If a want to simulate the orbital motion at 1PN for 10 years, which set of eqs. should I use?
