I think
s2 = Solve[eqs2[r, \[Theta], \[Phi]] == R^2, r][[2]] // Simplify;
is not appropriate. We are looking for an intersection of the paraboloid and the sphere. All these points are elements of the spherical surface and therefore have a r-value of R, so we must not look for r.
The intersection being a line means that it is a one-dimensional entity, depending only on ONE pararmeter. So I think fixing r to the radius of the sphere we should ask for theta as function of phi or phi as function of theta and insert this in the representation of a point:
For example
s3 = Solve[eqs1[R, \[Theta], \[Phi]] == R^2, \[Phi]]
Length[s3]
plot3 = ParametricPlot3D[{x1[R, \[Theta], \[Phi]],
x2[R, \[Theta], \[Phi]], x3[R, \[Theta], \[Phi]]} /.
s3[[5, 1]], {\[Theta], 0, Pi},
PlotStyle -> {Thick, Blue}]
Show[plot2, plot3]
or
s4 = Solve[eqs1[R, \[Theta], \[Phi]] == R^2, \[Theta]]
Length[s4]
plot4 = ParametricPlot3D[{x1[R, \[Theta], \[Phi]],
x2[R, \[Theta], \[Phi]], x3[R, \[Theta], \[Phi]]} /.
s4[[6, 1]], {\[Phi], 0, 2 Pi},
PlotStyle -> {Thick, Blue}]
Show[plot2, plot4]
As s3 and s4 yield several solutions it seems that the intersection is given piecewise. But put together it seems that both solutions are equivalent, as it should be:
plot5= ParametricPlot3D[
{x1[R, \[Theta], \[Phi]], x2[R, \[Theta], \[Phi]],
x3[R, \[Theta], \[Phi]]} /. # & /@ Flatten[s3],
{\[Theta], 0, Pi},
PlotStyle -> {Thick, Blue}]
plot6= ParametricPlot3D[{x1[R, \[Theta], \[Phi]], x2[R, \[Theta], \[Phi]],
x3[R, \[Theta], \[Phi]]} /. # & /@ Flatten[Drop[s4, 3]],
{\[Phi], 0, 2 Pi}, PlotStyle -> {Thick, Blue}]
Show[plot2, plot5, plot6]