Equations defining hyperbolic geodesics in $\mathbb C \setminus\{0,1\}$ Let $X=\mathbb C\setminus\{0,1\}$, equipped with the hyperbolic structure it inherits from Klein's modular $\lambda$ function $\lambda:\mathbb H \to X$. In each (non-peripheral and nontrivial) free-homotopy class of loop $S^1\to X$, there is a unique hyperbolic representative.
For example, the figure-eight curve $\alpha$ with winding number $+1$ around $1$ and $-1$ around $0$, pictured below, has a unique hyperbolic geodesic representative.

In fact, the curve drawn in the above linked picture is precisely the path traced by the hyperbolic geodesic representative of $\alpha$, which is given by the lemniscate with equation $$16(x^2+y^2)((x-1)^2+y^2)=1.$$
Proof sketch: The conformal transformation $w(z)=(2z-1)^2$ is an orbifold covering map from $X$ to a hyperbolic orbifold $Y$, where $Y$ is $\mathbb C \setminus\{1\}$ with a single orbifold point at $0$ with angle $\pi$. The geodesic representative of $\alpha$ projects to a geodesic $\eta$ that begins on the orbifold point, travels once counterclockwise around $1$ and hits the orbifold point, then travels once clockwise around $1$ before finishing again at the orbifold point. The anticonformal automorphism $w\mapsto \frac {\bar{w}}{\bar{w}-1}$ is an orientation-reversing isometry of $Y$ that preserves the homotopy class of $\eta$, so it's not hard to see that $\eta$ must have image given by the fixed-point set of this isometry, namely $|w-1|^2=1$. Putting the pieces together, we find that the points in the image of the geodesic representative of $\alpha$ satisfy $|(2z-1)^2-1|^2=1$, which simplifies to the equation above.
Question: Are all closed hyperbolic geodesics on $X$ algebraic curves?
I admit that it is a bit rash to ask the question in this way — one could more modestly ask which hyperbolic geodesics have polynomial defining equations — but I do want to emphasize that I am not aware of any negative results about the defining equations (the other two figure-eight curves are easy to work out as above). Of course, it's possible that the lemniscate equation arises as some low-dimensional accident due to the special symmetry present in the figure-eight curve, but it also seems entirely plausible that all hyperbolic geodesics on $X$ are somehow very convenient from the viewpoint of the canonical complex coordinate of $X$.
I've tried just looking at other examples using Mathematica to see if I might guess some equation for them (Mathematica's built-in ModularLambda function makes it easy to draw hyperbolic geodesics on $X$), but perhaps I am not fluent enough in solution sets of rational equations in the plane to make anything work.
Two other comments:

*

*I've tagged this question with "modular forms", because one could think of the above polynomial equation as a polynomial identity relating the real and imaginary parts of the modular function $\lambda$ along the hyperbolic geodesics in $\mathbb H$ that project to closed curves.


*There is an obvious analogy with this very recent MathOverflow question, which is a knot theory version of the same idea.
 A: Let me rephrase the proof of @Ian Agol as I understand it.
