I'm not entirely sure if you're asking for a closed expression that gives the points, or just a description of these points and am going to assume the latter.
The first step in solving this problem is to change it into an easier one. We may rotate the sphere in such a way the point $P$ becomes the north pole $N=(0,0,1)$. If we call this rotation (which may be expressed as a $3\times3$ matrix) $R$, a point $Q$ is on the circle with radius $r$ around $P$ iff $R(Q)$ is on the circle with radius $r$ around $N$. This is due to the fact that a rotation of the sphere does not change distances on the sphere.
The next step is to consider these particular circles of radius $r$ around $N$. The best way to see this is to draw a sketch of the intersection of the sphere with the $yz$-plane, i.e. a circle of radius 1 around the origin. We can also see to points $P_1, P_2$ on the circle at distance $r$ from $N$. (I'll assume $r$ is small enough for this to make sense) These points should be such that $P_1$ is obtained from $P_2$ by reflecting around the $z$-axis. The spherical coordinate $\theta$ corresponding to these points can now easily be obtained: this angle is given as an angle wrt the $z$-axis, while $r$ is precisely the arc length corresponding to this angle on a circle with radius 1. This leads to the conclusion that $r=\theta$ in this special case.
The next thing to notice, is that all points at distance $r$ from $N$ have the same $z$-coordinate. By projecting either $P_1$ or $P_2$ at this $z$-coordinate (which again can be done insightful in the same picture), one finds this height to be $z_0=\cos(\theta)=\cos(r)$. Moreover, yet again from considering the same picture, the required circle is a circle parallel to the $xy$-plane with center $(0,0,z_0)$ and radius $\sin(\theta)=\sin(r)$.
We may thus conclude that the points at distance $r$ from $N$ are of the form $(\sin(r)\cos(\phi),\sin(r)\sin(\phi),cos(r))$ with $\phi\in[0,2\pi)$. Using the rotation matrix $R$ from before, we thus get the set of points satisfying your original question:
$$\{R^{-1}(\sin(r)\cos(\phi),\sin(r)\sin(\phi),cos(r))\mid \phi\in[0,2\pi)\}.$$
$\DeclareMathOperator{\proj}{proj}$This can be handled easily using vector algebra. Let $A$ and $B$ be the specified points on the unit sphere, viewed as vectors in space, and assume $A \neq \pm B$ and $\cos r < A \cdot B$ ($A$ is outside the circle $C$ of radius $r$ centered at $B$). Let $U_{1}$ be the unit vector orthogonal to $B$ and tangent to the "short" great circle arc from $B$ to $A$, and let $U_{2} = B \times U_{1}$.
Theorem: If
$$
\theta = \arccos\left[\frac{(A \cdot B)\tan r}{A \cdot U_{1}}\right],
$$
then the points $T_{1}$ and $T_{2}$ are
$$
(\cos r)B + \sin r\bigl((\cos\theta) U_{1} \pm (\sin\theta) U_{2}\bigr).
$$
Computational note:
- $U_{1}$ is the second vector obtained from the Gram-Schmidt algorithm on the ordered basis $\{B, A\}$; that is,
$$
U_{1} = \frac{A - \proj_{B}A}{\|A - \proj_{B}A\|}.
$$
Proof of Theorem: Every point on the circle $C$ has the form
$$
T = (\cos r)B + \sin r\bigl((\cos\theta) U_{1} \pm (\sin\theta) U_{2}\bigr)
\tag{1}
$$
for some real $\theta$.
If $T$ is an arbitrary point of $C$, let $A' = A - T$ and $B' = B - T$ denote the displacements to $A$ and $B$. The projections of these vectors to the tangent plane at $T$ are
$$
A'' = A' - (A' \cdot T)T,\qquad
B'' = B' - (B' \cdot T)T,
$$
and $T$ is a point of tangency if and only if $A'' \perp B''$. Expanding the dot product, using the fact that $T$ is a unit vector $(T \cdot T = 1$) lying on the circle of radius $r$ centered at $B$ ($B \cdot T = \cos r$), gives (omitting most of the algebra)
$$
0 = A'' \cdot B''
= A' \cdot B' - (A' \cdot T)(B' \cdot T)
= A \cdot B - (\cos r)A \cdot T,
$$
or
$$
A \cdot T = \frac{A \cdot B}{\cos r}.
\tag{2}
$$
Substituting (1) in (2), using the fact that $A \cdot U_{2} = 0$, and solving for $\theta$ gives the value in the theorem.
Best Answer
Construct a parametrization of the circle $d$ in $3$ stages of increasing generality.
1. Suppose that $B$ is the north pole (colatitude $\varphi = 0$).
The circle is the parallel (line of constant colatitude $\varphi = \alpha$), which has $z = \cos \alpha$ and $r = \sin \alpha$; hence, it is parametrized in rectangular coordinates by $t \mapsto \big( x(t), y(t), z(t) \big)$, where $$ \left\{ \begin{align} x &= \sin \alpha \cos t \\ y &= \sin \alpha \sin t \\ z &= \cos \alpha \end{align} \right. \qquad \text{for } 0 \le t < 2\pi. $$
2. Suppose that $B$ is on the prime meridian (longitude $\theta = 0$) but colatitude has some value $\varphi = \beta$, where $0 < \beta \le \pi$.
We take the coordinates $(x, y, z)$ of the circle with center at the north pole (from 1.) and rotate them through an angle $\beta$ along the great circle that includes the prime meridian. Using standard formulas for rotating coordinates, we have $$ \begin{bmatrix} x \\ y \\ z \end{bmatrix} \mapsto \begin{bmatrix} x \cos \beta + z \sin \beta \\ y \\ -x \sin \beta + z \cos \beta \end{bmatrix}. $$
3. Suppose that $B$ is anywhere on the sphere with spherical coordinates $(\varphi, \theta) = (\beta, \gamma)$. Take the resulting coordinates from 2. and rotate them about the polar axis of the sphere through an angle of $\gamma$: $$ \begin{bmatrix} x \\ y \\ z \end{bmatrix} \mapsto \begin{bmatrix} x \cos \gamma - y \sin \gamma \\ x \sin \gamma + y \cos \gamma \\ z \end{bmatrix}. $$
Putting these together, the general case for a circle $d$ with center at $(\varphi, \theta) = (\beta, \gamma)$ and central angle $\alpha$ (that determines the radius) has rectangular coordinates: $$ \left\{ \begin{align} x &= \phantom{-}( \sin \alpha \cos \beta \cos \gamma ) \cos t + ( \sin \alpha \sin \gamma ) \sin t - (\cos \alpha \sin \beta \cos \gamma ) \\ y &= -( \sin \alpha \cos \beta \sin \gamma ) \cos t + ( \sin \alpha \cos \gamma ) \sin t + (\cos \alpha \sin \beta \sin \gamma ) \\ z &= \phantom{-}( \sin \alpha \sin \beta ) \cos t + \cos \alpha \cos \beta. \end{align} \right. $$