By a Blaschke product, I guess we mean a function of the form
$$B(z) = \beta \prod_{k=1}^n\frac{z-a_k}{1-\overline{a_k}z},$$
where $|\beta|=1$ and $|a_k|<1$ for each $k$. These are exactly the $n$-to-1 functions on the unit disk with specified zeros in the disk. We can define one using Mathematica in terms of a $\beta$ and a list of the $a$s as follows:
term[a_] := (z - a)/(1 - Conjugate[a] z);
B[As_List, beta_][z_] := beta*Times @@ (term /@ As);
Now, a given point on the unit circle has $n$ preimages that are also on the unit circle. The paper states that all the $n$ choose $2$ hyperbolic geodesics joining the pairs of these preimages are tangent to a particular algebraic curve in the disk. The curve may be visualized as the envelope of a large number of these geodesics.
To implement this, we'll need a simple way to generate the geodesics. Given a pair of points on the unit circle, the geodesic joining them is simply the unique circle that passes through both of them and is perpendicular to the unit circle at those points. When the two points are opposite one another, that circle will actually be a line through the origin.
toR2[z_] := {Re[z], Im[z]};
arc[{z1_, z2_}, z0_] := Module[
{theta1, theta2}, theta1 = Arg[z1 - z0];
theta2 = Arg[z2 - z0];
If[Abs[N[theta1 - theta2]] > N[Pi],
If[N[theta1] < N[theta2],
theta1 = N[theta1 + 2 Pi],
theta2 = N[theta2 + 2 Pi]]
];
Circle[toR2[z0], Abs[z1 - z0],
Sort[{theta1, theta2}, N[#1] < N[#2] &]]];
geodesic[{z1_, z2_}] := Module[{center},
If[Arg[z2] - Arg[z1] == Pi,
Line[toR2 /@ {z1, z2}],
center = Exp[I (Arg[z2] + Arg[z1])/2]/Cos[(Arg[z2] - Arg[z1])/2];
arc[{z1, z2}, center]]];
Using this, we could now write a function pic[As_,beta_,m_]
to generate the picture using m
initial points on the circle as follows:
preimages[As_List, beta_][lambda_] :=
z /. NSolve[B[As, beta][z] == lambda, z];
pic[As_List, beta_, m_] :=
Module[{allPreimages, preimagePairs, geodesics},
allPreimages = Table[preimages[As, beta][Exp[2 k*Pi*I/m]], {k, 1, m}];
preimagePairs = Flatten[Subsets[#, {2}] & /@ allPreimages, 1];
geodesics = geodesic /@ preimagePairs;
Graphics[{{Thickness[0.01], Circle[]}, {Opacity[0.6], geodesics}}]]
When $n=3$, for example, we generate an object called a hyperbolic ellipse:
pic[{1/2, I/2, (1 + I)/3}, 1, 12]

Here's an example for the $n=4$ case:
pic[{1/2, I/2, (1 + I)/3, -1/2 - I/3}, 1, 12]

It's odd that the paper you cite has no pictures. The author has a presentation here with a couple of images on page 5 that seem to agree with the images above.