Let $\mathbb{H}^2$ be the upper sheet of the hyperboloid defined by $x^2+y^2-z^2=-1$ in three-dimensional Minkowski space $(\mathbb{R}^3, g_M)$, where $g_M = \text{diag}(1,1,-1)$. In other words, consider the hyperboloid model of the hyperbolic plane.
Given two points $p,q$ on the hyperboloid, I would like to derive a parametric equation for the geodesic connecting $p,q$.
So far, I was able to specify the plane containing $p,q$ and the origin by $\pi : (p \times q) \cdot x=0$, because $(p \times q)$ is a normal vector to the plane and the plane must contain $(0,0,0)$. The equation of the geodesic is now the intersection of this plane with the hyperboloid. Since $(p \times q)$ can't be parallel to the $z$-axis, we may assume w.l.o.g. that the $x$ coefficient of the plane is non-zero. Solving for $x$ using the plane equation and substituting this into the hyperboloid equation gives a complicated equation in $y$ and $z$.
I also thought of using the geodesic characterisation $p\cosh(t)+v\sinh(t)$ for a unit tangent vector $v$, but here the trouble is finding the correct $v$ so that the geodesic will pass through $q$.
Any help is appreciated.
Thanks to Andrew's comment, I have figured it out.
Suppose we have points $p = \begin{pmatrix} \sinh(v)\cos(\theta)\\\sinh(v)\sin(\theta)\\\cosh(v)\end{pmatrix}$ and $q = \begin{pmatrix} \sinh(u)\cos(\phi)\\\sinh(u)\sin(\phi)\\\cosh(u)\end{pmatrix}$.
We first rotate $p$ into the $x$-$z$-plane by applying the isometry
$R_{-\theta} = \begin{pmatrix} \cos(-\theta) & -\sin(-\theta) & 0\\ \sin(-\theta) & \cos(-\theta) & 0 \\ 0 & 0 & 1 \end{pmatrix}$
Next, we apply the isometry
$S_{-v} = \begin{pmatrix} \cosh(-v) & 0 & \sinh(-v)\\ 0 & 1 & 0 \\ \sinh(-v) & 0 & \cosh(-v) \end{pmatrix}$,
which is essentially a shift along the hyperboloid.
Lets call the composite isometry $T^{-1} = S_{-v}R_{-\theta}$.
Observe that this isometry acts on p as $p \mapsto T^{-1}p = \begin{pmatrix} 0\\0\\1\end{pmatrix}$.
We need to figure out its action on $q$, which is easily calculated to be
$q \mapsto T^{-1}q = \begin{pmatrix} \sinh(u)\cosh(-v)\cos(\phi-\theta) + \cosh(u)\sinh(-v) \\ \sinh(u)\sin(\phi-\theta) \\ \sinh(u)\sinh(-v)\cos(\phi-\theta) + \cosh(u)\cosh(-v) \\ \end{pmatrix}$
The tangent plane $T_{(0,0,1)}\mathbb{H}^2$ has orthonormal basis $\{(1,0,0),(0,1,0)\}$. Let $\Pi = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0\end{pmatrix}$ be the projection onto this tangent plane.
Let $v = \frac{\Pi T^{-1}q}{\lvert\lvert \Pi T^{-1}q \rvert\rvert }$, where the norm is taken with respect to the Minkowski metric. Thus, $v$ is a unit vector in $T_{(0,0,1)}\mathbb{H}^2$.
It suffices to find the geodesic starting at $(0,0,1)$ in the direction of $v$. By the formula given in the original question it is given by $\gamma(t) = \begin{pmatrix} 0\\0\\1\end{pmatrix}\cosh(t) + v\sinh(t)$. It is easy to see that this curve is unit-speed. It is also easy to calculate that, in order for the curve to terminate at $T^{-1}q$, we must have $t \in [0, \text{arccosh}(-p\cdot q)]$, where the dot product is with respect to the Minkowski metric, so we have derived the distance between points $p$ and $q$ along the way.
The final answer is given by inverting the isometry:
$T\gamma$ is the cuve we are looking for.