I figured it out, since the sphere is centered at the origin, the cross product between $\vec{N}^{\,}$ and $\vec{Z}^{\,}$ will result in $\vec{R}^{\,}$.
The tangent lines to the sphere form a cone with apex at the given point $\mathbf p$ and axis the line through $\mathbf p$ and the center $\mathbf c$ of the sphere. (To keep things clearer, I’ll mark inhomogeneous Cartesian coordinate vectors with a tilde; if there’s no tilde, this indicates a homogeneous coordinate vector.) The points of tangency lie on the intersection of the sphere with the polar plane of $\mathbf p$. This will be a circle in this polar plane, which is perpendicular to the cone’s axis, with center $\mathbf c'$ at the intersection of the plane and axis. To parameterize this circle, find the radius $r'$ via the Pythagorean theorem, and find a pair of unit vectors $\mathbf u$ and $\mathbf v$ that are orthogonal to the axis and each other. The circle is then $$\mathbf q = \mathbf c'+r'\mathbf u\cos\theta+r'\mathbf v\sin\theta.$$ The line through this point be found in any of the usual ways, such as the parametrization $(1-\lambda)\mathbf p+\lambda\mathbf q$.
To keep things clearer in the following, I’ll mark variables that represent inhomogeneous Cartesian coordinate vectors with a tilde; no tilde indicates a homogeneous coordinate vector. Representing the sphere by the matrix $$S = \left[\begin{array}{c|c} I_3 & -\tilde{\mathbf c} \\ \hline -\tilde{\mathbf c}^T & \|\tilde{\mathbf c}\|^2-r^2\end{array}\right],$$ the polar plane to $\mathbf p$ is $\mathbf\pi=S\mathbf p$. (The components of this vectors correspond to the coefficients of the implicit Cartesian equation of the plane.) Use your favorite method to find its intersection with the line through $\mathbf p$ and $\mathbf c$. Here, I’ll use the Plücker matrix of the line: $$\mathbf c' = (\mathbf p\mathbf c^T-\mathbf c\mathbf p^T)\mathbf\pi = (\mathbf c^T\mathbf\pi)\mathbf p - (\mathbf p^T\mathbf\pi)\mathbf c.$$ (The parenthesized quantities on the r.h.s. are just the dot products of $\mathbf p$ and $\mathbf c$ with $\mathbf\pi$.) The radius of the circle is found from $$r'^2 = r^2-\|\tilde{\mathbf c}-\tilde{\mathbf c}'\|^2.$$ For the unit vectors $\mathbf u$ and $\mathbf v$ of the parametrization, choose one of the cross products of $\tilde{\mathbf c}-\tilde{\mathbf p}$ with the unit basis vectors (at least two of them are guaranteed to be nonzero) and normalize it for $\mathbf u$, and then take $\mathbf u\times(\tilde{\mathbf c}-\tilde{\mathbf p})$, normalized, for $\mathbf v$. You can instead use $\tilde{\mathbf c}-\tilde{\mathbf c}'$ or $\tilde{\mathbf c}'-\tilde{\mathbf p}$ for this calculation; I’m not sure which will provide the best numerical stability.
For your example, I get $$S = \begin{bmatrix} 1&0&0&-100\\0&1&0&50\\0&0&1&-75\\-100&50&-75&17500\end{bmatrix}$$ from which $\mathbf\pi=[400:150:-75:-27500]$. Using the above formula for the axis-plane intersection gives $$\mathbf c' = [400:150:-75:-27500],$$ which corresponds to the inhomogeneous Cartesian coordinates $$\tilde{\mathbf c}' = \left({30500\over301},{14900\over301},{22500\over301}\right)$$ and the circle’s radius works out to be $r'=250\sqrt{\frac3{301}}$.
For the unit vector $\tilde{\mathbf u}$, I'll take the cross product that has the greatest norm, $(-150,400,0)^T$, so $\tilde{\mathbf u} = \frac1{\sqrt{73}}(-3,8,0)$ and (according to Mathematica; the radicals are starting to get a bit ugly) $\tilde{\mathbf v}=\left({24\over\sqrt{21973}},{9\over\sqrt{21973}},2\sqrt{{73\over301}}\right)^T$. The parametrization of the circle finally comes out to be approximately $$x = 101.329 - 8.76349 \cos\theta + 4.04095 \sin\theta \\ y= -49.5017 +
23.3693 \cos\theta + 1.51536 \sin\theta \\ z = 74.7508 +
24.5825 \sin\theta.$$
Depending on what it is you want to do with this, it might be more convenient to work with the tangent cone as a whole. It turns out that its matrix is easily computed from $\mathbf p$ and $C$: $$C_{cone} = (\mathbf p^TC\mathbf p)C - (C\mathbf p)(C\mathbf p)^T.$$ Note that the second term isn’t a dot product; it’s the tensor product of $C\mathbf p$ with itself.
Best Answer
Approach: Intersection of tangent plane to the sphere at P with Y axis, call it Q. Line that contains PQ.
The equation of the tangent plane is:
$$x+9y+\sqrt{18} z+D=0$$
And $$1+81+18+D=0$$
The plane is then $$x+9y+\sqrt{18} z-100=0$$
And Q is $(0, 100/9,0)$
PQ is then $( -1,-9+100/9,-\sqrt{18})$
And the desired line is
$$P + \lambda PQ$$
That is to say:
\begin{align} x&= 1-\lambda \nonumber \\ y&= 9 + \frac{19}{9} \lambda. \nonumber \\ z&= \sqrt{18} -\sqrt{18} \lambda. \nonumber \end{align}
which is the equation of the line in parametric coordinates.
In cartesian coordinates:
$$1-x=\frac{9}{19}(y-9) = \frac{\sqrt{18}-z}{\sqrt{18}}$$
Pictorially: