I am updating this answer to try to address the edited version of the question.
A nice thing about the conventional $(x,y,z)$ Cartesian coordinates is everything works the same way. In fact, everything works so much the same way using the same three coordinates in the same way all the time in Cartesian coordinates--points in space, vectors between points, field vectors--that it may be surprising that in just about any other coordinate system different things sometimes work differently from each other.
In Cartesian coordinates, you can identify the components of a point's coordinates with the components of the point's position vector.
And you can get the vector sum of two of those vectors by adding the coordinates:
$$
\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}
+ \begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix}
= \begin{pmatrix} 1+0 \\ 0+1 \\ 0+0 \end{pmatrix}
= \begin{pmatrix} 1 \\ 1 \\ 0 \end{pmatrix}.
$$
But if you try to describe a vectors by treating them as position vectors and using the spherical coordinates of the points whose positions are given by the vectors, the left side of the equation above becomes
$$
\begin{pmatrix} 1 \\ \pi/2 \\ 0 \end{pmatrix}
+ \begin{pmatrix} 1 \\ \pi/2 \\ \pi/2 \end{pmatrix},
$$
while the right-hand side of the equation becomes
$$
\begin{pmatrix} \sqrt2 \\ \pi/2 \\ \pi/4 \end{pmatrix}.
$$
So spherical coordinates are really awkward to do vector operations in. You just can't do it like you can with Cartesian coordinates.
Another thing that's going on with these formulas is we have to make sense of things like $\frac\partial{\partial\theta}.$
If you have something that can vary as $r,$ $\theta,$ and $\varphi$ vary, then
$\frac\partial{\partial\theta}$ is looking for what happens when we hold $r$ and $\phi$ constant but vary $\theta.$
What happens to the location defined by the coordinates $(r, \theta, \varphi)$
if we hold $r$ and $\varphi$ constant but increase $\theta$ is that we move in a "downward" direction perpendicular to the position vector.
Given a starting point $(r, \theta, \varphi)$ for this motion, we define
$\hat{\boldsymbol \theta}$ as the unit vector in that direction.
Similarly, there is a unit vector $\hat{\boldsymbol \varphi}$ in the direction we would travel if we hold $r$ and $\theta$ constant and increase $\varphi.$
There is also a unit vector $\hat{\mathbf r}$ in the direction of constant $\theta$ and $\varphi$ but increasing $r.$
I name these vectors $\hat{\mathbf r},$ $\hat{\boldsymbol \theta},$ and
$\hat{\boldsymbol \varphi}$ because that's how they're named on the first Wikipedia page the question linked to. These are the same vectors named
$\hat e_r,$ $\hat e_\theta,$ and $\hat e_\varphi$ elsewhere.
Unlike the unit vectors $\hat{\mathbf x},$ $\hat{\mathbf y},$ and
$\hat{\mathbf z},$ which are always parallel to the $x,$ $y,$ and $z$ axes,
the unit vectors $\hat{\mathbf r},$ $\hat{\boldsymbol \theta},$ and
$\hat{\boldsymbol \varphi}$ depend on the coordinates of a point
$(r, \theta, \varphi)$ relative to which those vectors are defined.
Specifically,
$$\begin{align}
\hat{\mathbf r} &= \sin\theta \cos\varphi \hat{\mathbf x}
+ \sin\theta\sin\varphi \hat{\mathbf y} + \cos\theta \hat{\mathbf z}, \\
\hat{\boldsymbol \theta} &= \cos\theta \cos\varphi \hat{\mathbf x}
+ \cos\theta \sin\varphi \hat{\mathbf y} - \sin\theta \hat{\mathbf z}, \\
\hat{\boldsymbol \varphi} &= -\sin\varphi \hat{\mathbf x} + \cos\varphi \hat{\mathbf y}.
\end{align}
$$
You can visualize a physical model for these vectors by supposing that the Earth is spherical and is centered at the origin of the coordinate system. Then at any point on the surface of the Earth, $r$ is the Earth's radius,
$\theta$ is the co-latitude (like latitude, but it goes from zero at the north pole to $180$ degrees at the south pole), and $\varphi$ is the longitude.
Then $\hat{\boldsymbol \theta}$ is in the direction of increasing co-latitude, so it points due south;
$\hat{\boldsymbol \varphi}$ is in the direction of increasing longitude, so it points due east;
and $\hat{\mathbf r}$ is in the direction of increasing radius, so it points directly away from the center.
In particular, $\hat{\mathbf r}$ is a unit vector in the exact same direction as the position vector $\vec r,$
and the other two unit vectors are orthogonal to $\hat{\mathbf r}$ and to each other.
This clearly does not work as the basis of a set of position vectors for all the points in a three-dimensional space. If you wanted the vector describing the direction and distance from point $P$ to point $Q,$ do you use the vectors
$\hat{\mathbf r},$ $\hat{\boldsymbol \theta},$ and
$\hat{\boldsymbol \varphi}$ according to how they are defined relative to $P$ or according to how they are defined relative to $Q$?
The answer is (usually) neither, you don't typically use those vectors for that kind of calculation at all.
What those vectors are good for is as an orthonormal basis for things that occur locally at a point: a vector field at that point, the gradient of that field, the divergence of that field, and so forth.
Now if you have a vector field with the value $\vec A$ at some point with spherical coordinates
$(r, \theta, \varphi),$ then we can break that vector down into orthogonal components exactly as you do:
\begin{align}
A_r &= \vec A\cdot\hat{\mathbf r}, \\
A_\theta &= \vec A\cdot\hat{\boldsymbol \theta}, \\
A_\varphi &= \vec A\cdot\hat{\boldsymbol \varphi}.
\end{align}
Now consider the case where $\vec A = \vec r.$ Then $\vec A$ is in the exact same direction as
$\hat{\mathbf r},$ and therefore
$$ A_r = \vec A\cdot\hat{\mathbf r} = \lVert \vec A\rVert \lVert \hat{\mathbf r}\rVert = r. $$
On the other hand, $\vec A$ is orthogonal to each of the vectors
$\hat{\boldsymbol \theta}$ and $\hat{\boldsymbol \varphi},$ so
$$ \vec A\cdot\hat{\boldsymbol \theta} = \vec A\cdot\hat{\boldsymbol \varphi} = 0.$$
In summary, $A_r = r$ and $A_\theta = A_\varphi = 0.$
Moreover, no matter which way you go, $A_\theta$ and $A_\varphi$ both remain zero, so
$$
\frac{\partial}{\partial\theta}(A_\theta\sin\theta)
= \frac{\partial A_\varphi}{\partial\varphi} = 0.
$$
If you are measuring the field at the point $(1,1,1),$ then $\hat{\mathbf r}$ at that point is
$\hat{\mathbf r} = \frac1{\sqrt3}\hat{\mathbf x} +
\frac1{\sqrt3}\hat{\mathbf y} + \frac1{\sqrt3}\hat{\mathbf y},$
and therefore
$\sqrt3\hat{\mathbf r} = \hat{\mathbf x} + \hat{\mathbf y} + \hat{\mathbf z},$ or in the usual $(x,y,z)$ Cartesian system, the vector $(1,1,1).$
Plugging in the correct values for all three components, the divergence formula gives the result $3.$
Now if you had a field that was not simply a radial field always going directly outward from the center or directly inward toward the center, that is, if there were some field vectors making other angles with the radial vectors, then you would see some non-zero values of $A_\theta$ and $A_\phi,$ and there would be a reason you needed the rest of the formula for divergence. You just happen to have chosen an example in which you have eliminated those components.
If you want a more interesting example, you might try a vector field such as
$\vec A = x \hat{\mathbf x},$ which has non-zero derivatives in all three directions at most points.
Best Answer
Making some progress with the algebra, for the problem in spherical coordinates I've managed to obtain the expression:
$$ \tan \left ( \varphi_\infty - \varphi_0 \right ) = \frac{r_0 \sin \theta_0 \left . \frac{\partial P_\varphi}{\partial P_r} \right |_{t=0}}{ \sin \theta_0 + r_0 \cos \theta_0 \left . \frac{\partial P_\theta}{\partial P_r} \right |_{t=0}} $$
If we restrict our particle to equatorial motion, taking $\theta_0 = \pi / 2$, this formula recovers the formula in the OP for polar coordinates:
$$ \tan \left ( \varphi_\infty - \varphi_0 \right ) = r_0 \left . \frac{\partial P_\varphi}{\partial P_r} \right |_{t=0}$$
After a significant amount more algebra, I believe I have also arrived at a formula for $\theta_\infty$:
$$ \cos \theta_\infty = \frac{ \cos \theta_0 - r_0 \sin \theta_0 \left . \frac{\partial P_\theta}{\partial P_r} \right |_{t=0} } { \sqrt { 1 + \left ( r_0 \left . \frac{\partial P_\theta}{\partial P_r} \right |_{t=0} \right ) ^2 + \left ( r_0 \sin \theta_0 \left . \frac{\partial P_\varphi}{\partial P_r} \right |_{t=0} \right ) ^2 } } $$
As of yet I haven't checked the correctness of these formulas. Has anyone come across such formulas, e.g. when computing the Jacobian for a translation in spherical coordinates?