For some reason, I would need to know what are the timelike geodesic congruences which cover the entire exterior region of the Schwarschild spacetime. In fact the only thing I really need is the geodesic vector fields defining these congruences.

I have tried something, but being an absolute non-expert in this field I'm not quite sure about it. Here is my attempt:

The derivative with respect to an affine parameter of a generic timelike geodesic is

$$\left\{\matrix{\dot t&=&(1-{r_s\over r})^{-1}E\cr \dot r&=&\pm \left(E^2-(1-{r_s\over r})(1+{L^2+Q\over r^2})\right)^{1/2}\cr \dot\theta&=&\pm\left(Q-L^2\mathrm{cotan}^{2}\theta\right)^{1/2}r^{-2}\cr \dot\varphi&=&{L\over r^2\sin^2\theta} }\right.$$

where $E$, $Q$, $L$ are constants of the motion. Seeing the RHS as a vector field, we obtain a geodesic congruence of each value of the triple of constants. However, we must have $L=0$ for the field to be defined everywhere in the outside region. Hence the solution to our problem is given by the familly of vector fields above with $L=0$ and $E,Q$ arbitrary.

Is this correct or am I missing something ?

Aside: The texts I have found concerning the Schwarzschild geodesics always focus on the equatorial ones for which $\theta$ is constantly equal to $\pi/2$. I found the formula above in the Wikipedia entry for the Kerr metric, and I set $a$ to $0$. I would appreciate any reference for this formula in the context of the Schwarzschild metric.

Edit: I realize that the situation is more complex than I thought. However, I think I have an argument proving that in 1+2 dimensions, the only everywhere-defined geodesic vector fields become nearly radial at infinity.

Here it is: Point 1: Since the geodesic vector fields are not explicitly $t$-dependent, their intergral curves in $t=$constant planes form a partition.

Point 2: The partitions into straight lines of ${\mathbb R}^2\setminus B$, where $B$ is a ball are either radial or made of all parallel lines (this goes astray in $3$ space dimensions).

Point 3: A timelike geodesic congruence of the Schwarzschild spacetime will go to a partition of the above kind in the limit of small $r_S/r$. However, the parallel case is excluded, since in a family of geodesics which are nearly parallel straight lines at infinity we can always find two which are symmetrically deviated by the mass at the origin and will therefore cross somewhere.

Can this argument be pushed to prove that the congruence is exactly radial ? What of the $1+3$-dimensional case ?