Yes, good, you've saved me some typing along these lines. (Well, in retrospect, maybe not that much :D )
In typical vacuum solutions of the Einstein Field Equations the metric is what determines the available geodesics[0]. Putting an isolated object into a vacuum solution causes it to "select" an appropriate geodesic from the ones available, and if left alone, that object will at every time be somewhere on that geodesic. A "test object" is pointlike and massless, so adding it to the vacuum does not change the background metric at all. But if we make it heavier, or bigger and rotating or oscillating, then the spacetime is no longer exactly modelled by the vacuum solution. At that point one might use perturbation theory, and consider the vacuum metric as a background and the metric the not-quite-a-test-object itself generates as a perturbation field overlaid on that.
The algebra usually looks like g_{\mu\nu} = \eta_{\mu\nu} + h_{\mu\nu} where the greek subscripts are the usual indices running 0,1,2,3 in 4-dimensional Lorentzian spacetimes (like ours, and most that you'll ever run into), \eta is the chosen background field, h is the perturbation field, and g is the "true" metric. One can also add perturbation fields or play around with higher-order contributions from self-interactions (say if the not-a-test-object is a black hole binary, or a solid fragment of a supernova (an Earth-massed blob of hot metal?) moving relativistically).
If we start with an expanding Robertson-Walker spacetime vacuum (no matter, not even dark matter), and add a single galaxy cluster, we need to figure out the metric to see how it evolves (or to see what geodesics are available for test particles we throw in as probes).
We can do this in a couple of ways: either using perturbative methods like above, e.g. g = {RW} + {galaxy}, or start with a single g for the whole spacetime from a known set of appropriate solutions like Schwarzschild-de Sitter or Kottler, or by engaging in "inside"/"outside" stitching together of metrics like I described earlier (swiss-cheese). (There are other approaches too!)
Each has advantages and drawbacks.
The middle option has the problem that a galaxy cluster will only look like a Schwarzschild source from so far away that it shrinks to a point. Closer observers will see (optically, even!) that it's lumpy, and if they have sufficiently sensitive gravimeters or a probe like Synge's "five-point curvature detector" they will see that Schwarzschild is not quite accurate. In particular, the geodesics the components of the gravimeter or probe "find" are not the ones that would be generated in the model metric.
The first option is hard to do and usually leads one into numerical relativity. That's not such a bad thing these days, thanks to modern codes and supercomputers. However, it's often hard to extract intuitions from numerical solutions, and even harder to boil down into an explanation for others ("just run this code and you'll see" is not very satisfying).
The last option has been around for decades and is a bit half-way. Its advantage is that you generally can say intelligent things about what's going on far from the thin shell boundary, and can develop a good map between the evolution of the configurations of matter within the "inside" metric and that outside. Its disadvantage is that it is a painful amount of work that computers aren't very good at helping with yet.
Which one chooses will depend on some combination of the problem being studied, personal preference, and trade-offs between manual work, accuracy, and comprehensibility. One can combine the perturbation approach with the inside/outside approach, so that one has (dropping the indices) g_{outside} = \eta_{outside} + h_{outside} and g_{inside} = \eta_{inside} + h_{inside} and an Israel-Darmois junction between g_{inside} and g_{outside} which is at least in principle something one can think about functionally (a test object crossing the junction will have its momentum altered; if it just drifts across, the junction is something like a "kink" linking two otherwise smooth geodesics).
Now, back to the standard cosmology, using the "inside"/"outside" approach. In the "outside", the Robertson-Walker vacuum does just what you say:
> [as the] universe [is] empty and expanding, then some light non-interacting test dust would indeed get pulled apart everywhere and at all length scales
and then
> But when we add matter and let it collapse
... we can no longer describe the whole spacetime as a Robertson-Walker vacuum that we're simply probing. Our matter perturbs the RW metric, so the true metric g must be the combination of RW and whatever metric the introduced mass generates. Since the introduced mass is collapsing, it's probably approximately described some sort of Tolman dust. So we select e.g. the Lemaître-Tolman-Bondi metric and perturb that in turn, so g = RW + LTB + O(h) + O(h^2) + ... where h is the deviation of the matter we introduced from exact LTB, and h^n is higher-order perturbations.
This g, like any other metric, will generate a set of geodesics throughout the spacetime. If our collapsing matter was a dust, then each of the individual particles will select its own metric. Flashes of light will find their own null geodesics determined by the metric, and so forth. The results are different from Robertson-Walker, at least close to the matter that we introduced. At even fairly short distances the higher order terms in h fall into irrelevance, and at large distances h itself becomes negligible.
I want to write a bit about extended objects, but that'll have to wait until a bit later. The key feature there is that the individual components of such objects are stuck together, so they don't find their own individual geodesics. The behaviour of extended objects in curved spacetime is fun!
> may eventually reach a kind of equilibrium where the effects of the expansion and the collapsing matter cancel out locally leading to a non-expanding region
You can move the thin shell of the junction around to explore that kind of thing. Or taking a purely perturbational approach you can look for regions of spacetime where geodesics have the characteristics you want (drop down test objects near point in spacetime near a suspected "cancelling-out" region and see if the inertially-moving test objects eventually collide or separate).
How you represent what's happening there is up to you. You could use your colloquial expression, a line from the Karen Masters link or something like it, or talk about the mathematical structures. The theory underdetermines the precise mathematical expression (there are lots of exactly equivalent ways of writing it down, and lots of so close it doesn't matter approximations); the English has it even worse. :/
However, the key thing is that exact vacuum solutions stop being exact when you add matter. Matter rarely cooperates with exact non-vacuum solutions (they typically are maximally symmetric: spherical arrangements of matter of uniform density, with no multipole moments from internal motions, etc. which is far from what we see measuring Earth's or Moon's gravitational fields with e.g. GOCE and GRAIL, and solar prominences big and obvious, galaxies can be spiral, a tiny black hole orbiting a large one will induce tides on the large one, and so on.)
- --
[0] I'm deliberately going to avoid talking about non-geodesic trajectories through spacetime for now. May come back to that when discussing extended objects in a followup.
- --
ETA: I should add that the metric expansion is not really an expansion of space, but rather an extra term tacked onto the metric g. So following the perturbation notation above we could start with flat spacetime \eta and say g = \eta + {expansion} and end up with g = {expanding Robertson-Walker}.
This is quite common: the "true" background is Minkowski flat-space, and then the vacuum eternal spherically-symmetric black hole is a perturbation on that, and then there may be other perturbations. (Instead of starting with the "true" background of Schwarzschild.) This kind of thing leads to some insights about asymptotic flatness.
So for an expanding cosmos, we can start with a perfect flat non-expanding one and perturb that, expanding the metric.
(This is being really loose with terminology, but I think gets the point across reasonably. It's the metric expansion of space because there is a time parameter in the line-element of the metric; spatial distances are determined by when (coordinate-time) you are. g = \eta + f(t), where f is a function taking a time coordinate and generating a perturbation field; maybe better to write g(t) = \eta(t) + f(t) where g(t) is the metric's spatial distance functions for all events on a slice of space that where everything is at the same time coordinate t; but this dives into 3+1 spacetime splittings which is another big can of worms...).
- --
ETA2: the Overview section of this puts it well in a different way. (No claims about the other sections, I haven't read them yet.) https://www.cs.mcgill.ca/~rwest/wikispeedia/wpcd/wp/m/Metric...
