This isn't yet a complete question; rather, I'm looking for a qual-level question and answer describing a gravitationally bound system in an expanding universe. Since it's qual level, this needs a vastly simplified model, preferably one in a Newtonian framework (if this is even possible) that at least shows the spirit of what happens. If it isn't possible, and there is some important principle that any Newtonian model will miss, that is probably the first thing to point out. Otherwise, …
I'm thinking of a setup that starts something like this:
A body of mass $m$ orbits another of mass $M \gg m$ (taken to be the origin) in a circle of radius $R$ (and consequently with a period $T = 2 \pi \sqrt{ R^3 / GM }$). At time $t = 0$, space begins to slowly expand at a constant Hubble rate $H$ (where "slowly" means $ HT \ll 1 $, so that the expansion in one period is negligible compared to the orbital radius).
Now here is where I'm not sure how to phrase the problem. I'm thinking instead of something like the "bead on a pole" problems. Then the bead has generalized Lagrangian coordinates relative to a fixed position on the pole, but the pole is expanding according to some external function $a(t)$, a la FLRW metric (For a constant Hubble rate, $H = \dot a / a$ and so $ a(t) = e^{Ht} $). So, we can call the bead's distance from the origin in "co-moving units" the Lagrangian coordinate $q$. Then the distance of the bead from the origin is $d(t) = q(t) a(t)$, giving us a time dependent Lagrangian.
Basically, I need a 3D version of this that will work for the gravitationally bound system described. I also hope that the problem can be solved as an adiabatic process, which might mean changing coordinates to be around the original position relative to the more massive body at the origin, rather than co-moving coordinates.
EDIT 1: I won't claim that I fully understood the details of Schirmer's work, but I think one of the big takeaway points is that the cosmological expansion damages the stability of circular orbits and causes them to decay. I finished my very hand-waivey "Newtonian" model of a solar system in an expanding universe, and I don't think it captures the essence of this part of the GR solution. The model does have several sensical properties:
- There is a distance at which bound solutions are impossible and the two bodies are guaranteed to expand away from each other
- There is still a circular orbit (I haven't checked that it is stable) whose radius is modified by a term involving the Hubble constant.
- Using solar parameters, this shift is entirely negligible; it is larger (though still small) on a galactic scale.
Here is the model:
Without this scaling factor, the Lagrangian would be
$$ \mathcal{L} = \frac{1}{2} m \left( \dot r^2 + r^2 \dot \theta^2 + r^2 \sin^2 \theta \, \dot \phi^2 \right) + \frac{GMm}{r} $$
Adding in the scale factor, all the kinetic terms will pick up a factor $e^{2Ht}$ and the potential a factor $e^{-Ht}$:
$$ \mathcal{L} = \frac{1}{2} m \left(\dot r^2 + r^2 \dot \theta^2 + r^2 \sin^2 \theta \, \dot \phi^2 \right) e^{2Ht}+ \frac{GMm}{r} e^{-Ht} $$
(I have not added in any derivatives of the scaling factor; this would mean I was just changing coordinates and not adding in an externally controlled expansion of space). Since the problem is spherically symmetric, the total angular momentum is conserved and the motion lies in a plane (the expansion of space does not change this) at $\theta = \pi/2$. This reduces the effective Lagrangian to
$$ \mathcal{L} = \frac{1}{2} m \left( \dot r^2 + r^2 \dot \phi^2 \right) e^{2Ht} + \frac{GMm}{r} e^{-Ht} $$
Now let's change coordinates to $r = \eta e^{-Ht}$. This represents a point that is stationary with respect to the origin (i.e. the massive body that the smaller body orbits). Then $\dot r e^{Ht}= \dot \eta – H \eta $, so the Lagrangian becomes
$$ \mathcal{L} = \frac{1}{2} m \left((\dot \eta – H \eta)^2 + \eta^2 \dot \phi^2 \right) + \frac{GMm}{\eta} $$
This looks like a Lagrangian that has been modified slightly by the parameter $H$. After these manipulations, $\phi$ is still a cyclic coordinate, and its conjugate momentum $p_\phi = m\eta^2 \dot \phi = L$ is still conserved. The momentum conjugate to $\eta$ is
$$ p_{\eta} = m (\dot \eta – H \eta)$$
So, the equation of motion for $\eta$ is
$$ m \frac{d}{dt} (\dot \eta – H \eta) = m (\ddot \eta – H \dot \eta) = – mH(\dot \eta – H \eta) + m \eta \dot \phi^2 – \frac{GMm}{\eta^2} $$
Canceling terms and substituting in the equation of motion for $\phi$, this becomes
$$ \ddot \eta = H^2 \eta + \frac{L^2}{m^2 \eta^3} – \frac{GM}{\eta^2} $$
So, the expansion of space shows up as a force term proportional to $\eta$. For large $\eta_0$ at time $t = 0$, there is a solution $\eta(t) = \eta_0 e^{Ht}$, where $\dot \eta = H \eta$ agrees with Hubble's law (the corresponding comoving coordinate is $r(t) = \eta_0$). There is also an a highly unstable orbital solution for large $\eta$, where gravitational force just balances cosmological expansion. This new term also shifts the location of the "stable orbit" slightly (I haven't checked that it actually is stable in this situation). The location of the stable circular orbit when $H = 0$ is
$$ \eta_0 = \frac{L^2}{GMm^2} $$
Then let
$$ \eta = \eta_0 + \delta \eta = \eta_0 \left( 1 + \frac{\delta \eta}{\eta_0} \right)$$
To lowest order, the shift in location of the stable orbit is
$$ \frac{\delta \eta}{\eta_0} = \frac{H^2}{GM/\eta_0^3 – H^2} $$
For the orbit of the earth around the sun, this shift would be entirely negligible- about 15 pm. Taking the mass to be the mass of the Milky way and the distance to the distance of the Sun from the center, the shift would be a bit larger – a fractional amount of about 2e-7, which is several hundred AU – but still largely small.
Best Answer
I wrote an unpublished paper on this once, and presented the poster at an AAS meeting. You can actually do some relatively simple model building with this if you look at Schwarzschild-de Sitter solutions of Einstein's equation.${}^{1}$
Once you have the exact solution, you can then look at circular orbits, and do the usual stability analysis that you use to derive the $r > 6M$ limit that you get in the Schwarzschild solution. What you will find is that, in the Schwarzschild-de Sitter solution, you will get a cubic equation. One solution will give you a $ r < 0$ value, which is unphysical, but you will also get an innermost stable circular orbit and an outermost stable circular orbit. The latter represents matter getting pulled off of the central body by the cosmology.
There are also exact solutions that involve universes with nonzero matter densities other than the cosmological constant with a central gravitating body. These solutions generically don't have a timelike killing vector nor globally stable orbits. I ran numerical simulations of these orbits and published them in my dissertation's appendices.
EDIT: stack exchange seems to be killing hotlinks, but search for my name on the arxiv if you're interested.
${}^{1}$EDIT 2: you can find the Schwarzschild-de Sitter solution by replacing $1 - \frac{2M}{r}$ everywhere with $1-\frac{2M}{r} - \frac{1}{3}\Lambda r^{2}$. It should be easy enough to prove that that solution satisfies Einstein's equation for vacuum plus cosmological constant.