Let $\Gamma(2)$ be the congruence subgroup
of level $2$ (it consists of all $2\times 2$ integer matrices $A$ with
${\mathrm{det}\, }A=1$ and $A\equiv I\; (\mathrm{mod}\; 2)$.
Each $A=(a_{ij})\in\Gamma(2)$ with $|{\mathrm{tr}\, }A|>2$ defines the geodesic
by the rule: let $\gamma$ be the circle orthogonal
to the real axis, and passing through the fixed points of the
linear-fractional transformation $z\mapsto (a_{11}z+a_{12})(a_{21}z+a_{22})$. Then $\lambda(\gamma)$ is a closed geodesic and
they are all obtained in this way. Here $\lambda:H\to C\backslash\{0,1\}$
is the universal covering automorphic with respect to $\Gamma(2)$. Here $H$ is the upper half-plane.
Reflection in $\gamma$ is given by the matrix
$B=2A-({\mathrm{tr}\, }A)I$:
$$z\mapsto \sigma(z):=
\frac{b_{11}\overline{z}+b_{12}}{b_{21}\overline{z}+b_{22}},\quad B=(b_{ij}).$$
Consider the group
$$\Gamma'=B^{-1}\Gamma(2)B=d^{-1}B'\Gamma(2)B,$$
where $d={\mathrm{det}\, }B$ and $B'B=dI$.
Claim. $\Gamma_0:=\Gamma'\cap\Gamma(2)$ is of finite index in both
$\Gamma'$ and $\Gamma(2)$.
Proof. If $C\in\Gamma(2d)=\{ A\in\Gamma(2):A\equiv I\; (\mathrm{mod}\; 2d)\}$,  then $E=d^{-1}BCB'$ is in $\Gamma(2)$, so $C=B^{-1}EB\in\Gamma'$. So both groups
contain $\Gamma(2d)$ with finite index.
Now let $\lambda$ be the modular function, $\lambda:H\to C\backslash\{0,1\}$,
is it invariant under $\Gamma(2)$, and let
$\mu(z)=\overline{\lambda(\sigma(z))}$, which is invariant with respect to $\Gamma'$. Let
$X:=H/\Gamma_0$. Then $X$ is a
Riemann surface of finite type whose only boundary elements are punctures,
and both $\lambda,\mu$ are meromorphic functions on it.
It is easy to see that their singularities at the punctures are at most
poles, therefore they are algebraically dependent. So we have
a polynomial identity $F(\lambda,\mu)=0$ which holds in $H$.
Now notice that for $z\in\gamma$ we have $\lambda(z)=\overline{\mu(z)}$
therefore $(\lambda+\mu)/2={\mathrm{Re}\, }\lambda(z)$
and $(\lambda+\mu)/(2i)={\mathrm{Im}\, }\lambda(z),\; z\in\gamma$.
By the uniqueness theorem this implies that there is a polynomial equation
$G({\mathrm{Re}\, }\lambda,{\mathrm{Im}\, }\lambda)=0$ which holds for $z\in\gamma$.
Remark. Except in the two trivial examples,
when $\Gamma'\cap\Gamma(2)\in SL(2,Z)$,
the degree of this polynomial seems to be
very large, so it does not help in plotting the geodesics, while they can
be easily plotted directly by using an explicit expression for $\lambda$.
A: I wondered this myself, I made some similar pictures to approximate the stable lamination of a pseudo-Anosov map (the blue curve is a geodesic, the other colors horocycles).

Every geodesic on the modular orbifold (or $\mathbb{H}^2/\Gamma(2)$ that you're considering) is fixed by a hyperbolic element $A\in \Gamma(2) < PSL_2(\mathbb{Z}), |tr(A) |> 2$ (note: I’m writing matrix representatives in $GL$ of representatives in $PGL$).
Then the matrix $A'=2A-tr(A)I \in PGL_2(\mathbb{Q})$ is a matrix in the commensurator fixing the geodesic fixed by $A$ ($tr(A')=0, det(A')<0$).
One interesting analogy to make is that if $B\in PGL_2(\mathbb{Q})$ is an elliptic involution, so $tr(B)=0, det(B)>0$, then the fixed points $b,\overline{b} \in \mathbb{C}-\mathbb{R}$ of $B$ in the upper and lower half hyperbolic planes are quadratic imaginary points. Hence $\lambda(b)$ is an algebraic number by the Kronecker Jungendtraum, lying in an abelian extension of the quadratic imaginary number field generated by that point (this can be seen by the relation between $\lambda$ and $j$).
I think this extends to your case. Consider $\Gamma(2)\cap A'\Gamma(2)(A')^{-1}$, then this is a finite-index subgroup of $\Gamma(2)$. In turn this gives a finite étale congruence cover which is an algebraic curve $C$ (see the Belyi Theorem). Then $A'$ lifts to an anti-holomorphic involution of $C$, and hence the fixed set should be a real-algebraic subset of this curve with a component projecting to the geodesic fixed by $A$. Projecting down to $\mathbb{H}^2/\Gamma(2) = \mathbb{C}-\{0,1\}$, I think the image downstairs should also be a component of a real algebraic curve containing the desired geodesic. However, it seems to me that a priori there could be other components (as immersed curves).
