Convert both points to Cartesian coordinates, calculate the cross product a×b. The direction of a×b is the axis and the magnitude is

|a×b| = |a| |b| sin(α)

where α is the angle between a,b (i.e. the rotation angle which maps a to b).

Conversion to Cartesian coordinates depends upon the interpretation of θ,φ.

If θ is the angle relative to the positive Z axis, θ∈[0,π), and φ is 0 on the X-Z plane and π/2 on the Y-Z plane, φ∈[0,2π), then

x = r sin(θ) cos(φ)

y = r sin(θ) sin(φ)

z = r cos(θ)

Different conventions for θ,φ will produce similar equations but with x,y,z permuted, θ,φ permuted, and/or sin,cos permuted.

Using r=1 => |a|=|b|=1 => |a×b| = sin(α)

[x1,y1,z1] × [x2,y2,z2] = [y1 z2 - y2 z1, x2 z1 - x1 z2, x1 y2 - x2 y1]

Caveat: any value of sin(α) gives you two possible values for α. One will rotate a to b, the other will rotate a to -b (so the vectors will be parallel but in opposite directions). You can use the dot product:

[x1,y1,z1] · [x2,y2,z2] = x1 x2 + y1 y2 + z1 z2

to determine which to use. If a·b is positive, the angle between them is less than π/2, and α should be the principal arcsine (i.e. the value which asin() will return). If a·b is negative, the angle is greater than π/2 so you need to use π-asin(α).

Another caveat: magnitude is always positive so sin(α) and α will always be positive. If you swap the order of the vectors, b×a = -(a×b), i.e. you'll still get a positive rotation but the axis will be pointing in the opposite direction.

Final caveat: if the vectors are exactly perpendicular, floating-point rounding error can result in |a×b| being slightly larger than 1, which will result in asin() reporting a domain error (because sin(x)∈[-1,1]), so you might want to check for that.