Visualization of Riemann Surfaces of Algebraic Functions

The author brings the symbolic, numeric and graphic capabilities of Mathematica to bear on Riemann surfaces. The result is an illuminating catalog ranging from simple to extremely complicated images.

by Michael Trott

We begin by considering an analytic function z (w) such z ((Subscript[w,  1])) = z ((Subscript[w,  2])) for (Subscript[w,  1]) ≠ (Subscript[w,  2]). (Think of z (w) = (w^( 2)) as a very simple example where z (1) = z (-1).) The inverse function w (z) must be a multivalued function. Properly glueing together copies of the complex z-plane for all values of (Subscript[w,  i]) (z) at points where we have degenerate solutions, we get the Riemann surface of the function. Riemann surfaces of analytic functions are very interesting and useful and have been investigated for more than 100 years. In many mathematics departments one finds plaster or wood models of them (for a nice collection of photographs of such surfaces see [Fischer 1986]). With Mathematica's symbolic and numeric capabilities, it is possible to recreate Riemann surfaces that are sketched in most textbooks on the theory of functions of one complex variable and related books [Brill 1925; Jung 1923; Mathews and Howell 1977; Monna 1975; Osgood 1923] as well as to visualize much more complicated ones. In this article we will focus on the visualization of Riemann surfaces of algebraic functions. (For the more mathematical side of Riemann surfaces of algebraic functions, such as the genus of a Riemann surface or the Riemann-Roch theorem see [Rauch and Lebowitz 1973]. For information concerning the visualization near branch points see [Shihara and Sasaki 1996].) Given an (irreducible) polynomial p in two variables \(TraditionalForm
      \`p = ∑\_\(j, \(k\ \) \( = \ \) 0\)\%n a\_\(i  j\)\ \(z\^\(\ j\)\) 
        w\^\(\ k\)\), the equation p = 0 defines w (z) implicitly as a generically multivalued function. In fact this class of functions exhausts all algebraic functions, including elliptic and hyperelliptic ones. Graphing the real or imaginary part of w (z) creates a faithful representation of the Riemann surface of w (z). 
The semi-numerical algorithm for the construction of the Riemann surface given below will be used to visualize the Riemann surfaces of a variety of algebraic functions. 
As a warm-up before dealing with the general case, let us look at some simple Riemann surfaces.

Some simple examples

Lets us start with some simple examples of the form w = f (z) where f is a function with branch points and branch cuts that are easy to find and easy to exclude in graphic output. 
Because we are here interested in the visualization of the Riemann surfaces, from time to time in the following discussion we will explicitly exclude branch points, branch cuts and poles and start an ε(ε > 0) away from such points with the polygons of the surface. 
Our first example will be the simplest possible of all, that is (w^2) - z = 0 or w (z) = (Sqrt[z]). The only branch point is z = 0, and the branch cut is usually chosen to extend along the negative real axis. (Here and in the following we refer to branch points in the finite part of the complex z-plane. Possible branch points at ∞ will be discussed in the last section.)

Plot3D[Im[(Sqrt[x + I y])], {x, -2, 2}, {y, -2, 2}, 
      PlotPoints -> 30];

We choose the imaginary part of (Sqrt[z]) as a representation of the Riemann surface. Rewriting z in polar form as r  (e^i)φ we have the parametric representation \(TraditionalForm
      \`{r, φ} -> {r\ cos \((φ)\), r\ sin 
        \((φ)\), \@r\ sin \((φ\/2)\)}\):

ParametricPlot3D[{r Cos[φ], r Sin[φ], 
        (Sqrt[r]) Sin[(φ/2)], 
        SurfaceColor[Hue[(φ/(4 π))]]}, {r, 0, 1}, {
        φ, 0, 4 π}, PlotPoints -> {20, 60}, 
      Boxed -> False, Axes -> False];

Looking at the cross section y = 0 we see the typical branching behaviour of a second order branch point.

ParametricPlot[
      Evaluate[({r Cos[#1], (Sqrt[r]) Sin[#(1/2)]}&)/@{0, π, 2 π, 
            3 π, 4 π}], {r, 0, 1}, Frame -> True, 
      Axes -> False];

A more algebraic way to generate the surface of the imaginary part would be the following: 
We write z = x + i  y, w = u + i  v, where x, y, u and v are purely real quantities. 
Then we solve the pair of equations 
            \(TraditionalForm
      \`ℜ(z) = 
        \(ℜ(w\^2) = \(x = u\^2 - v\^2\)\)\)
            \(TraditionalForm
      \`ℑ(z) = 
        \(ℑ(w\^2) = \(y = 2  u\ v\)\)\)
for v, eliminate u and then display v = v(x, y). (We will come back to this approach later.)

Simplify[Solve[{x == (u^2) - (v^2), y == 2 u v}, {v}, {u}]]

The two real solutions for v for real x and y then yield the same surface as above. (The real part u (x, y) is used for the coloring.)

\(\(Show[
      \((Plot3D[{#1, Hue[Abs[y\/#1]]}, {x, \(-2\), 2}, {y, \(-2\), 2}, 
              DisplayFunction -> Identity, PlotPoints -> 38]&)\)/@{
          \@\(\@\(x\^2 + y\^2\) - x\)\/\@2, 
          \(-\(\@\(\@\(x\^2 + y\^2\) - x\)\/\@2\)\)}, 
      DisplayFunction -> $DisplayFunction, BoxRatios -> {1, 1, 1}]; 
    \)\)

Our second example is (w^3) - z = 0 or w(z) = (z^(1/3)). The only branch point is z = 0 and the branch cut is typically chosen again to extend along the negative real axis. We use again the imaginary part as a faithful representation of the Riemann surface. The following graphic shows the principle sheet—the branch cut along the negative real axis is clearly visible.

Plot3D[Im[(((x + I y))^(1/3))], {x, -2, 2}, {y, -2, 2}, 
      PlotPoints -> 30];

Rewriting z again in polar form as r  (e^i)φ we have the parametric representation \(TraditionalForm
      \`{r, φ} 
        -> {r\ \(cos(φ)\), r\ \(sin(φ)\), \) (r^(1/3))sin((φ/3))}:

ParametricPlot3D[{r*Cos[\CurlyPhi]], r*Sin[φ], r^(1/3)*Sin[φ/3],
    SurfaceColor[Hue[\CurlyPhi]/(6*Pi)]]}, {r, 0, 1}, {φ, 0, 6*Pi},
   PlotPoints -> {20, 90}, Boxed -> False, Axes -> False];

In a similar way one can treat (w^( 6)) - z, or equivalently w(z) = (z^(1/6)). Now we have a Riemann surface with six sheets.

ParametricPlot3D[{r Cos[φ], r Sin[φ], 
        (r^(1/6)) Sin[(φ/6)], {
          SurfaceColor[RGBColor[(φ/(12 π)), 0, 0], 
            RGBColor[(φ/(12 π)), 0, 0], 1.45], 
          EdgeForm[{Thickness[0.0001], 
              RGBColor[0, 1 - (φ/(12 π)), 0]}]}}, {r, 0, 
        0.6}, {φ, 0, 12 π}, PlotPoints -> {20, 180}, 
      Boxed -> False, Axes -> False];

The form of the Riemann surface for w(z) = (z^(1/n)) is now obvious. Next let us look at the function (w^2) + (z^2) - 1 = 0 which can be solved with respect to w in the form w(z) = (Sqrt[1 - (z^2)]). This function has two branch points at ±1. In the literature the branch cut is usually chosen to extend along the real axis, starting from the branch points and going to ±∞. Here the real and imaginary parts of the principle sheet are shown:

Show[
      GraphicsArray[
        (Plot3D[#1[(Sqrt[1 - (((x + I y))^2)])], {x, -2, 2}, {y, 
                  -2, 2}, PlotPoints -> 37, 
                DisplayFunction -> Identity]&)/@{Re, Im}]];

Here both sheets of the real and of the imaginary parts are shown:

With[{ε = (1/(10^12))}, 
      Show[GraphicsArray[
          Function[reIm, 
              Show[Function[f, 
                    Apply[Plot3D[reIm[f], {x, -2, 2}, {y, ##1}, 
                          PlotPoints -> {37, 27}, 
                          DisplayFunction -> Identity]&, {{
                          ε, 2}, {-ε, 
                          -2}}, {1}]]/@{(Sqrt[1 - (((x + I y))^2)]), 
                    -(Sqrt[1 - (((x + I y))^2)])}, 
                BoxRatios -> {1, 1, 1}, 
                ViewPoint -> {4, -2, 0.7}]]/@{Re, Im}]]];

The next function we will treat is (w^2) + (z^4) - 1 = 0 which can be solved with respect to w in the form w (z) = (Sqrt[1 - (z^4)]). This function has four branch points at ±1 and ±i, and the branch cuts start at the branch points and extend radially outwards. Here the real and imaginary parts of the principle sheet are shown:

Show[
      GraphicsArray[
        (Plot3D[#1[(Sqrt[1 - (((x + I y))^2)])], {x, -2, 2}, {y, 
                  -2, 2}, PlotPoints -> 37, 
                DisplayFunction -> Identity]&)/@{Re, Im}]];

The following picture shows all four sheets of the imaginary part (we omit the polygons from the third quadrant to have a better look at the inner parts) of the function (w^4) = 1 - (z^4).

rsf = 
      With[{ε = (1/(10^12))}, 
        Show[Function[f, 
              Apply[Plot3D[Re[f], {x, #1, #2}, {y, #3, #4}, 
                    PlotPoints -> 17, 
                    DisplayFunction -> Identity]&, {{ε, 2, 
                    ε, 2}, {-ε, -2, 
                    ε, 2}, {-ε, -2, 
                    -ε, -2}}, {1}]]/@{
              (Sqrt[1 - (((x + I y))^4)]), -(Sqrt[1 - (((x + I y))^4)]), 
              I (Sqrt[1 - (((x + I y))^4)]), 
              -I (Sqrt[1 - (((x + I y))^4)])}, 
          DisplayFunction -> $DisplayFunction, 
          BoxRatios -> {1, 1, 1}, ViewPoint -> {4, -2, 1.7}]];

For Riemann surfaces with more than two or three sheets we typically encounter the problem that it is hard to look 'inside'. To have a better look in the inner parts of the surface we define a function holePolygon which cuts a hole of 'size' f in a polygon.

holePolygon[Polygon[l_], f_] := 
      Module[{m = Plus@@(l/Length)[l], innerPoints}, 
        innerPoints = (m + f (#1 - m)&)/@l; {
          MapThread[
            Polygon[Join[#1, Reverse[#2]]]&, {
              Partition[(Append[#1, First[#1]]&)[l], 2, 1], 
              Partition[(Append[#1, First[#1]]&)[innerPoints], 2, 1]}], 
          (Line[Append[#1, First[#1]]]&)[innerPoints]}]

The following picture shows the previous Riemann surface with 'holed polygons'.

Show[
      Graphics3D[{EdgeForm[], Thickness[0.0001], 
          SurfaceColor[Hue[0.8199999999999999], Hue[0.22], 
            1.209999999999999], 
          Cases[rsf, _Polygon, ∞] /. 
            p_Polygon :> holePolygon[p, 0.8]}], 
      BoxRatios -> {1, 1, 1.399999999999999}, 
      ViewPoint -> {4, -2, 1.7}];

A slightly more complicated example is given by the function (w^2) - (z^2) + 1 = 0, which can be solved with respect to w in the form w (z) = (Sqrt[(z^2) - 1]). Now we have two branch points at ±1. In the literature the branch cut is usually chosen to extend along the real axis between the two branch points.

Plot3D[
      Im[(Sqrt[(((x + I y))^2) - 1])], {x, -2, 2}, {y, -2, 2}, 
      PlotPoints -> 37];

Using (Sqrt[(z^2)]) - 1 and the corresponding roots of the other sheets given by \(TraditionalForm
      \`\(e\^\(j\ 2/5\ i\ π\)\) 
        \(\@z\^2 - 1\)\), j = 1, 2, 3, 4 strictly inside the first, second, third and fourth quadrant we get pieces of the Riemann which have no discontinuities and which fit together smoothly.

rsf = 
      With[{ε = (1/(10^12)), ϵ = -1}, 
        Show[Graphics3D[{EdgeForm[], Thickness[0.0001], 
              SurfaceColor[Hue[0.8100000000000001], Hue[0.62], 1.83], 
              Cases[Apply[
                    Table[ParametricPlot3D[{r Cos[φ], 
                            r Sin[φ], 
                            Im[(ϵ^j) 
                                (Sqrt[(((r Exp[I φ]))^2) - 1])]}, {
                            r, ε, 2}, {φ, ##1}, 
                          PlotPoints -> {33, 16}, 
                          DisplayFunction -> Identity], {j, 0, 1}]&, {{
                        ε, (π/2) - ε}, {
                        (π/2) + ε, 
                        π - ε}, {π + ε, 
                        ((3 π)/2) - ε}, {
                        ((3 π)/2) + ε, 
                        2 π - ε}}, {1}], _Polygon, 
                  ∞] /. 
                p_Polygon :> holePolygon[p, 0.74]}]]];

Looking at the cross section y = 0 shows the two branch points nicely.

Show[rsf, PlotRange -> {All, {0, 2}, All}, ViewPoint -> {0.3, -1.8, 0.5}];

As the next example of this introductory section let us take ((((w^( 2)) - 1))^2) - z = 0 or w(z) = (Sqrt[1 + (Sqrt[(z^(-1))])]). Now we have a singularity (in the sense that  | w(z) | (Underscript[⟶, z -> 0])∞) at z = 0. Proceeding in a manner similar to the last example gives the following picture of the Riemann surface.Until now we have looked only at the imaginary part of the functions this time we show both the real and the imaginary parts.

g[reIm_] := 
      With[{ε = (1/(10^8))}, 
        Graphics3D[{EdgeForm[], Thickness[0.0001], 
            SurfaceColor[Hue[0.28], Hue[0.36], 1.45], 
            Cases[Function[f, 
                    Apply[ParametricPlot3D[{r Cos[φ], 
                            r Sin[φ], reIm[f]}, {r, 
                            ε, 1}, {φ, ##1}, 
                          PlotPoints -> {30, 20}, 
                          DisplayFunction -> Identity]&, {{
                          ε, π - ε}, {
                          π + ε, 
                          2 π - ε}}, {1}]]/@{
                    (Sqrt[1 + 
                        (Sqrt[(1/((r/(E
                                    ^(-(-(-I φ))
                                        )))))])]), 
                    -(Sqrt[1 + 
                          (Sqrt[(1/((r/(E
                                      ^(-(-(-I φ))
                                        )))))])]), 
                    (Sqrt[1 - 
                        (Sqrt[(1/((r/(E
                                    ^(-(-(-I φ))
                                        )))))])]), 
                    -(Sqrt[1 - 
                          (Sqrt[(1/((r/(E
                                      ^(-(-(-I φ))
                                        )))))])])}, _Polygon, 
                ∞] /. 
              p_Polygon :> holePolygon[p, 0.7599999999999999]}, 
          PlotRange -> {-2.2, 2.5}]];

Show[GraphicsArray[{g[Im], g[Re]}]];

Sometimes expressions one expects to simplify in Mathematica don't; for example, z - (Sqrt[(z^2)]) does not `simplify' to 0. The reason is again the branch cut of the Power function. The following picture demonstrates the difference between the slightly more complicated example (Sqrt[(z^5)]) and ((((Sqrt[z])))^5) in Mathematica.

