[Edit: the following answer addresses the case of the line through the center of the sphere.]
Note: there is an ambiguity due to the fact that "angle $\theta$" does not specify which way the rotation is going (you can call the North Pole the South Pole, and suddenly the Earth rotate in the other way!)
What I will assume is that you are given the center $(a,b,c)$ of the sphere and a vector $(\alpha,\beta,\gamma) \neq (0,0,0)$ so that $x_2 = a + \alpha$, $y_2 = b + \beta$, and $z_2 = c + \gamma$. We shall assume the rotation is right handed relative to the vector $(\alpha,\beta,\gamma)$, that is, if you look along the direction given by $(\alpha,\beta,\gamma)$, the rotation is clockwise. Note that for any $\lambda \neq 0$, $(\alpha,\beta,\gamma)$ and $(\lambda \alpha,\lambda\beta,\lambda\gamma)$ determine the same line. But if $\lambda < 0$ their rotational directions are opposite.
For convenience we will require that the vector $(\alpha,\beta,\gamma)$ is a unit vector, that is $\alpha^2 + \beta^2 + \gamma^2 = 1$. We can always get this by dividing the vector by its length.
Then we can use this formula here plus a translation: a point $(x,y,z)$ is sent to
$$ \begin{pmatrix} x \\ y \\ z\end{pmatrix} \mapsto
\begin{pmatrix} a \\ b \\ c\end{pmatrix} +
\begin{pmatrix}
\cos \theta + (1 - \cos \theta)\alpha^2 & \alpha\beta(1-\cos\theta) - \gamma \sin\theta & \alpha\gamma(1-\cos\theta) + \beta \sin\theta \\
\alpha \beta(1-\cos\theta) + \gamma \sin\theta & \cos\theta + (1-\cos\theta)\beta^2 & \beta\gamma(1-\cos\theta) - \alpha\sin\theta\\
\alpha\gamma(1-\cos\theta) - \beta \sin\theta & \beta\gamma(1-\cos\theta) + \alpha\sin\theta & \cos\theta + (1-\cos\theta) \gamma^2
\end{pmatrix}
\begin{pmatrix} x - a \\ y - b \\ z - c\end{pmatrix}
$$
As a first step, I'd move $P_1$ to the origin, so that the points become $P_1=(0,0,0)$, $P_2=(x_2-x_1, y_2-y_1, z_2-z_1)$, $P_3=(x_3-x_1, y_3-y_1, z_3-z_1)$. From there it's just a question of applying a rotation matrix.
A rotation matrix $\mathbf{R}$ which preserves all distances and shapes etc is
given by
$$
\mathbf{R}=\left(
\begin{array}
[c]{ccc}%
e_{1}^{x} & e_{2}^{x} & e_{3}^{x}\\
e_{1}^{y} & e_{2}^{y} & e_{3}^{y}\\
e_{1}^{z} & e_{2}^{z} & e_{3}^{z}%
\end{array}
\right)
$$
where the vectors $\left( e_{1}^{\alpha},e_{2}^{\alpha},e_{3}^{\alpha
}\right) $ for $\alpha\in\left\{ x,y,z\right\} $ are orthogonal and of unit
length, i.e.
$$
\sum_{i=1}^{3}e_{i}^{\alpha}e_{i}^{\beta}=\left\{
\begin{array}
[c]{ccc}%
1 & & \alpha=\beta\\
0 & & \alpha\neq\beta
\end{array}
\right. .%
$$
Applying $\mathbf{R}$, your rotated points $P_{i}^{\prime}=\left(
x_{i}^{\prime},y_{i}^{\prime},z_{i}^{\prime}\right) $ are given by
$$
\left(
\begin{array}
[c]{c}%
x_{i}^{\prime}\\
y_{i}^{\prime}\\
z_{i}^{\prime}%
\end{array}
\right) =\mathbf{R}\left(
\begin{array}
[c]{c}%
x_{i}\\
y_{i}\\
z_{i}%
\end{array}
\right)
$$
I suggest making the top row of $\mathbf{R}$ the unit
vector in the direction from $P_{1}$ to $P_{2}$, i.e.
$$
\left(
\begin{array}
[c]{c}%
e_{1}^{x}\\
e_{2}^{x}\\
e_{3}^{x}%
\end{array}
\right) =\frac{1}{\sqrt{\left( x_{1}-x_{2}\right) ^{2}+\left( y_{1}%
-y_{2}\right) ^{2}+\left( z_{1}-z_{2}\right) ^{2}}}\left(
\begin{array}
[c]{c}%
x_{1}-x_{2}\\
y_{1}-y_{2}\\
z_{1}-z_{2}%
\end{array}
\right)
$$
which would mean that the line from $P_{1}$ to $P_{2}$ gets mapped into the
$x$-axis.
Working out what you could use for the second and third rows of $\mathbf{R}$ is now just a question of solving some simple linear equations, to make sure that the line from $P_{1}$ to $P_{3}$ has no z component.
To solve for $\left( e_{1}^{z},e_{2}^{z},e_{3}^{z}\right)$, the vector
$e_{i}^{z}$ must have a zero dot-product with both $e_{i}^{x}$ and $\left(
x_{3}-x_{1},y_{3}-y_{1},z_{3}-z_{1}\right) $, so the equations are
$$
\begin{align*}
e_{1}^{x}e_{1}^{z}+e_{2}^{x}e_{2}^{z}+e_{3}^{x}e_{3}^{z} & =0\\
\left( x_{3}-x_{1}\right) e_{1}^{z}+\left( y_{3}-y_{1}\right) e_{2}%
^{z}+\left( z_{3}-z_{1}\right) e_{3}^{z} & =0
\end{align*}
$$
Hence
$$
\left(
\begin{array}
[c]{c}%
e_{1}^{z}\\
e_{2}^{z}\\
e_{3}^{z}%
\end{array}
\right) =\lambda_{z}\left(
\begin{array}
[c]{c}%
\left( z_{3}-z_{1}\right) e_{2}^{x}-\left( y_{3}-y_{1}\right) e_{3}^{x}\\
\left( x_{3}-x_{1}\right) e_{3}^{x}-\left( z_{3}-z_{1}\right) e_{1}^{x}\\
\left( y_{3}-y_{1}\right) e_{1}^{x}-\left( x_{3}-x_{1}\right) e_{2}^{x}%
\end{array}
\right)
$$
for some $\lambda_z$, which should be determined so that $e_{i}^{z}$ is a vector
of unit length. Finally $e_{i}^{y}$ can be determined as the vector which is
perpendicular to both $e_{i}^{x}$ and $e_{i}^{z}$, i.e.
$$
\left(
\begin{array}
[c]{c}%
e_{1}^{y}\\
e_{2}^{y}\\
e_{3}^{y}%
\end{array}
\right) =\lambda_{y}\left(
\begin{array}
[c]{c}%
e_{3}^{z}e_{2}^{x}-e_{2}^{z}e_{3}^{x}\\
e_{1}^{z}e_{3}^{x}-e_{3}^{z}e_{1}^{x}\\
e_{2}^{z}e_{1}^{x}-e_{1}^{z}e_{2}^{x}%
\end{array}
\right)
$$
where again $\lambda_{y}$ is determined so that $e_{i}^{y}$ is a vector of
unit length.
Best Answer
Let the vectors be denoted by $v_1$ and $v_2$ (my preference :)).
Get the cross product of the two: $c = v_1 \times v_2$
$c$ is orthogonal to both $v_1$ and $v_2$.
Get the angle between $c$ and $k=(0,0,1)$: $$ \cos\theta = \frac{c\cdot k}{|c|} $$
Rotate both vectors about $a = c \times k$ by $\theta$ using Rodrigues' rotation formula
You have to normalize $a$ to get a unit vector before applying it to the formula: $a' = a/|a|$.
The pair of resulting vectors is one instance of possible rotations and both are perpendicular to $k$ ($z$-axis).
The rotation matrix can be obtained from Rotation matrix from axis and angle