You can use the 2D case to compute this: take one circle of radius $r_1$ centered at the origin, and a second circle of radius $r_2$ centered at $(d, 0)$. Then you get two points of intersection (computed using Wolfram Alpha):
$$\left(\frac{d^2+r_1^2-r_2^2}{2d},
\pm\sqrt{r_1^2-\frac{(d^2+r_1^2-r_2^2)^2}{4d^2}}\right)$$
In your 3D setup, you set $d$ to the distance between the two centers of the spheres, then do the above computation. The first coordinate of the result represents the distance from the point $p_1$ in the direction of $p_2$. The second coordinate, choosen positively, represents the radius of the circle of intersection.
Because of the radial symmetry of the problem, it’s not substantially different from the two-dimensional one. Indeed, if you restrict your attention to a plane that passes through the centers of both spheres, then you have precisely that two dimensional problem of finding the intersections of two circle. By symmetry, one can see that the intersection of the two spheres lies in a plane perpendicular to the line joining their centers, therefore once you have the solutions to the restricted circle intersection problem, rotating them around the line joining the sphere centers produces the other sphere intersection points. Actually, because of the symmetry of the 2-d problem, you really only need one intersection point since the other circle intersection, if any, will generate the same circle in 3-d.
The real complications arise when you try to describe this intersection circle with equations. A single implicit three-dimensional Cartesian equation generally describes a surface, not a curve; you need two such equations to describe a curve. Using the equations of the two spheres takes you right back to step one, so you’d need to find a different pair of equations to describe this same circle. A reasonably descriptive choice might be the intersection of the plane on which the circle lies with a rectangular cylinder that’s perpendicular to that plane, but even that’s going to be rather opaque.
A more immediately understandable description is via a parametric vector equation of the form $$\mathbf c+r\mathbf u\cos t+r\mathbf v\sin t.$$ (This is where those trigonometric functions that you mention are likely coming into the picture.) Here, $\mathbf c$ is the circle’s center, $r$ its radius and $\mathbf u$ and $\mathbf v$ an orthogonal pair of orthogonal vectors that are parallel to the circle’s plane. Geometrically, this is just the parametric equation $r\cos t(1,0,0)+r\sin t(0,1,0)$ of a circle in the $x$-$y$ plane of radius $r$ centered at the origin that’s been rotated and translated into position. To tie this to the first paragraph, $\mathbf c+r \mathbf u$ is a solution to a reduced planar circle intersection problem (in the plane through the sphere centers and $\mathbf c+\mathbf u$), and the above formula describes all possible rotations of this point around the line through the sphere centers. I’m sure that one can find other parameterizations of the intersection circle, but I think this one is the most transparent.
Best Answer
HINT.
See below a section of the spheres, passing through their centers $A$ and $B$. They intersect each other orthogonally if radii $AC$ and $BC$ are perpendicular.
It follows that $ABC$ is a right triangle with legs $r_1$ and $r_2$. And its altitude $CH$ is the radius of the intersection circle.