If by "cosmological constant dominated universe" you mean $Ω_Λ=1$, i.e., de Sitter space, then this is pretty easy to solve. You can use de Sitter static coordinates,
$$ds^2 = (1-r^2/R^2)\,dt^2 - (1-r^2/R^2)^{-1}dr^2$$
in which there is a natural notion of constant distance, namely, both galaxies being at some constant $r$ independent of $t$. (This agrees with radar distance, and you can even maintain it with a rigidly moving strut.) Rather than solve the equations of motion, you can use the fact that $(dr/ds)^2-(r/R)^2$ is a constant of motion, with the first term being in effect kinetic energy and the second gravitational potential. From this you can easily find the "escape velocity" from one galaxy to the other, given their positions (if both are in line with the origin).
If you mean something like our present-day parameters, with $Ω_m \gg 0$, then the problem is much more complicated, and I don't think it's even well defined.
Your second paper (Clavering) correctly points out that the notion of constant distance in the first paper (Davis et al) is dubious since it doesn't reduce to constant distance in Minkowski space. Radar distance does avoid that problem, but it has another problem: it isn't reciprocal, i.e., if B is at a constant radar distance from A, then A isn't generally at a constant radar distance from B. At least, it doesn't appear to be reciprocal. Directly verifying that is difficult because the worldlines obtained from this definition are so complicated, which is another problem with it.
Also, while radar distance avoids the issue of blueshifted receding galaxies in Minkowski space, it doesn't avoid it in general. Even in de Sitter space, galaxies stationary at different $r$ in static coordinates will see each other as red/blueshifted, and if you give them small nonzero speeds you can create situations in which they're receding in radar distance but blueshifted, or vice versa. Zero-redshift distance avoids that by construction in one direction, but isn't reciprocal and doesn't match the natural notion of distance in de Sitter space.
The root problem is unfortunately not solvable: there is no good definition of constant distance in FLRW spacetimes.
My inclination is to solve the problem in a toy Newtonian model, because the GR "corrections" to the Newtonian result are not obviously even meaningful.
In the Newtonian model, there is a uniform sphere of matter centered at the origin with a radius of $a(t)R$, where $R$ is large enough that the sphere includes the whole experiment. The universality of gravity immediately implies that freefall motion of a test particle satisfies $x''(t)/x(t) = a''(t)/a(t)$. This is the same as Davis et al's equation (14). Newtonian distance also matches the distance measure used in Davis et al. (That's not to say I think they were right to use it in a general-relativistic context.) You should be able to plug in a particular $a(t)$ and solve the problem easily enough.
Miscellaneous comments:
Both galaxies are very likely in accelerating reference frames.
Note that because of the background Hubble flow, there isn't any relativity of motion here. There are distinct cases where one galaxy has zero peculiar velocity and the other accelerates to stay at a fixed distance (however defined), and where both have a peculiar velocity and both accelerate, symmetrically or not. You can see this in the de Sitter case.
(Also, "in accelerating reference frames" is meaningless; just say that the galaxies are accelerating or noninertial.)
This can be understood by the fact that the peculiar velocities of both galaxies would normally trend towards zero, and so expansion would pull them apart. Therefore, a force must act on them to keep them together.
Peculiar velocities do tend to zero at late times in an expanding universe, but that doesn't mean that the galaxies will tend to move apart immediately given these initial conditions. Davis et al correctly point out that if $Ω_m \gg 0$ then the distance between them is likely to decrease at first if they move geodesically. They call that counter-intuitive in the abstract, but the reason for it is just, well, gravity. There's matter between the galaxies, and it attracts them, so they accelerate inward. To counter that you need a non-gravitational outward acceleration.
Note that peculiar velocities tending to zero is not a space-expansion effect. It happens even in the Newtonian model, essentially because an object with a peculiar velocity will move somewhere else, "stopping" only when its velocity matches the background matter's.
It's tempting to speculate that $v_1=v_2$ because of conservation of energy, but in an expanding universe governed by GR conservation of energy is not really well defined.
Energy conservation is violated even in the Newtonian model, because it treats the Hubble flow as a fixed background. It's true that defining a conserved energy is problematic in GR, but you could do much better than Davis et al by including the backreaction on the FLRW matter, at the cost of making the analysis enormously more complicated.
