limit of integration is a vector, but I can rewrite the integrand as a scalar. What do I do? I have to solve the following integral $$\Phi(\boldsymbol r) =- \int_{\boldsymbol r_0}^{\boldsymbol r}\boldsymbol{E}(\tilde{\boldsymbol  r}) \cdot d\boldsymbol s$$
This comes from Electrostatics (I'm computing the electric potential of a charge Q at the origin, with reference point $\boldsymbol r$ being asymptotic infinity)
We know that in general for a charge at the origin$$\boldsymbol E(\boldsymbol r) = \frac{1}{4\pi\epsilon_0}\frac{Q}{r^2}\hat{\boldsymbol e}_r$$
My trial
Clearly we use polar coordinates $\{r, \theta, \phi\}$ because of the spherical symmetry.Then the line element becomes $d\boldsymbol s = dr\hat{\boldsymbol e}_r+rd\theta\hat{\boldsymbol e}_{\theta}+r\sin{\theta}d\phi\hat{\boldsymbol e}_{\phi}$ so then 
$$-\int_{\boldsymbol r_0}^{\boldsymbol r}\boldsymbol{E}(\tilde{\boldsymbol  r}) \cdot d\boldsymbol s = -\frac{Q}{4\pi\epsilon_0}\int_{\boldsymbol  r_0}^{\boldsymbol r}\frac{1}{\tilde{r}^2}\hat{\boldsymbol e}_{\tilde r}\cdot d\boldsymbol s =-\frac{Q}{4\pi\epsilon_0}\int_{\boldsymbol r_0}^{\boldsymbol r}\frac{1}{\tilde{r}^2}d\tilde{r}$$
Because the spherical polar basis is an orthogonal basis, hence they are all perpendicular to each other and the result of the dot product is $1$ coming from the radial vector. But now I have a scalar integral and I will have to evaluate the result at a vector?

Can someone explain to me (possibly explaining each passage and with some strong mathematical rigor) how someone would finish this calculation?

 A: Note that in general we have
$$\Phi(\vec r)=\Phi(\vec r_0)-\int_{\vec r_0}^{\vec r} \vec E(\vec r')\cdot d\vec \ell'$$
The notation for the line integral can be more explicitly written by parametrizing the curve traversed from $\vec r_0$ to $\vec r$.  Let the curve be described parametrically as $\vec r(t)$ where $\vec r=\vec r(1)$ and $\vec r_0=\vec r(0)$.  Then, we can write
$$\int_{\vec r_0}^{\vec r} \vec E(\vec r')\cdot d\vec \ell'=\int_0^1 \vec E(r(t),\theta(t),\phi(t))\cdot \left(\hat r\frac{dr(t)}{dt}+\hat r\frac{d\theta(t)}{dt}+\hat r\frac{d\phi(t)}{dt}\right)\,dt \tag 1$$
If $\vec E$ has only a radial component that depends on $r(t)$ only, then $(1)$ becomes
$$\begin{align}
\int_{\vec r_0}^{\vec r} \vec E(\vec r')\cdot d\vec \ell'&=\int_0^1 E_r(r(t))\frac{dr(t)}{dt}\,dt\\\\
&=\int_{r(0)}^{r(1)}E_r(r)\,dr \tag 2
\end{align}$$
Note that $r(1)=|\vec r|$ and $r(0)=|\vec r_0|$.  Thus, we can write the integral in $(2)$ as 
$$\int_{\vec r_0}^{\vec r} \vec E(\vec r')\cdot d\vec \ell'=\int_{|\vec r_0|}^{|\vec r|} E_r(r)\,dr \tag3$$
Finally, using $E_r=\frac{Q}{4\pi \epsilon_0 r^2}$ we find from $(3)$ that 
$$\Phi(\vec r)=\Phi(\vec r_0)+\frac{Q}{4\pi \epsilon_0 |\vec r|}-\frac{Q}{4\pi \epsilon_0 |\vec r_0|}$$
Letting $|\vec r_0|\to \infty$ and setting $\Phi(\infty)=0$, we obtain the coveted result
$$\Phi(\vec r)\frac{Q}{4\pi \epsilon_0 |\vec r|}$$
