Are there unique geodesics in the NIL and SOL geometry? Is there a unique geodesics between any two points in the NIL (resp. SOL) geometry?
If so, is there a nice way of parametrizing them? For example geodesics in $S^3$ can be parametrized using the embedding in $\mathbb{R}^4$ and $\sin , \cos$ functions. Geodesics in hyperbolic space can be parametrized using the hyperboloid model and the functions sinh,cosh.
 A: Despite their non-uniqueness, a lot is known about geodesics of left-invariant metrics on Heisenberg groups. For example, Jang-Park [Conjugate points on 2-step nilpotent groups, Geom. Dedicata 79 (2000), no. 1, 65--80] describe conjugate points for all geodesic passing through the identity element of a simply-connected 2-step nilpotent Lie group with 1-dimensional center. Earlier Kaplan and Eberlein obtained explicit equations for geodesics (see references in the above paper).
A: The geodesics in SOL that are not contained in a vertical (hyperbolic) foliation tend to come in families meeting periodically and are thus highly non unique. They are precisely described in my paper L’Horizon de SOL.  Expos. Maths 1998.
In this paper, the SOL geodesics are explicitly parametrized using elliptic functions.
A: The geodesics between points are not unique in both cases. Moreover the following is true: if $M$ is a universal cover of a compact Riemannian manifold whose fundamental group is virtually solvable but not virtually abelian, then there are conjugate points on some geodesics in $M$ and hence geodesics between some points are not unique.
See Croke and Schroeder "The fundamental group of compact manifolds without conjugate points", Comment. Math. Helv.  61  (1986),  no. 1, 161--175, MR847526, for the case when the metric is analytic, and Lebedeva "On the fundamental group of a compact space without conjugate points", here, for the general case.
A: There are not unique geodesics between points in Sol geometry. The Sol metric on $\mathbb{R}^3$ may be given as $e^{-z}dx^2+e^zdy^2+dz^2$. The claim is that there is not a unique geodesic between $(0,0,0)$ and $(t,t,0)$ for $t$ large enough. 
There is a rotational isometry $(x,y,z)\overset{\varphi}{\mapsto}(y,x,-z)$. This leaves invariant the line $l=\{(t,t,0)|t\in \mathbb{R}\}$,
and therefore $l$ is a geodesic. For $t$ small enough, $(t,t,0)$ will lie in a normal coordinate patch about $(0,0,0)$, so $l$ will be the unique shortest geodesic between these points.
However, for $t$ large, there are much shorter paths connecting $(0,0,0)$ and $(t,t,0)$. The length of the geodesic $l$ is linear in $t$. But one may take a piecewise geodesic path, starting as a geodesic in the hyperbolic plane $x=0$ connecting $(0,0,0)$ and $(0,t,0)$, and then a geodesic in the hyperbolic plane $y=t$ connecting $(0,t,0)$ and $(t,t,0)$. Since the path $(0,u,0), 0\leq u\leq t$ is a horocycle in the hyperbolic plane $x=0$, and similarly $(u,t,0),0\leq u\leq t$ is a horocycle in the hyperbolic plane $y=t$, 
the length of this piecewise geodesic path is on the order of $C\log(t)$. Since the metric space is complete, there is a minimal length geodesic connecting $(0,0,0)$ and $(t,t,0)$ which is not invariant under $\varphi$, and thus there are at least two minimal length geodesic paths (and at least three geodesics) connecting the two points.  
