How to plot spiral that goes around circular paraboloid

I have a 3D circular paraboloid surface and I would like to plot a spiral that starts from an arbitrary point on the surface and goes down while "hugging" the surface.

This is my attempt so far:

import numpy as np
import matplotlib.pyplot as plt

fig = plt.figure()
ax = plt.axes(projection='3d')

# Surface ------------------
# Create the mesh in polar coordinates and compute corresponding Z
r0 = 5
r = np.linspace(0, r0, 50)
p = np.linspace(0, 2*np.pi, 50)
R, P = np.meshgrid(r, p)
Z = -R**2 + r0**2

# Express the mesh in the cartesian system
X, Y = R*np.cos(P), R*np.sin(P)

# Plot the surface
ax.plot_surface(X, Y, Z, linewidth=0, antialiased=False, alpha=0.2)

# Spiral -------------------
u = np.arange(0, 29, 0.1)
x = 0.17*u*np.cos(u)
y = 0.17*u*np.sin(u)
z = -0.15*u/np.pi*(x**2 + y**2) + r0**2

# Plot spiral
ax.plot3D(x, y, z, 'r')

plt.show()

However, my spiral is not actually following the surface.

I also tried this:

x = []
y = []
z = []
for i in range(50):
    x.append(X[i,i])
    y.append(Y[i,i])
    z.append(-(X[i,i]**2 + Y[i,i]**2) + r0**2)
ax.plot3D(x, y, z, 'b')

which is going around the surface but I don't know how to make it do more circles around the surface. Any ideas?

The formula in the second attempt is correct. I get what you want if I use the same formula in your first attempt.

The line z = -0.15*u/np.pi*(x**2 + y**2) + r0**2 needs to be replaced with -(x**2 + y**2) + r0**2.

For reproducibility:

%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt

fig = plt.figure()
ax = plt.axes(projection='3d')

# Surface ------------------
# Create the mesh in polar coordinates and compute corresponding Z
r0 = 5
r = np.linspace(0, r0, 50)
p = np.linspace(0, 2*np.pi, 50)
R, P = np.meshgrid(r, p)
Z = -R**2 + r0**2

# Express the mesh in the cartesian system
X, Y = R*np.cos(P), R*np.sin(P)

# Plot the surface
ax.plot_surface(X, Y, Z, linewidth=0, antialiased=False, alpha=0.2)

# Spiral -------------------
# Attempt 1
u = np.arange(0, 29, 0.1)
x = 0.17*u*np.cos(u)
y = 0.17*u*np.sin(u)
z = -(x**2 + y**2) + r0**2
# z = -0.15*u/np.pi*(x**2 + y**2) + r0**2
# Plot spiral
ax.plot3D(x, y, z, 'r')


plt.show()

Output is as shown below: