How to find 4-velocity components in a perturbed metric? Beginning with a metric with small perturbations
\begin{eqnarray}
g_{00} &=& 1-2\frac{U}{c^2} + \mathcal{O}(c^{-4}) \\
g_{0i} &=& \mathcal{O}(c^{-3}) \\
g_{ij} &=& -\delta_{ij}\left(1+2\frac{U}{c^2} + \mathcal{O}(c^{-4})\right)
\end{eqnarray}
I'm trying to show that the components of the 4-velocity vector $u^0$ and $u^i$ can be written as 
\begin{eqnarray}
u^0 &=& 1 + \frac{1}{2}\frac{v^2}{c^2} + \frac{U}{c^2} + \mathcal{O}(c^{-4}) \\
u^i &=& \frac{v^i}{c}\left(1 + \frac{1}{2}\frac{v^2}{c^2} + \frac{U}{c^2}\right) + \mathcal{O}(c^{-5})
\end{eqnarray}
I want to explicitly keep powers of $c$ in the calculation to help learn about different orders of smallness in PPN. The method I am attempting to use is to write $u^\mu u_\mu = u^0 u_0 + u^i u_i = +1$. This may be rearranged as
\begin{eqnarray}
u^0 u^0 g_{00} &=& 1 - u^i u^j g_{ij} - 2u^0 u^i g_{0i} \\
\left(u^0\right)^2 \left[ 1-2\frac{U}{c^2}\right] &=& 1 + \left(\frac{dx^i}{ds}\right)^2\left[ 1+2\frac{U}{c^2} \right]
\end{eqnarray}
I dropped the $g_{0i}$ term because it is $\mathcal{O}(c^{-3})$. I found that if I write $\frac{dx^i}{ds} = \frac{1}{c}\frac{dx^i}{dt}\frac{dt}{ds} = \frac{v}{c}$, solve the equation above for $u^0 = \sqrt\frac{...}{...}$, and then taylor expand about the small terms $\left(\frac{v}{c}\right)^2$ and $\frac{U}{c^2}$, I get the correct form of $u^0$.
My questions are:


*

*Is this the best approach to solve for $u^0$ and $u^i$, or is there a better/more highly recommended method to use?

*I don't understand if or why it is correct to write $\left(\frac{dx^i}{ds}\right)^2 = \left(\frac{v}{c}\right)^2$. The term $\frac{dx^i}{ds}$ should have no units (like $\frac{v}{c}$), but when I use the chain rule to expand $\frac{dx^i}{ds}$, shouldn't there be a factor of c next to dt in the numerator above as well? Could it have something to do with thinking of ds as proper time, and then writing $dt\approx ds$ for slow moving objects?

*Also, I read that metric components of $g_{0i}$ will have odd powers of $c^{-1}$, while metric components of $g_{00}$ and $g_{ij}$ will have even powers of $c^{-1}$, but have never seen a justification for that claim.

*Ultimately my goal is to understand the PPN formalism. Are there any textbooks or papers where the author works through expanding a perturbed metric while keeping factors of $c^{-1}$ explicitly in the calculation to show different orders of the small terms (e.g. does not set $c=1$)? So far, the report I linked below is the only paper I have ever seen that does not work in units where $c=1$.
My question is about how to arrive at eq. (15) in this link (the link will automatically download a pdf of the paper). 
If you don't want to click on a link that automatically downloads a paper, the title of the paper is "The post-Newtonian formalism" by Rene Michelsen.
 A: Express four-velocities $u^\alpha$ and coordinate velocities $v^i$ according to
$$ u^\alpha = \frac{dx^\alpha}{d\tau}, \quad v^i = \frac{dx^i}{dt}, $$
where proper time and coordinate time is given by $\tau$ and $t$ respectively. Hereafter, Latin ($i,j = 1,2,3$) and Greek ($\alpha,\beta = 0,1,2,3$) indices account for spatial and spacetime variables respectively.
I will answer your question $(1)$ in terms of arbitrary metric tensor components, where I will outline the procedure you need to follow in order to get the desired results. 
A timelike particle obeys the normalisation condition given by 
$$ g_{\alpha \beta} u^\alpha u^\beta = -c^2,  \tag{1} \label{eqn1}$$ 
where $c$ is the speed of light. Eq.  (\ref{eqn1}) is given explicitly by 
$$ \left( u^0 \right)^2 \left( g_{00} + 2 g_{0i} \frac{u^i}{u^0} + g_{ij} \frac{u^i}{u^0} \frac{u^j}{u^0} \right) = -c^2. \tag{2} \label{eqn2} $$
Now, we know $x^0 = ct$. Hence $d x^0 = cdt. $ Using this we can express $u^i/u^0$ as the following 
$$ \frac{u^i}{u^0} = \frac{dx^i}{d\tau} \left(\frac{d x^0}{d\tau} \right)^{-1} = \frac{dx^i}{dt} \frac{dt}{d\tau} \left( \frac{dx^0}{dt} \frac{dt}{d\tau} \right)^{-1} = \frac{v^i}{c}.   \tag{3} \label{eqn3}$$ 
(Note: I hope that addresses your confusion regarding question $2$. ) Returning to Eq. (\ref{eqn2}) we have 
$$ u^0 = c \left[ - g_{00} - 2g_{0i} \left( \frac{v^i}{c} \right) - g_{ij} \left(\frac{v^i}{c}\right) \left( \frac{v^j}{c} \right) \right]^{-1/2}.$$
Substituting in the appropriate metric tensor components and expanding the square root to PN order will give the desired $u^0$. Consulting Eq. (\ref{eqn3}) gives the appropriate expression for the spatial part of the four-velocity namely 
$$ u^i = \frac{v^i}{c} u^0. $$
So to answer your question $(1)$, I would consult the normalisation condition to easily obtain the components of $u^\alpha$. By following the method outlined above you will obtain the expressions you have quoted for $u^0$ and $u^i$. 
I should note that there is a big difference between the post-Newtonian and parametrised post-Newtonian formalisms which I would advise you familiarise yourself with. I assume you mean you want to understand the PN formalism in your question $(4)$. Works by Chandrasekhar on the PN formalism generally retain the factor of $1/c$ for convenience. However, you should know that this is simply a matter of correctly determining the appropriate dimensions and reinstating. Being able to do this is more beneficial in the longrun. 
