What's the default NormalsFunction for ParametricPlot3D?

For Plot3D, I can use the method here to compute vertex normals,

With[{f = Function[{x, y}, x^2 + y^2]},
Plot3D[f[x, y], {x, -2, 2}, {y, -2, 2}, Mesh -> None,
NormalsFunction ->
Function[{x, y}, Evaluate[Cross @@ Transpose[D[{x, y, f[x, y]}, {{x, y}}]]]]]
]


I wonder how to compute vertex normals for ParametricPlot3D, I've tried

With[{para = {Sin[u] Cos[v] , Sin[u] Sin[v], Cos[u]}},
ParametricPlot3D[para, {u, 0, Pi}, {v, 0, 2 Pi}, Mesh -> None,
NormalsFunction ->
Function[{u, v}, Evaluate[Cross @@ Transpose[D[para, {{u, v}}]]]]]
]


It looks bad. Do you know how the default NormalsFunction is defined, or how can I define a suitable function?

Your math is correct, but a problem with the code arises because With[] causes the literally-appearing parameters of Function[] to have a \$ added to them, but not those in para. The result is that the normal is the zero vector at all points. Also check the arguments to Function for a NormalsFunction for ParametricPlot3D[].
With[{para = {Sin[u] Cos[v], Sin[u] Sin[v], Cos[u]}},