What kind of effect would cosmic expansion have on planetary motion or keplerian orbit? The question is, if two following phenomena are thought to effect to the motion of the planet that orbits a star:


*

*Celestial body is pulled by some massive star by ordinary Newtonian gravity:


$$ m \frac{d^2r}{dt^2} = -  G \frac{Mm}{r^2}   $$


*There is additional radial velocity component


$$ \frac{dr}{dt} = kr\ ,\qquad k << 1 $$
 or (linear approximation):
$$ r(t) = r_0 (1+ kt) , k << 1\ , $$
that pulls the body $m$ away from the mass $M$ - and that is independent from the gravity (that means it happens also in zero gravity), what kind of solutions would these 2 rules give for orbits and can the orbits be stationary in some situations?

The Newtonian gravitation law only, when the mass $M \gg m$, would result a Keplerian orbit: $$ r(\theta) = \frac{a(1-e^2)}{1 + e\cos(\theta)}$$
that is circle, ellipse, hyperbola, or parabola.

Partial answers:
A. At least the circular motion equilibrium can be maintained at the start simply by demanding that the initial velocity of the body to be not exactly tangential, but have a small radial component inward:
$$ v_{initial} = v_{tan} + v_{r}\ , $$ where
$$ v_{r} = -kr  $$
Now the radial component of initial velocity would cancel the radial velocity $kr$ at initial time. What would yield to the familiar equilibrium equation for circular motion:
$$ v^2/r = GM/r^2 $$
B. Also, when the planet is put into aphelion point of the orbit, where the Newtonian gravitation is in then minimum and according Newtonian gravitation law its velocity would be perpendicular to $r$ and in the minimum, the initial velocity should be slightly inward from the tangential velocity $v$, if this point would be the true aphelion point. From this initial condition, the first guess is that the resultant orbit or trajectory should be then somewhere between the aphelion circle $r_{max}$ and perihelion circle $r_{min}$ -at least during the first cycle.

The first guess for the orbit equation is simply:
$$ r(\theta) = \frac{a(1-e^2)}{1 + e\cos(\theta)}(1+kt)$$
or
$$ r(\theta) = \frac{a(1-e^2)}{1 + e\cos(\theta)}e^{kt}\ ,$$
if nothing else (e.g gravity of other planets, orbital decay, etc) effect on the motion.