With[{z = r Exp[I φ]}, 
      Show[GraphicsArray[
          (ParametricPlot3D[{r Cos[φ], r Sin[φ], 
                    Im[#1]}, {r, 0, 3}, {φ, 0, 2 π}, 
                  PlotPoints -> {20, 60}, BoxRatios -> {1, 1, 0.4}, 
                  DisplayFunction -> Identity]&)/@{(Sqrt[(z^5)]), 
              ((((Sqrt[z])))^5)}]]];

Both are two-valued functions. Here both sheets are shown—the two functions represent the same Riemann surface.

With[{z = r Exp[I φ], ε = (1/(10^12))}, 
      Show[GraphicsArray[
          (Graphics3D[
                  ({#1, Map[{1, 1, -1} #1&, #1, {-2}]}&)[
                    Cases[Table[
                        ParametricPlot3D[{r Cos[φ], 
                            r Sin[φ], Im[#1]}, {r, 0, 3}, {
                            φ, 
                            π + ((i 2 π)/5) + ε, 
                            π + (1/5) (i + 1) 2 π - 
                              ε}, PlotPoints -> {20, 15}, 
                          BoxRatios -> {1, 1, 0.4}, 
                          DisplayFunction -> Identity], {i, 0, 
                          4}], _Polygon, ∞]], 
                  BoxRatios -> {1, 1, 0.4}]&)/@{(Sqrt[(z^5)]), 
              ((((Sqrt[z])))^5)}]]];

Now we will take a look at a function with two branch points and one pole: The two branch points are located at -1 and i and the pole is at (3/2) (i - 1). The polynomial form \(TraditionalForm
      \`∑\_\(j, k = 0\)\%n a\_\(i  j\)\ \(z\^j\) w\^k = 0\) is this time much more complicated than in the earlier examples, we use Mathematica's GroebnerBasis to calculate it.

𝓅 = 
      With[{p = ((3 (I - 1))/2)}, 
        Factor[GroebnerBasis[
              Numerator[
                Together[
                  Apply[Subtract, {z + 1 == (aux1^2), I z + 1 == (aux2^3), 
                      w == aux1 + aux2 + (1/(z - p))}, {1}]]], {w, z}, {
                aux1, aux2}, MonomialOrder -> EliminationOrder]
              [[1]]]]

Here is the principal sheet of \(TraditionalForm
      \`\@\(z + 1\) + \@\(i\ z + 1\) + 
        \((z - \(3\/2\) \((i - 1)\))\)\^\(-1\)\).

Plot3D[
      Im[(Sqrt[x + I y + 1]) + (((I (x + I y) + 1))^(1/3)) + 
          (1/(x + I y - ((3 (I - 1))/2)))], {x, -3, 3}, {y, 
        -3, 3}, PlotPoints -> 60];

The implicit representation is a polynomial of order 6 in w,

Exponent[𝓅, w]

so we have now six sheets which we can calculate by using all possible combinations of the two sheets of (Sqrt[z + 1]) and the three sheets of (Sqrt[i  z + 1]).

rsf = 
      With[{ε = (1/(10^12))}, 
        Show[Apply[
            Function[{ϵ1, ϵ2}, 
              Apply[Plot3D[
                    Evaluate[
                      Im[ϵ1 (Sqrt[x + I y + 1]) + 
                          ϵ2 (((I (x + I y) + 1))^(1/3)) + 
                          (1/(x + I y - ((3 (I - 1))/2)))]], {
                      x, #1, #2}, {y, #3, #4}, PlotPoints -> 30, 
                    DisplayFunction -> Identity]&, {{ε, 3, 
                    ε, 3}, {-ε, -3, 
                    ε, 3}, {ε, 3, 
                    -ε, -3}, {-ε, 
                    -3, -ε, -3}}, {1}]], 
            Flatten[Outer[
                List, {-1, 1}, {1, (((-1))^(1/3)), 
                  (((-1))^(2/3))}], 1], {1}], 
          DisplayFunction -> $DisplayFunction, 
          PlotRange -> {-6, 6}]];

Here is the image with the polygons from the second quadrant taken out. We see the two branch points at -1 and i.

rsf2 = 
      Show[Graphics3D[{EdgeForm[Thickness[0.0001]], 
            SurfaceColor[Hue[0.09], Hue[0.42], 1.5], 
            Select[Cases[rsf, _Polygon, ∞], 
              (Max[#1] < 0 && Min[#2] > 0&)@@
                  Transpose[#1[[1
                      ]]]&]}], 
        BoxRatios -> {1, 1, (3/2)}, Boxed -> False, 
        ViewPoint -> {2.299999999999999, -2, 0.9}];

Here are the same surfaces from another viewpoint with some polygons removed:

Show[
      Graphics3D[{EdgeForm[Thickness[0.0001]], 
          SurfaceColor[Hue[0.33], Hue[0.66], 1.11], 
          Select[Cases[rsf2, _Polygon, ∞], 
            Function[p, 
              ((((-1)
                          )^(Round[
                            (1/3) #1[[1
                                ]] 29] + 
                          Round[-(1/3) #1[[2
                                ]] 29])) == 1&)[
                (1/4) Plus@@p[[1]]]]]}], 
      BoxRatios -> {1, 1, 2}, Boxed -> False, 
      ViewPoint -> {2.299999999999999, -2, 2.9}];

The function (z^2) (w^2) - 2  w - 1 = 0 has a pole at z = 0 on one sheet, but not on the other.

sol = Solve[(z^2) (w^2) - 2 w - 1 == 0, w]

(Series[#1, {z, 0, 2}]&)/@(w /. sol)

Here the two sheets are shown:

Show[
      GraphicsArray[
        (Plot3D[Im[#1], {x, -3, 2}, {y, -2, 2}, 
                PlotPoints -> 50, DisplayFunction -> Identity]&)/@
          (w /. sol /. z -> x + I y)]];

Here the two-sheeted Riemann surface of (z^2) (w^( 2)) - 2  w - 1 = 0 is shown:

Show[
      Graphics3D[{EdgeForm[], Thickness[0.0001], Hue[0.7], 
          SurfaceColor[Hue[0.22], Hue[0.21], 2.9], 
          Cases[With[{ε = (1/(10^12))}, 
                  Apply[Function[f, 
                          Graphics3D[
                            Plot3D[Im[f], ##1, PlotPoints -> 25, 
                              DisplayFunction -> Identity]]]/@
                        (w /. sol /. z -> x + I y)&, {{{x, 
                          -2, -ε}, {y, -2, 
                          -ε}}, {{x, ε, 2}, {
                          y, -2, -ε}}, {{x, -2, 
                          -ε}, {y, ε, 2}}, {{
                          x, ε, 2}, {y, ε, 2}}}, {
                      1}]], _Polygon, ∞] //. 
              x_?(Abs[#1] > 5&) :> 5 Sign[x] /. 
            p_Polygon :> holePolygon[p, 0.65]}, 
        PlotRange -> {-3, 3}]];

For one of our final examples in this section, we take the sum of three square roots with branch points at ±i and 0.

Plot3D[
      Im[(Sqrt[-I + x + I y]) + (Sqrt[x + I y]) + (Sqrt[I + x + I y])], {x, 
        -2, 2}, {y, -2, 2}, PlotPoints -> 30];

The implicit representation is again more complicated than the explicit one.

𝓅 = 
      Factor[GroebnerBasis[
            Apply[Subtract, {z - I == (aux1^2), z == (aux2^2), 
                z + I == (aux3^2), w == aux1 + aux2 + aux3}, {1}], {w, z}, {
              aux1, aux2, aux3}, MonomialOrder -> EliminationOrder]
            [[1]]]

Using both sheets for each of the square roots, we get all eight sheets of the Riemann surface.

With[{ε = (1/(10^12))}, 
      Show[Graphics3D[{EdgeForm[], Thickness[0.0001], 
            SurfaceColor[Hue[0.85], Hue[0.15], 1.31], 
            Cases[Show[
                  Apply[Function[{ϵ1, ϵ2, ϵ3}, 
                      Apply[Plot3D[
                            Evaluate[
                              Im[ϵ1 (Sqrt[-I + x + I y]) + 
                                  ϵ2 (Sqrt[x + I y]) + 
                                  ϵ3 (Sqrt[I + x + I y])]], {x, 
                              -2, 2}, {y, ##1}, 
                            PlotPoints -> {32, 8}, 
                            DisplayFunction -> Identity]&, {{-2, 
                            -1 - ε}, {
                            -1 + ε, -ε}, {
                            ε, 1 - ε}, {
                            1 + ε, 2 - ε}}, {
                          1}]], Flatten[
                      Outer[List, {-1, 1}, {-1, 1}, {-1, 1}], 
                      2], {1}]], _Polygon, ∞] /. 
              p_Polygon :> holePolygon[p, 0.8]}], 
        BoxRatios -> {1, 1, 1.199999999999999}, 
        ViewPoint -> {-2, 0, 0}]];

In the last example w (z) = (Sqrt[(z^3)]) + (z^(-3)) (or 1 + (w^3) (z^3) - (z^6) = 0). In this section the determination of the form of the branch cuts becomes a bit more complicated—the discontinuities will no longer be located on a straight line. Here the principal sheet is shown:

Plot3D[
      Evaluate[Im[((((z^3) + (1/(z^3))))^(1/3)) /. z -> x + I y]], {x, 
        -2, 2}, {y, -2, 2}, PlotPoints -> 80];

Using a ContourPlot allows us a nice view of the branch cuts. They form clusters of contour lines.

ContourPlot[
      Evaluate[Im[((((z^3) + (1/(z^3))))^(1/3)) /. z -> x + I y]], {x, 
        -(3/2), (3/2)}, {y, -(3/2), (3/2)}, 
      PlotPoints -> 200, ColorFunction -> (Hue[Random[]]&), 
      Contours -> 30, 
      ContourStyle -> {{Thickness[0.0001], GrayLevel[0]}}];

Again all branch cuts are caused by the branch cut of the Power function, and occur when (z^3) + (z^(-3)) is equal to a negative real number ν.

sol = Solve[z^3 + 1/z^3 == ν], z]

Analyzing these solutions in detail shows that they form rays originating from -1, (((-1))^(1/3)), and (((-1))^(2/3)) projecting radially outwards and intersecting arcs of the unit circle between {π /6, π /2}, {5π/6, π}, and {3π/2, 11π/6}. The branch points are at 0,-i,i,-(((-1))^(1/6)),(((-1))^(1/6)),-(((-1))^(5/6)),"\!\(\*StyleBox[\"\\\"and\\\"\", \"TI\"]\)" (((-1))^(5/6)).

Show[
      (ParametricPlot[
              Evaluate[({Re[#1], Im[#1]}&)/@(z /. sol)], #1, 
              DisplayFunction -> Identity]&)
          /@{{ν, -8, -2}, {ν, -2, 0}}, 
      DisplayFunction -> $DisplayFunction, PlotRange -> All, 
      Frame -> True, Axes -> False, AspectRatio -> Automatic];

Knowing the location of the branch cuts now allows us to split the complex z-plane into regions such that inside each region we can use ParametricPlot3D to get a smooth surface.

\(\(With[{ε = 10\^\(-10\), 
        ϵ = E\^\(2/3\ I\ π\)}, 

      Show[Graphics3D[{EdgeForm[], Thickness[0.0001], 
            SurfaceColor[Hue[0.26], Hue[0.33], 1.88], 
	
            Cases[Table[
                  Apply[ParametricPlot3D[{r\ Cos[φ], 
                          r\ Sin[φ], 
                          Im[\(ϵ\^j\) 
                              \@\(\((r\ E\^\(I  
                                        φ\))\)\^3 + 
                                  \((r\ E\^\(I  
                                        φ\))\)\^\(-3
                                      \)\)\%3]}, #1, #2, 
                        PlotPoints -> {10, 22}, 
                        DisplayFunction -> Identity]&, 
                    
{\ {{r, ε, 1 - ε}, {
                          φ, π/3 + ε, 
                          π - ε}}, 

                      		{{r, ε, 1 - ε}, {
                          φ, \(-π\) + ε, 
                          \(-π\)/3 - ε}}, 

                      	\ \ {{r, ε, 1 - ε}, {
                          φ, \(-π\)/3 + ε, 
                          π/3 - ε}}, 

                      		{{r, 1 + ε, 2}, {φ, 
                          π/3 + ε, 
                          π - ε}}, 

                      		{{r, 1 + ε, 2}, {φ, 
                          \(-π\) + ε, 
                          \(-π\)/3 - ε}}, 

                      		{{r, 1 + ε, 2}, {φ, 
                          \(-π\)/3 + ε, 
                          π/3 - ε}}}, {1}], 

                  						{j, 0, 2}], _Polygon, ∞] /. 
              p_Polygon :> holePolygon[p, \ 0.75]}], 
        PlotRange -> {\(-3.4\), 3.4}, BoxRatios -> {1, 1, 1.9}, 
        Boxed -> False]]; \)\)

Let's stop here with making 3D pictures of simple Riemann surfaces. We could go on and construct Riemann surface with many more sheets, such as the following sum of five square roots \(TraditionalForm
      \`w \((z)\) = \@\(z - 2\) + \@\(z - 1\) + \@z + \@\(z + 1\) + 
        \@\(z + 2\)\). 
The implicit form this time is rather large—it is a polynomial with 81 terms of degree 32 in w.

Short[
      Factor[GroebnerBasis[
            Numerator[
              Together[
                Apply[Subtract, {z - 2 == (aux1^2), z - 1 == (aux2^2), 
                    z == (aux3^2), z + 1 == (aux4^2), z + 2 == (aux5^2), 
                    w == aux1 + aux2 + aux3 + aux4 + aux5}, {1}]]], {w, z}, {
              aux1, aux2, aux3, aux4, aux5}, 
            MonomialOrder -> EliminationOrder][[1
          ]]], 8]

Here the cross section y = 0 is shown.

With[{ε = (1/(10^12))}, 
      Plot[Evaluate[
          Flatten[Table[
              Table[Im[
                  (∑_(k = -2)^2) (Subscript[ϵ, k]) 
                      (Sqrt[k + x + α I ε])], 
                Evaluate[
                  Sequence@@
                    Table[{(Subscript[ϵ, k]), -1, 1, 2}, {k, -2, 2}]]], {
                α, -1, 1, 2}]]], {x, -3, 3}, 
        Frame -> True, Axes -> False, PlotRange -> All, 
        PlotPoints -> 200, PlotStyle -> {Thickness[0.005]}]];

Obstacles occurring in dealing with more complicated polynomials

Let us now try to deal with a slightly more complicated example, say

𝓅 = -(z^4) - 5 (z^3) - 2 (z^2) w - (z^2) (w^2) + (w^3)

p is a polynomial of degree 3 in w, so the Riemann surface has 3 sheets. Actually we can compute the three solutions for p = 0 explicitly:

(Subscript[sol, 1]) = Solve[𝓅 == 0, w]

But trying to directly plot these three solutions reveals the fact the they are not continuous functions of z. Here the real part of each of the three solutions is shown:

Show[
      GraphicsArray[
        (Plot3D[Re[#1], {x, -4, 4}, {y, -4, 4}, 
                DisplayFunction -> Identity, PlotPoints -> 25]&)/@
          (w /. (Subscript[sol, 1]) /. z -> x + I y)]];

The imaginary part shows similar discontinuities:

Show[
      GraphicsArray[
        (Plot3D[Im[#1], {x, -4, 4}, {y, -4, 4}, 
                DisplayFunction -> Identity, PlotPoints -> 25]&)/@
          (w /. (Subscript[sol, 1]) /. z -> x + I y)]];

The lines of discontinuity form a fairly complicated pattern. They arise ultimately from the branch cut of the Power function. We can get a picture showing just these lines of discontinuity by taking all arguments of the Power function with rational second argument and checking when this argument is purely a negative real. The branch points are determined as the z-values where the arguments of the Power function becomes 0.

powerArguments𝓅 = {
        643 + 450 z + 59 (z^2) + 20 (z^3) + 4 (z^4), 
        135 (z^3) + 45 (z^4) + 2 (z^6) + 
          3 (Sqrt[3]) (z^3) (Sqrt[643 + 450 z + 59 (z^2) + 20 (z^3) + 4 (z^4)])};

bps𝓅 = 
      z /. Flatten[N[(Solve[#1 == 0, z]&)/@powerArguments𝓅], 1]

To visualize the branch cuts, we first make a ContourPlot by showing only the lines where the imaginary part vanishes and then delete the points where the real part is positive. The following picture shows the branch cuts as well as the branch points:

Show[
      Graphics[{
          Function[rootArg, 
              Map[If[Re[rootArg /. z -> #1 . {1, I}] > 0, 
                    Sequence@@{}, #1]&, 
                Cases[Graphics[
                    ContourPlot[
                      Evaluate[Im[rootArg /. z -> x + I y]], {x, 
                        -6, 6}, {y, -6, 6}, PlotPoints -> 250, 
                      Contours -> {0}, ContourShading -> False, 
                      DisplayFunction -> Identity]], _Line, 
                  ∞], {-2}]]/@powerArguments𝓅, {
            PointSize[0.03], Hue[0], 
            (Point[{Re[#1], Im[#1]}]&)/@bps𝓅}}], 
      PlotRange -> All, Frame -> True, 
      AspectRatio -> Automatic];

The last picture shows that using the closed form solution for the roots and avoiding traversing the branch cut will be not easy. Trying the same with Root objects gives similar discouraging results:

(Subscript[sol, 2]) = {ToRules[Roots[𝓅 == 0, w, Cubics -> False]]}

Show[
      GraphicsArray[
        (Plot3D[Im[#1], {x, -4, 4}, {y, -4, 4}, 
                DisplayFunction -> Identity, PlotPoints -> 25]&)/@
          (w /. (Subscript[sol, 2]) /. z -> x + I y)]];

The explicit locations of the lines of discontinuity are unfortunately not easy to calculate for Root objects either. (Actually they are more complicated to calculate than the above example involving radicals, because this time the discontinuities are caused by the switch of the root numbers.) So we must conclude that our first naive trial to generate a Riemann surface of a generic algebraic function failed. (It's possible to determine the lines of discontinuity, but it would turn out that calculating the branch cuts and constructing appropriate meshes is not easy.) 

A second possibility which comes to mind is the following strategy: 
p = 0 means that simultaneously the following two equations must hold

ℜ((∑_(i, j = 0)^n)(Subscript[a, ij]) (((x+I y))^j)(((u+I v))^k))=0

ℑ((∑_(i, j = 0)^n)(Subscript[a, ij]) (((x+I y))^j)(((u+I v))^k))=0

These two equations in the four real variables x, y, u, v allow us to eliminate one variable, say u, to obtain one implicit equation in x, y, v. Here this done for the polynomial p from above.

ReIm𝓅 = 
      Expand[𝓅 /. {z -> x + I y, w -> u + I v}]

Im𝓅 = Plus@@Expand[Cases[ReIm𝓅, _Complex _]/I]

Re𝓅 = Expand[ReIm𝓅 - I Im𝓅]

(Subscript[𝓅, 2]) = Resultant[Re𝓅, Im𝓅, u]

Using the function ContourPlot3D from the package Graphics`ContourPlot3D` we can display the surface defined by the implicit equation (Subscript[((𝓅)), 2]) = 0.

Needs["Graphics`ContourPlot3D`"]

Show[Graphics3D[{EdgeForm[], Thickness[0.0001], SurfaceColor[Hue[0.2], Hue[0.45], 1.34],
     Cases[ContourPlot3D[Evaluate[N[Subscript[𝓅, 2]]], {x, -4, 4}, {y, -4, 4},
        {v, -4, 4}, PlotPoints -> 25, MaxRecursion -> 0, DisplayFunction -> Identity,
        PlotRange -> All], _Polygon, Infinity] /. p_Polygon :> holePolygon[p, 0.8]}]];

Although in principle this method works—but with PlotPoints  ->  25 in the last picture one still can see some regions where a finer division would be appropriate—for higher order polynomials the Resultant calculation will become much more expensive. Also the generation of a high resolution ContourPlot3D requires much CPU time, and this method will be unable to deal, for instance, with the following example of a two sheeted hyperelliptic curve with 100 branch points. (In the following contour plot of the imaginary part of the square root of a random polynomial 𝓅 of degree 100, the branch points are shown in blue, the branch cuts are shown in red, the black areas are characterized by "\[GothicCapitalI]"(𝓅) < 0 and the white ones by "\[GothicCapitalR]"(𝓅) < 0.)

SeedRandom[222]Module[{branchPoints = Table[Random[Complex], {100}], hyperellipticF, 
        cp, gr, branchCuts}, 
      hyperellipticF[z_] = Times@@(z - branchPoints); 
      cp = ContourPlot[
          Evaluate[Im[hyperellipticF[x + I y]]], {x, -0.1, 1.1}, {y, 
            -0.1, 1.1}, Contours -> {0}, PlotPoints -> 400, 
          DisplayFunction -> Identity]; 
      branchCuts = 
        Map[If[Re[hyperellipticF[{1, I} . #1]] > 0, Sequence@@{}, #1]&, 
          Cases[Graphics[cp], _Line, ∞], {-2}]; 
      Show[{cp, 
          Graphics[{{Thickness[0.004], Hue[0], branchCuts}, {PointSize[0.01], 
                Hue[0.7], (Point[{Re[#1], Im[#1]}]&)/@branchPoints}}]}, 
        DisplayFunction -> $DisplayFunction]];

Other methods which come to mind to construct the Riemann surface of an algebraic function include the Weierstrass continuation method [Henrici 1974] (for an example see [Trott 1997]), representing the sheets as hypergeometric functions [Appell and Kampé de Feriet 1926; Eisenstein 1844; Trott and Adamchik 1994] or as Siegel Theta functions [Lindemann 1892; Lindemann 1884; Umemura 1984]. Due to their complexity, all of these methods are limited to treating relatively small examples. 
That's why in the following sections we will take another approach.

Outline of the construction

The last section showed that the main obstacle in making a smooth picture of all the sheets is due to branch cuts of radicals and/or of Root objects. To avoid discontinuities we would have to carefully restrict the regions where we plot the function. On the other hand, branch cuts can be chosen to a large extent arbitrarily, and only the location of the branch points represents a property inherent and well defined for an algebraic function. (Branch cuts have to connect branch points, but how exactly a branch cut is chosen in the principle sheet between the branch points doesn't matter). That's why we will avoid using closed form radical or Root formulas and branch cuts at all and instead use NDSolve to calculate smooth patches inside branch-point-free regions of the complex z-plane. This treatment avoids completely dealing with branch cuts (at least on a Riemann surface). The patches will be chosen in such a way that the branch points (which must be treated carefully, where two or more sheets are glued together) are always at their edges . This also allows a good spatial resolution near the branchpoints, which is important because typically w (z) varies most near a branch point. This will result in a smooth, in 3D self-intersecting surfaces representing the Riemann surface of the polynomial under consideration. 
To construct the Riemann surface we will proceed with the following, partially symbolic, partially numerical steps:

1. We use the discriminant to calculate all potential branch points 𝓅.

2. We divide a circular part of the complex z-plane radially and azimuthally into sectors generated by the outer product of the radial and azimuthal coordinates of the branch points, such that every sector has no branch point inside.

3. We derive a differential equation for (w^′)(z) from 𝓅.

4. By starting (very) near the branch points we solve the differential equation radially outwards for all sheets of every sector.

5. Starting from the radial solutions of the differential equations we solve the differential equation azimuthally inside the sectors constructed in (3.). This results in parametrized (by the radius r and the angle φ) parts of all sheets for each sector.

6. We generate polygons from the parametrizations for all sheets and all sectors and display all things together (and maybe cut holes into the polygons for better visibility of the inner parts).

Implementing the construction

1. At a branch point w(z) has a multiple solution. This means 𝓅 and (∂_w) 𝓅 must both be zero. (For a detailed discussion of the local behaviour of algebraic functions near branch points see [Cano 1995; Chudnovsky and Chudnovsky 1987; McDonald 1995; Teissier 1995; Walker 1950].)

discriminant[poly_, w_] := 
      Factor[Resultant[poly, (∂_w) poly, w]]

branchPoints[poly_, {w_, z_}] := 
      Union[z /. 
          NSolve[discriminant[poly, w] == 0, z, WorkingPrecision -> 25]]

As an example for the 𝓅 let us take the following (we will use this polynomial throughout this section for exemplifying the various functions defined below):

𝓅 = -z - 3 (w^2) (z^2) + (w^6) - (z^3) w + 4;

This function has 18 potential branch points:

Short[bPs = branchPoints[𝓅, {w, z}], 6]

Here their location in the complex z-plane is shown.

Show[
      Graphics[{PointSize[0.02], (Point[{Re[#1], Im[#1]}]&)/@bPs}], 
      PlotRange -> All, AspectRatio -> Automatic, 
      Frame -> True, Axes -> False];

2. Given the branch points we calculate their distances from the origin and their phase and divide the r-φ-plane into radial and azimuthal sectors which are formed between radially and azimuthally adjacent branch points.

sectors[poly_, {w_, z_}] := 
      Module[{ε = (1/(10^5)), bps = branchPoints[poly, {w, z}], 
          arg, rMax, rList, φList}, 
        arg[ξ_] := If[ξ == 0. , 0``25, Arg[ξ]]; 
        rMax = (If[#1 == 0. , 1, ((12 #1)/10)]&)[
            Max[Abs[rList = Abs/@bps]]]; 
        rList = Prepend[Append[rList, rMax], 0] ⋃ 
            (SameTest -> (Abs[#1 - #2] < ε&)); 
        φList = 
          arg/@bps ⋃ 
            (SameTest -> (Abs[#1 - #2] < ε&)); 
        φList = 
          Append[φList, First[φList] + 2 π]; 
        (N[#1, 25]&)[
          Flatten[Table[{
                rList[[i]] + 
                  ε, 
                φList[[j]] + 
                  ε, 
                rList[[i + 1]] - 
                  ε, 
                φList
                    [[j + 1]] - 
                  ε}, {i, Length[rList] - 1}, {j, 
                Length[φList] - 1}], 1]]]

Here the sectors for our example polygon 𝓅 are calculated:

Short[s1 = sectors[𝓅, {w, z}], 12]

Let's visualize the sectors together with the branch points by coloring each sector with a randomly chosen color:

sectorPolygon[{r1_, φ1_, r2_, φ2_}] := 
      With[{pp = 12}, 
        Polygon[Join[
            Table[r1 {Cos[φ], Sin[φ]}, {φ, 
                φ1, φ2, 
                ((φ2 - φ1)/pp)}], 
            Table[r {Cos[φ2], Sin[φ2]}, {r, r1, r2, 
                ((r2 - r1)/pp)}], 
            Table[r2 {Cos[φ], Sin[φ]}, {φ, 
                φ2, φ1, 
                ((φ1 - φ2)/pp)}], 
            Table[r {Cos[φ1], Sin[φ1]}, {r, r2, r1, 
                ((r1 - r2)/pp)}]]]]

Show[
      Graphics[{({Hue[Random[]], sectorPolygon[#1]}&)/@s1, 
          PointSize[0.02], (Point[{Re[#1], Im[#1]}]&)/@bPs}], 
      PlotRange -> All, Frame -> True, 
      AspectRatio -> Automatic];

3. Now let's derive a differential equation equivalent to the polynomial under consideration.

The solution of the differential equation we are looking for should reproduce the original polynomial. If we integrate the differential equation and avoid branch points along the integration path, then this solution will be a smooth function representing the original polynomial. Let us take a simple example to demonstrate this: w(z) = (Sqrt[z]). This equation is equivalent to (((w(z)))^3) = z. Differentiating this equation with respect to z gives the differential equation (w^′)(z) (((w(z)))^2) = 1. Solving this differential equation along the unit circle (after changing the independent variable from z to z -> φ with \(TraditionalForm
      \`z = e\^\(i\ φ\)\) gives:

nsol = 
      NDSolve[{ψ'[φ]/(I Exp[I φ]) 
              ψ[φ]^2 == 1, ψ[0] == 1}, 
        ψ, {φ, 0, 8  π}]

Because the function w (z) = (Sqrt[z]) has a Riemann surface with three sheets, we need three revolutions around the origin to come back to the starting value.

\(\(Plot[
      Evaluate[Re[ψ[φ] /. nsol]], {φ, 0, 8\ π}, 
      \ 
	\ \ \ \ \ \ \ \ 
      Ticks -> {Table[j\ π, {j, 8}], Automatic}]; \)\)

Now let us return to the general case. Also here we will derive a first-order nonlinear differential equation. (We could also generate higher-order linear differential equations with polynomial coefficients [see [Cockle 1860; Cockle 1862; Harley 1862; Lachtin 1898] for details], but for higher-order equations getting the initial conditions is no longer so easy—for first-order cases one just solves numerically the polynomial under consideration for a given z with respect to w.) 

If we view w as a function of z, we can differentiate the polynomial with respect to z and solve for (w^′)(z).

ode[poly_, {w_, z_}]:=Equal@@Solve[(∂_z)(poly /. w -> w[z]) == 0,w'[z]][[1, 1]]

The nonlinear first-order differential equation we get this way for our example polynomial is the following:

ode𝓅 = ode[𝓅, {w, z}]

Later on we will need to solve this differential equation in azimuthal and radial directions. 
Transforming the independent variable from z to φ, we get the differential equation in azimuthal direction:

odeφ[ode_, {w_, z_}, {ψ_, φ_}, r_]:=ode//.{w'[z]==rhs_->ψ'[φ]==rhs 
                  (∂_φ)(r Exp[I φ]),w[z] -> ψ[φ],z -> r Exp[I φ]}

ode𝓅φ = 
      odeφ[ode𝓅, {w, z}, {ψ, φ}, r]

Transforming the independent variable from z to r, we get the differential equation in radial direction dir.

oder[ode_, {w_, z_}, {ψ_, ρ_}, dir_]:=ode//.{w'[z]==rhs_->ψ'[ρ]==rhs (∂_ρ)(ρ dir),w[z] -> ψ[ρ],z -> dir ρ}

ode𝓅r = oder[ode𝓅, {w, z}, {ψ, ρ}, dir]

The points (4.), (5.) and (6.) we will carry out not sequentially, but rather for every sector at once. The advantage is that we do not have to store a large number of InterpolatingFunction objects (which can be quite memory intensive). 
For the numerical solution, we set some options in ND- Solve. (The following option settings will be suitable for the examples treated in the following—more complicated polynomials might need larger values of the maximum number of steps as well of the precision/accuracy goals.)

SetOptions[NDSolve, PrecisionGoal -> 12, 
      AccuracyGoal -> 12, MaxSteps -> 2000];

The function patch calculates polygons of all sheets of the polynomial poly using ppφ azimuthal and ppr radial plot points in a given sector. (We give the radial as well as the azimuthal differential equation as arguments to patch. Although they could be calculated from poly, later on we have to call patch on each of the sectors, and we would like to avoid recalculating the differential equation for every sector.)

patch[
        reIm_, {poly_, {w_, z_}}, {odeφ_, {ψ_, φ_}, 
          r_}, {oder_, {ψ_, ρ_}, dir_}, {r1_, φ1_, r2_, 
          φ2_}, {ppφ_, ppr_}] := 
      Module[{radialStartingValues, radialIFs, azimuthalIFs, points, ra, 
          φa}, 
        radialStartingValues = 
          w /. Solve[poly == 0 /. {z -> r1 Exp[I φ1]}, w]; 
        radialIFs = 
          (NDSolve[{oder /. dir -> Exp[I φ1], 
                      ψ[r1] == #1}, ψ, {ρ, r1, r2}]
                    [[1, 1, 2]]&)/@
            radialStartingValues; 
        azimuthalIFs = 
          Table[ra = r1 + ((i (r2 - r1))/ppr); 
            NDSolve[{odeφ /. r -> ra, 
                  ψ[φ1] == 
                    radialIFs[[j]][
                      ra]}, ψ, {φ, φ1, φ2}]
                [[1, 1, 2]], {j, 
              Length[radialIFs]}, {i, 0, ppr}]; 
        points = 
          Table[ra = r1 + ((i (r2 - r1))/ppr); 
            φa = 
              φ1 + 
                ((k (φ2 - φ1))/pp)φ; {
              ra Cos[φa], ra Sin[φa], 
              reIm[azimuthalIFs[[j, 
                    i + 1]][φa]]}, {j, 
              Length[radialIFs]}, {i, 0, ppr}, {k, 0, ppφ}]; 
        Function[s, 
            Apply[Polygon[Join[#1, Reverse[#2]]]&, 
              Transpose/@Partition[(Partition[#1, 2, 1]&)/@s, 2, 1], {
                2}]]/@points]

Here the six sheets inside two of the sectors of our example polynomial are calculated. The right-hand example contains a clearly visible branch point.

Show[
      GraphicsArray[
        (Graphics3D[
                patch[Im, {𝓅, {w, z}}, {
                    ode𝓅φ, {ψ, φ}, r}, {
                    ode𝓅r, {ψ, ρ}, dir}, 
                  s1[[#1]], {8, 8}], 
                BoxRatios -> {1, 1, 2}, ViewPoint -> {2, 1, 0.6}]&)
            /@{6, 111}]];

For the whole Riemann surface we want the mesh as uniform as possible, and at the same time we want the branch points always to be (near) mesh points. The function patchPlotPoints calculates the number of radial and azimuthal points inside a given sector, provided the whole number of radial and azimuthal points should be around ppr and ppφ. The actual number of plot points might slightly deviate from ppr and ppφ.

\(patchPlotPoints[{r1_, φ1_, r2_, φ2_}, 
        rMax_, {ppφ_, ppr_}] := 
      \((Max[#1, 1]&)\)/@
        Round[\({ppφ, ppr}\ {φ2 - φ1, r2 - r1}
                \)\/{2\ π, rMax}]\)

Now we have all pieces together to implement the function RiemannSurface. Its first argument poly is a bivariate polynomial in w and z, the second argument is the function to be shown (in most cases Re or Im) and the last argument is the number of radial and azimuthal plot points. The function RiemannSurface returns a list of polygons of all sheets.

RiemannSurface[poly_, {w_, z_}, reIm_, {ppφ_, ppr_}] := 
      Module[{s, rMax, odePoly, odePolyφ, odePolyr}, 
        s = sectors[poly, {w, z}]; 
        rMax = Max[(#1[[3]]&)/@s]; 
        odePoly = ode[poly, {w, z}]; 
        odePolyφ = 
          odeφ[odePoly, {w, z}, {ψ, φ}, r]; 
        odePolyr = oder[odePoly, {w, z}, {ψ, ρ}, dir]; 
        Table[patch[
            reIm, {poly, {w, z}}, {odePolyφ, {ψ, φ}, 
              r}, {odePolyr, {ψ, ρ}, dir}, 
            s[[i]], 
            patchPlotPoints[s[[i]], 
              rMax, {ppφ, ppr}]], {i, Length[s]}]]

Here one half of the Riemann surface of 𝓅 is shown. A cross-sectional view along the plane y = 0 shows nicely the branching at the six branch points along the real axis.

Show[
      Graphics3D[{EdgeForm[Thickness[0.0001]], 
          SurfaceColor[Hue[0.45], Hue[0.53], 1.78], 
          RiemannSurface[𝓅, {w, z}, Im, {50, 36}]}], 
      PlotRange -> {All, {0, 6}, All}];

By using only one plot point in the azimuthal direction and replacing polygons by lines, we can make a picture showing cross section along all directions which contain branch points. One sees nicely the splitting of the branches near branch points. (Generically one expects most branch points to be of square root type.) Both the real as well as the imaginary part are shown.

Show[
      GraphicsArray[
        Function[l, 
            Function[reIm, 
                Map[If[Head[#1] === List, MapAt[reIm, #1, 3], #1]&, 
                  l, {-2}]]/@{Re, Im}][
          Graphics3D[{Thickness[0.0001], 
              RiemannSurface[𝓅, {w, z}, Identity, {1, 30}] /. 
                Polygon[{p1_, p2_, p3_, p4_}] 
                  :> {Line[{p2, p3}], Line[{p4, p1}]}}, 
            ViewPoint -> {-2, -1, 0.7}]]]];

Some examples

In this section we will use the function RiemannSurface to visualize a collection of Riemann surfaces of algebraic functions. 
We start with the simple function w (z) = (Sqrt[1 - (z^2)]) or (w^2) + (z^2) - 1 = 0, but this time we show the phase instead of the real or imaginary part.

𝓅 = (w^2) + (z^2) - 1;

Show[
      Graphics3D[{EdgeForm[], Thickness[0.0001], GrayLevel[0], 
            SurfaceColor[Hue[0.5], Hue[0.5], 2], 
            RiemannSurface[𝓅, {w, z}, Arg[#1]^2&, {60, 30}]} /. 
          p_Polygon :> holePolygon[p, 0.8]], 
      PlotRange -> All, BoxRatios -> {1, 1, 1}, 
      ViewPoint -> {1.9, -2.2, 0.6}];

Two quartics (symmetric in z and w):

Show[
      GraphicsArray[
        (Graphics3D[{EdgeForm[], Thickness[0.0001], RGBColor[0, 1, 0], 
                    SurfaceColor[Hue[0.78], Hue[0.22], 1.06], 
                    RiemannSurface[#1, {w, z}, Im, {60, 30}]} /. 
                  p_Polygon :> holePolygon[p, 0.7], 
                PlotRange -> {All, {0, 1.399999999999999}, All}, 
                ViewPoint -> {0, -2, 0}, 
                BoxRatios -> {1, 1, (3/2)}]&)/@{
            1 + z + (z^2) + (z^3) + (z^4) + w + (w^2) + (w^3) + (w^4), 
            1 + z + (z^2) + (z^3) + (z^4) - w - (w^2) - (w^3) - (w^4)}]];

Here is an elliptic curve.

𝓅 = w^2 + z^3 + 3*z^2 - 1;

Show[
      Graphics3D[{EdgeForm[], Thickness[0.0001], GrayLevel[0], 
            SurfaceColor[Hue[0.3], Hue[0.8199999999999999], 1.02], 
            RiemannSurface[𝓅, {w, z}, Im, {60, 30}]} /. 
          p_Polygon :> holePolygon[p, 0.8]], 
      PlotRange -> All, BoxRatios -> {1, 1, 1}];

The Riemann surface of a hyperelliptic curve. Coloring the individual patches differently makes it easier to spot the 8 branch points located at the unit circle.

𝓅 = (w^2) - (z^8) - 1;

To change the appearance of these pictures, we implement another form for making holes in the surface.

starPolygon[Polygon[{p1_, p2_, p3_, p4_}], f_] := 
      Module[{m = (1/4) (p1 + p2 + p3 + p4), q1 = ((p4 + p1)/2), 
          q2 = ((p1 + p2)/2), q3 = ((p2 + p3)/2), q4 = ((p3 + p4)/2), r1, 
          r2, r3, r4}, {r1, r2, r3, r4} = 
          (m + f (#1 - m)&)/@{q1, q2, q3, q4}; {
          Polygon/@{{r1, p1, r2}, {r2, p2, r3}, {r3, p3, r4}, {r4, p4, r1}}, 
          Line/@{{r1, r2, r3, r4, r1}, {p1, r2, p2, r3, p3, r4, p4, r1, p1}}}]

Show[Graphics3D[{EdgeForm[], Thickness[0.0001], GrayLevel[0],
      Map[{SurfaceColor[Hue[Random[]], Hue[Random[]], 3*Random[]], #1} & ,
       RiemannSurface[𝓅, {w, z}, Im, {60, 30}], {2}]} /.
     p_Polygon :> starPolygon[p, 0.5]], PlotRange -> All, BoxRatios -> {1, 1, 1}];

The Riemann surface of a generalized hyperelliptic function:

𝓅 = -1 + z - 4*z^2 + z^3 - z^4 + z^5 - z^6 + w^4;

Show[
      Graphics3D[{EdgeForm[], Thickness[0.0001], GrayLevel[1], 
            SurfaceColor[Hue[0.04], Hue[0.2399999999999999], 
              1.889999999999999], 
            RiemannSurface[𝓅, {w, z}, Im, {60, 30}]} /. 
          p_Polygon :> starPolygon[p, 0.3]], 
      PlotRange -> All, Background -> GrayLevel[0]];

The Riemann surface of the cubic (in w) from the section above. The polygons are colored according to their distance from the z-axis. The function holePolygon1 also cuts a hole into a polygon, but in a slightly different manner than holePolygon and starPolygon.

holePolygon1[Polygon[{p1_, p2_, p3_, p4_}]] := 
      Function[{q1, q2, q3, q4}, {
            Polygon/@{{q1, p1, q2}, {q2, p2, q3}, {q3, p3, q4}, {q4, p4, 
                  q1}}, Line[{q1, q2, q3, q4, q1}]}][((p4 + p1)/2), 
        ((p1 + p2)/2), ((p2 + p3)/2), ((p3 + p4)/2)]

𝓅 = -(z^4) - 5 (z^3) - 2 (z^2) w - (z^2) (w^2) + (w^3);

Show[
      Graphics3D[{EdgeForm[], GrayLevel[0.6], 
          ({(SurfaceColor[Hue[#(1/45)], Hue[#(1/45)], 2.299999999999999]&)
                      [(#1 . #1&)[
                      Drop[(1/4) 
                          Plus@@#1[[1]], 
                        -1]]], #1}&)/@
            DeleteCases[
              Flatten[RiemannSurface[𝓅, {w, z}, Im, {70, 30}] /. 
                  p_Polygon :> holePolygon1[p]], _Line]}], 
      PlotRange -> All, BoxRatios -> {1, 1, 1}];

Another cubic (in w):

𝓅 = 7 + 10 (w^2) - 10 (w^3) - 2 (z^3) - 15 (z^4);

Show[
      Graphics3D[{EdgeForm[Thickness[0.0001]], 
            SurfaceColor[Hue[0.32], Hue[0.4699999999999999], 2.24], 
            RiemannSurface[𝓅, {w, z}, Im, {60, 25}]} /. 
          p_Polygon :> starPolygon[p, 0.7]], 
      BoxRatios -> {1, 1, 1}, ViewPoint -> {-2, 2, 0.3}];

The Riemann surface of a trinomial sextic.

𝓅 = (w^6) - w - z;

Show[
      Graphics3D[{EdgeForm[], Thickness[0.0001], GrayLevel[0], 
            SurfaceColor[Hue[0.22], Hue[0.78], 1.52], 
            RiemannSurface[𝓅, {w, z}, Im, {60, 30}]} /. 
          p_Polygon :> starPolygon[p, 0.8]], 
      PlotRange -> All, BoxRatios -> {1, 1, 2}, 
      ViewPoint -> {2, -2, 0.7}];

The next example has a pole at z = 0:

𝓅 = ((w^3) + w - 1) (z^2) - 1;

Show[
      Graphics3D[{EdgeForm[], Thickness[0.0001], Hue[0], 
            SurfaceColor[Hue[0.52], Hue[0.65], 1.82], 
            RiemannSurface[𝓅, {w, z}, Im, {60, 30}]} /. 
          p_Polygon :> holePolygon1[p]], 
      BoxRatios -> {1, 1, 1}, Axes -> True, 
      PlotRange -> {-5, 5}];

Here is an example with five poles:

𝓅 = 
      ((w^4) + (w^2) - 1) ((z^3) - 1) (((3 (z^2) - 1))^2) - 1;

Show[
      Graphics3D[{EdgeForm[], Thickness[0.0001], Hue[0.22], 
            SurfaceColor[Hue[0.03], Hue[0.56], 1.07], 
            RiemannSurface[𝓅, {w, z}, Im, {80, 40}]} /. 
          p_Polygon :> starPolygon[p, 0.4]], 
      BoxRatios -> {1, 1, (3/2)}, PlotRange -> {-3, 3}];

A randomly chosen example with three poles.

𝓅 = 
      -11 (w^4) + 8 (w^2) (z^2) + 10 w (z^3) + 3 (w^4) (z^3) - 
        10 (z^4) + 5 w (z^5)

Show[
      Graphics3D[{EdgeForm[Thickness[0.0001]], 
            SurfaceColor[Hue[0.067], Hue[0.64], 1.04], 
            RiemannSurface[𝓅, {w, z}, Im, {60, 25}]} /. 
          p_Polygon :> starPolygon[p, 0.55]], 
      BoxRatios -> {1, 1, 1}, Axes -> True];

Compactified Riemann surfaces

Many more algebraic functions can be visualized using the function RiemannSurface, but we will end here. One might now go on and construct pictures where we do not show "\[GothicCapitalI]"(w) or "\[GothicCapitalR]"(w) over the complex z-plane, but rather over the Riemann sphere—this would allow to cover all z-values and result in pictures which show explicitly compact Riemann surfaces ('spheres with ears'). Given the Riemann sphere of radius R = 1/2 around {0, 0, 1/2} we could visualize "\[GothicCapitalI]"(w(x + i y)) as a point in direction of the image of x + i y on the Riemann sphere and distance radius \(TraditionalForm
      \`R + r\ \(arctan("ℑ"(w(x + i\ y)))\)\) (the arctan in the last formula allows to uniquely map the interval {-∞, ∞} to a finite interval). The function RiemannSpherePoint calculates the position of the point corresponding to the spherical coordinates {φ, θ} under the mapping of the function f on the 'stretched' Riemann sphere.

RiemannSpherePoint[{φ_, θ_}, f_] := 
      Module[{ξ, η, d, rMax = 0.2}, {ξ, η} = 
          N[Cot[(θ/2)] {Cos[φ], Sin[φ]}]; 
        d = {Cos[φ] Sin[θ], 
            Sin[φ] Sin[θ], Cos[θ]}; {0, 0, (1/2)} + 
          d ((1/2) + rMax ArcTan[f[ξ, η]])]

The function RiemannSpherePatch generates a patch for the function f in the specified part of the domain of the spherical coordinates.

RiemannSpherePatch[
        f_, {φ1_, φ2_}, {θ1_, θ2_}, {
          ppφ_, ppθ_}] := 
      Module[{ε = (1/(10^12))}, 
        points = 
          Table[RiemannSpherePoint[{φ, θ}, f], {φ, 
              φ1 + ε, φ2 - ε, 
              ((φ2 - φ1 - 
                  2 ε)/pp)φ}, {θ, 
              θ1 + ε, θ2 - ε, 
              ((θ2 - θ1 - 2 ε)/pp)θ}]; 
        Function[s, 
            Apply[Polygon[Join[#1, Reverse[#2]]]&, 
              Transpose/@Partition[(Partition[#1, 2, 1]&)/@s, 2, 1], {2}]]
            [points]]

Here the function (Sqrt[1 - (z^2)]) shown in using the above defined functions. The two branch points at ±1 are mapped into the points {±(1/2), 0, (1/2)}. The north pole corresponds to complex infinity.

rsf = 
      With[{ϵ = -1}, 
        Show[Graphics3D[{EdgeForm[], Thickness[0.0001], GrayLevel[0], 
                SurfaceColor[Hue[0.22], Hue[0.78], 1.919999999999999], 
                Table[Map[
                    RiemannSpherePatch[
                        Function[{x, y}, 
                          Im[(ϵ^i) 
                              (Sqrt[(((x + I y))^2) - 1])]], #1, {0, π}, {
                          20, 30}]&, {{0, (π/2)}, {(π/2), π}, {
                        π, ((3 π)/2)}, {((3 π)/2), 
                        2 π}}, {1}], {i, 0, 1}]} /. 
              p_Polygon :> holePolygon[p, 0.8]], 
          Boxed -> False]];

Despite the holes in the polygons the inner parts are hard to see in the last picture, so lets cut the last surface along the x-y-plane, along the x-z-plane and along the y-z-plane.

Show[
      GraphicsArray[{
          Show[rsf, DisplayFunction -> Identity, 
            PlotRange -> {All, All, {0.5, -0.16}}, 
            ViewPoint -> {0.2, 0, 2.1}], 
          Show[rsf, DisplayFunction -> Identity, 
            PlotRange -> {{0, -0.6}, All, All}, 
            ViewPoint -> {2.1, 0.2, 0}], 
          Show[rsf, DisplayFunction -> Identity, 
            PlotRange -> {All, {0, -0.7}, All}, 
            ViewPoint -> {0, 2.1, 0}]}]];

As another example let us take the function (Sqrt[(z^2)]) - 1, shown on the Riemann sphere. The two branch points are again at ±1 and are mapped into the points {±(1/2), 0, (1/2)}.

rsf = 
      With[{ϵ = (E^(((2 I π)/3)))}, 
        Show[Graphics3D[
            {EdgeForm[], Thickness[0.0001], GrayLevel[0], 
                  Table[Map[
                      RiemannSpherePatch[
                          Function[{x, y}, 
                            Im[(ϵ^i) 
                                ((((((x + I y))^2) - 1))^(1/3))]], #1, {
                            0, π}, {20, 30}]&, {{0, (π/2)}, {(π/2), 
                          π}, {π, ((3 π)/2)}, {((3 π)/2), 
                          2 π}}, {1}], {i, 0, 2}]} /. 
                p_Polygon :> starPolygon[p, 0.4] /. 
              p_Polygon 
                :> {
                  (SurfaceColor[Hue[0.6 #1], Hue[0.45 #1], 1.37]&)[
                    ((Sqrt[#1 . #1])&)[
                      (1/4) Plus@@
                            p[[1]] - {0, 
                          0, (1/2)}]], p}], Boxed -> False]];

Again we cut the last surface along the x-y-plane, along the x-z-plane and along the y-z-plane. In distinction to the above pictures over the complex z-plane, we now also see the branch point at ∞.

Show[
      GraphicsArray[{
          Show[rsf, DisplayFunction -> Identity, 
            PlotRange -> {All, All, {0.5, -0.15}}, 
            ViewPoint -> {0.2, 0, 2.5}], 
          Show[rsf, DisplayFunction -> Identity, 
            PlotRange -> {{0, -0.62}, All, All}, 
            ViewPoint -> {2.5, 0.2, 0}], 
          Show[rsf, DisplayFunction -> Identity, 
            PlotRange -> {All, {0, -0.67}, All}, 
            ViewPoint -> {0, 2.5, 0.2}]}]];

Here is the function (Sqrt[(z^5)]) - 1 shown on the Riemann sphere. The five branch points lie on the equator. We view from above into the lower half of the 'sphere'.

With[{ϵ = (E^(((2 I π)/5)))}, 
      Show[Graphics3D[{EdgeForm[], SurfaceColor[Hue[0.07], Hue[0.27], 1.55], 
            Table[Map[
                RiemannSpherePatch[
                    Function[{x, y}, 
                      Im[(ϵ^i) 
                          ((((((x + I y))^5) - 1))^(1/5))]], #1, {
                      (π/2), π}, {32, 30}]&, 
                Partition[(1/5) π Range[1, 11, 2], 2, 1], {1}], {i, 0, 
                4}]}], Boxed -> False, 
        ViewPoint -> {0.1, 0.9, 3.6}]];

One could now go on and visualize more complicated Riemann surfaces on the Riemann sphere. To do this in the spirit of the last paragraph we would rewrite oder to a differential equation odeθ for the longitudinal coordinate of a spherical coordinate θ system according to

odeθ[oder_, {ψ_, ρ_}, {ϕ_, θ_}]:=oder//.{ψ'[ρ]==rhs_->ϕ'[θ]==(Cos[θ] - 1) rhs,ψ[ρ] -> ϕ[θ],ρ -> Cot[(θ/2)]}

together with corresponding modifications in patch and RiemannSurface. Such an extension is straightforward; we end now and leave further explorations in this direction to the reader. 
So far for the construction of Riemann surfaces of algebraic functions. The construction of Riemann surfaces of transcendental functions is more complicated, a few examples (inverse trigonometric functions, hypergeometric functions, the ProductLog function) will be discussed in parts II to IV (forthcoming) of this series on the visualization of Riemann surfaces.

References

Appell, P., J. Kampé De Feriet. Fonctions hyperéométriques et hypersphériques, Gauthier-Villars, Paris 1926. 

Brill, A.  Vorlesungen über ebene algebraische Kurven und algebraische Funktionen, Vieweg, Braunschweig, 1925. 

Cano, J. , in D. L. Te, K. Saito, Teissier, B. (eds.). Singularity Theory, World Scientific, Singapore, 1995. 

Chudnovsky, D. V., G. V. Chudnovsky. J. Complexity 3, 1 (1987). 

Cockle, J. Phil. Mag. XX, 145 (1860). 

Cockle, J. Phil. Mag. XXIII, 135 (1862). 

Dyck, W. (ed.). Katalog mathematischer und mathematisch physikalischer Modelle, Apparate und Instrumente, Georg Olms Verlag, Heidelberg, 1994. 

Eisenstein, G. J. reine angew. Math. 27, 81 (1844). 

Fischer, G. (ed.). Mathematical Models, Vieweg & Sohn, Braunschweig, 1986. 

Jung, H. W. Einführung in die Theorie der algebraischen Funktionen einer Veränderlichen, de Gruyter, Berlin, 1923. 

Harley, R. Quart. J. pure appl. Math. V, 337 (1862). 

Henrici, P. Applied and Computational Complex Analysis, v. 1, John Wiley & Sons, New York, 1974. 

Lachtin, L. K. Mat. Sbornik XXI, 393 (1898). 

Lindemann, F. Nachrichten der Königl. Gesell. der Wiss. zu Göttingen, 292 (1892). 

Lindemann, F. Nachrichten der Königl. Gesell. der Wiss. zu Göttingen, 245 (1884). 

Mathews, J. H., Howell, R. W., Complex Analysis for Mathematics and Engineering, Jones and Bartlett, Boston 1997. 

McDonald, J. J. Pure Appl. Math. 104, 213 (1995). 

Monna, A. F., Dirichlet's Principle, Ososthoek, Scheltema & Holkema, Utrecht 1975. 

Osgood, W. F., Lehrbuch der Funktionentheorie, Teubner, Leipzig, 1923. 

Rauch, H. E., Lebowitz, A., Elliptic Functions, Theta Functions and Riemann Surfaces, Williams & Wilkins, Baltimore 1973. 

Shihara, K., Sasaki, T. Japan. J. Indust. Appl. Math. 13, 107 (1996). 

Teissier, B., in D. L. Te, K. Saito, Tteissier B. (eds.). Singularity Theory, World Scientific, Singapore, 1995. 

Trott, M., Adamchik, v., Solving the Quintic with Mathematica, http://www.mathsource.com/MathSource/Applications/Mathematics/ 0207 

Trott, M. The Mathematica Guidebook, TELOS, 1998, forthcoming. 
umemura, h. in d. mumford. Tata Lectures on Theta II, Birkhäuser, Boston 1984. 

Walker, R. J. Algebraic Curves, Princeton University Press, Princeton, 1950.

About the Author

Michael Trott studied physics at the Humboldt-University of Berlin from 1981 to 1986. He got his Ph.D. in theoretical solid state physics in 1990 at the Technical University of Ilmenau. In 1994 he joined the Research and Development team at Wolfram Research, Inc. His current scientific interests are visualization in mathematics, elliptic functions, the application of computer algebra to solved and unsolved problems from physics and the foundational problems of quantum mechanics. 

Michael Trott
Wolfram Research, Inc.
100 Trade Center Drive
Champaign, IL 61820
mtrott@wolfram.com