The reason why I ask this question is that I would like to think what would happen if the cosmological expansion of the universe would effect on the celestial mechanics. I have read that this is not thought to be true however, cosmic expansion is not effecting the solar system. But if it would be true, the cosmic expansion would be described approximately by equations:
$$r = r_0(1 + Ht) $$ that is linear approximation when $t$ is small, or $$ \frac{dr}{dt} = Hr\ , $$ which is another approximation when $t$ is small, whose solution is exponential function. 
And the value for $H$ would be around $H = 2.20 \times 10^{-18}\, \mathrm{s}^{-1}$ or $H = 6.93 \times 10^{-11}\, \mathrm{year}^{-1}$, that is simply the Hubble constant is the units [1/s] and [1/year].
With this value for $k$, for example Earth-Moon distance 384000 km should  increase by 2.63 cm/year, and Earth-Sun distance $150\times10^6\, \mathrm{km}$ would increase by 10.4 m/year, if nothing cancels this effect. 
The observed values are 3.8 cm/year (radar measurements) and 10.4 cm/year (That is 100 times smaller. I don't know how the value of AU is actually measured). The Moon-Earth system is thought to separate due to tidal forces, but it is unclear to me how much this effect really contribute to the increase of Earth-Moon distance.
 A: Your question has some bearing on what some people interpret erroneously as the source of the Pioneer anomaly. As some people point out there is some issue with what happens with a solar system in a galaxy. Really the influence of the cosmological constant is most likely to occur on that scale instead of a stellar system of planets. I will set this up some and address some aspects of this question.
I will use the stationary metric for the de Sitter spacetime that contains a central gravitating body
$$
ds^2~=~\left(1~-~\frac{2m}{r}~-~\frac{\Lambda r^2}{3}\right)dt^2~-~\frac{1}{\left(1~-~\frac{2m}{r}~-~\frac{\Lambda r^2}{3}\right)}dr^2~-~r^2(d\theta^2~+~sin^2\theta d\phi^2).
$$
Here $m~=~GM/c^2$ and the Hubble parameter is $H^2~=~\Lambda/3c^2$ and $\Lambda~\simeq~10^{-52}m^{-2}$. Now take the weak field approximation, and set the orbit in the plane with $\theta~=~\pi/2$
$$
ds^2~=~\left(1~-~\frac{2m}{r}~-~\frac{\Lambda r^2}{3}\right)dt^2~-~dr^2~-~r^2d\phi^2.
$$
The orbit in coordinate time with $rd\phi~=~v_\phi dt$ is contained in the general Lorentz gamma factor
$$
c^2\Gamma^{-2}~=~\left(\frac{ds}{dt}\right)^2~=~c^2\left(1~-~\frac{2GM}{rc^2}~-~\frac{\Lambda r^2}{3}\right)~-~v_r^2~-~v_\phi^2,
$$
where constants are all restored.
In special relativity we know that the kinetic energy of a particle is given by $K~=~(\gamma~-~1)mc^2$. We now consider the motion of a particle of mass $m$, not to be confused with $m~=~GM/c^2$ to advance the Hamiltonian in this weak field limit 
$$
H~=~\frac{p_r^2}{2m}~+~\frac{p_\phi^2}{2m}~-~\frac{GMm}{r}~-~\frac{mc^2\Lambda r^2}{6},
$$
which one can now use to solve dynamics
$$
\dot p_r~=~-\frac{\partial H}{\partial r}~=~\frac{1}{2m}\frac{\partial p_\phi^2}{\partial r}~-~\frac{GMm}{r^2}~+~\frac{mc^2\Lambda r}{3}.
$$
For circular motion $p_\phi~=~\omega r$ and the left hand side is zero. This recovers the standard textbook equation modified with the harmonic oscillator force in the opposite direction from the standard version. The cosmological constant then produces a spring like force that increase in the direction of extension.
If one were to pursue this then perturbation methods might be appropriate. Clearly for the radius $r$ appropriate for a stellar system of planets $r~\sim~10^{10}$ to $10^{11}m$ that $|\frac{GM}{r^2}|~>>~|\frac{c^2\Lambda r}{6}|$. The difference is around $22$ orders of magnitude, which means the role of the cosmological constant on stellar systems of planets is negligable. We may then consider this for an entire galaxy. Here these two terms are fairly comparable. However, we know that galaxies have dynamics that is not Keplerian. This is due to the fact stars in a galaxy are immersed in some dark matter. Let $\rho~=~M_d/vol$ be the density of this dark matter. Gauss's law gives us as a boiler plate approximation a force due to this dark matter
$$
4\pi r^2 F~=~4\pi GM_d~=~16\pi G\rho r^3/3
$$
so this force has magnitude $F~=~4\pi G\rho r/3$. Now one has to include this term which gives
$$
H~=~\frac{p_r^2}{2m}~+~\frac{p_\phi^2}{2m}~-~\frac{GMm}{r}~+~\frac{(8\pi G\rho~-~mc^2\Lambda)r^2}{6}.
$$
This term largely dominates the dynamics. In fact the dark matter force is about an order of magnitude larger. In this case of course the orbital dynamics of stars around galaxies, except for very close to the central black hole, is not at all Keplerian.
A: To continue this as an addendum, I played around with some of this. It is clear that for a circular orbit that the cosmological constant will only slightly adjust the radius. There will be no change in the radius of the orbit with time. For an elliptical orbit this might be different over a very long period of time. I will continue with this dynamical equation modified by $mc^2\Lambda r^2/6$ potential.  I rewrite the Hamiltonian as
$$
H~=~\frac{p_r^2}{2m}~+~\frac{L^2}{2mr^2}~-~\frac{GMm}{r}~+~\frac{mc^2\Lambda r^2}{6}.
$$
We are interested in the radial dynamics and the Hamilton equation is
$$
\frac{dp_r}{dt}~=~-\frac{\partial H}{\partial r}~=~\frac{L^2}{mr^3}\frac{\partial p_\phi^2}{\partial r}~-~\frac{GMm}{r^2}~+~\frac{mc^2\Lambda r}{3}.
$$
We change variables with $u~=~1/r$ and substitute $dt~=~d\phi/\omega$ $\omega~=~L/mr^2$ so that this differential equation is
$$
\frac{d^2u}{d\phi^2}~+~u~+~\frac{\lambda}{3\ell^2 u^3}~=~\frac{m}{\ell^2},
$$
and $\ell$ angular momentum per mass.  This equation is a perturbed form of the standard equation for Keplerian dynamics with $u(t)~=~c_1~\sin\phi~+~c_2~\cos\phi~+~C$.
We treat the perturbation $\frac{\lambda}{3\ell^2 u^3}$ by using 
$$
u~=~\frac{m}{\ell^2}(1~+~e~\cos\phi),~m~=~GM/c^2.
$$
For small eccentricity we have $u^{-3}~\simeq~\frac{m}{\ell^2}(1~+~3e~\cos\phi))$
and special solutions are used $\phi\sin\phi$ to find the solution 
$$
u(t)~\simeq~\frac{m}{\ell^2}\left(1~+~e~\cos\phi~+~{3\lambda}{\ell^3}e\phi~\sin\phi\right)~=~\frac{m}{\ell^2}\left(1~+~e~\cos\left(1~-~{3\lambda}{\ell^3}\phi\right)\right)
$$
This will create a precession on the order of $10^{-10}$ times the precession of Mercury's orbit. The next order term to $e^2$ will be $\sim~\frac{3\lambda}{\ell^3}e^2\phi^2~\sin^2\phi$ This will produce a change in the orbital radius with the initial radius $R$
$$
\frac{\overline\delta r(t)}{R}~\sim~\frac{3\lambda}{\ell^3}e^2\phi^2~\simeq~\frac{3\lambda}{\ell^3}e^2t^2.
$$
For an orbit comparable to Earth $\frac{\overline\delta r(t)}{R}~\sim~10^{-33}$ to $10^{-34}s^{-2}\times t^2$. A period of one billion years will change the radius by a factor of $10^{-2}$ to $10^{-1}$, or for the Earth about $0.01\,\mathrm{AU}$ to $0.1\,\mathrm{AU}$. This implies that Earth has migrated outwards since it origin. Below is a graph of this process that is greatly exaggerated.

This is an interesting subject to pursue. The early Sun when life emerged on Earth was only $0.75$ as luminous as today. This means that Earth would have received about as much solar radiation as Mars does today. The average temperature on Earth would have been about $-100^\circ \mathrm{C}$ and even with a thicker $CO_2$ atmosphere this still would have left Earth as a very frozen place. However, if Earth were at $0.85\,\mathrm{AU}$ early on the Earth would have received close to $80\%$ of current solar radiation, which with $CO_2$ heat trapping would have kept the early Earth from being an ice ball.
