Gravitational force that a (thin) disk exerts on a particle 
Calculate gravitational force that a (thin) disk with mass $M$ uniformly distributed over disk radius $a$ exerts on a small particle with mass $m$. The particle is at distance $x$ above the disk center and perpendicular to the disk surface (see figure here).

I calculated the gravitational force in two different ways and got two different results. Below I give detailed procedures for both approaches. Can you please help me understand why am I getting two different results? I know there are other ways to solve this problem, e.g. by integrating over thin concentric rings, but I want to understand what am I missing by integrating directly over infinitesimally small mass $dm$. I have an idea what might be the problem, but I am not sure if I am missing some bigger picture here - see below the boldfaced question.
Apologies for such a long question - I wanted to include every detail of my attempt.

Gravitational force: Brief introduction
Gravitational force that a particle with mass $m_2$ exerts on a particle with mass $m_1$ is
$$\vec{F} = G \frac{m_1 m_2}{|\vec{r}|^2} \hat{r}$$
where $G$ is the gravitational constant, $\vec{r} = \vec{r}_2 - \vec{r}_1$ is the vector from particle 1 to particle 2, $\vec{r}_1$ and $\vec{r}_2$ are position vectors of particles 1 and 2, respectively, and $\hat{r} = \vec{r} / |\vec{r}|$ is the unit vector. Note that gravitational force is attractive, i.e. it always points from one particle to the other.

Approach I: Neglects disk thickness
Here I completely neglect the effect of the disk thickness $t$ on the total gravitational force. Let the disk be centered at the origin with its surface in the $ij$ plane, and let the particle be on the $k$ axis. The total gravitational force is
$$\vec{F} = \iiint_D G \frac{m dm}{|\vec{r}|^2} \hat{r} = \iiint_D G \frac{m (\lambda dV)}{|\vec{r}|^2} \frac{\vec{r}}{|\vec{r}|} = G m \lambda \iiint_D \frac{\vec{r}}{|\vec{r}|^3} \rho d\rho d\varphi dz$$
where $dm$ and $dV$ are infinitesimally small mass and volume on the uniform disk, respectively, $\lambda$ is the disk density, and vector $\vec{r}$ is
$$\vec{r} = \vec{r}_{dm} - \vec{r}_m = (\rho \cos\varphi, \rho \sin\varphi, 0) - (0, 0, x) = (\rho \cos\varphi, \rho \sin\varphi, -x)$$
The total gravitational force is then
$$\vec{F} = G m \lambda \int_{-t/2}^{t/2} \int_{0}^{2 \pi} \int_{0}^{a} \frac{\rho}{( \rho^2 + x^2 )^{3/2}} (\rho \cos\varphi \hat{\imath} + \rho \sin\varphi \hat{\jmath} - x \hat{k}) d\rho d\varphi dz$$
Obviously the $i$ and $j$ components will integrate to zero and we are left only with $k$-axis component
$$F_k = G m \lambda \cdot 2\pi \cdot t \int_{0}^{a} \frac{-\rho x}{( \rho^2 + x^2 )^{3/2}} d\rho = 2 G \frac{m M}{a^2} \int_{0}^{a} \frac{-\rho x}{( \rho^2 + x^2 )^{3/2}} d\rho$$
where $M = \lambda \cdot (a^2 \pi \cdot t)$ is the total mass of the disk. The integral is simply solved by substitution $u^2 = \rho^2 + x^2$ and $u du = \rho d\rho$
$$F_k = 2 G \frac{m M}{a^2} \Bigl( \frac{x}{\sqrt{a^2 + x^2}} - 1 \Bigr)$$
This result agrees with some available solutions done by integrating over thin rings. The sign also seems to be correct - for particle position $x > 0$ the gravitational force exerted by the disk on the particle is negative $F_k < 0$, which means it points towards the disk (origin).

Approach II: Does not neglect disk thickness
The procedure here is similar to the one presented in Approach I. The only difference is that I do not neglect effect of the disk thickness $t$ on the gravitational force
$$\vec{r} = \vec{r}_{dm} - \vec{r}_m = (\rho \cos\varphi, \rho \sin\varphi, z) - (0, 0, x) = (\rho \cos\varphi, \rho \sin\varphi, z-x)$$
Following the same procedure as in the Approach I, the gravitational force in $k$ axis is
$$F_k = 2 G \frac{m M}{a^2} \frac{1}{t} \int_{-t/2}^{t/2} \int_{0}^{a} \frac{\rho (z - x)}{( \rho^2 + (z - x)^2 )^{3/2}} d\rho dz$$
The integral over $d\rho$ is solved by substitution $u^2 = \rho^2 + (z-x)^2$ and $u du = \rho d\rho$
$$F_k = 2 G \frac{m M}{a^2} \frac{1}{t} \int_{-t/2}^{t/2} \Bigl( \frac{(z-x)}{\sqrt{(z - x)^2}} - \frac{(z-x)}{\sqrt{a^2 + (z - x)^2}} \Bigr) dz =$$
$$= 2 G \frac{m M}{a^2} \frac{1}{t} \int_{-t/2}^{t/2} \Bigl( 1 - \frac{(z - x)}{\sqrt{a^2 + (z - x)^2}} \Bigr) dz$$
The only thing that comes in mind is that $\sqrt{(z-x)^2}$ resolves into $\pm (z-x)$, where $+(z-x)$ gives wrong end result, while $-(z-x)$ gives correct end result. Is that the fallacy in this procedure - when resolving square roots I should always solve for both positive and negative signs and later check which one does physically make sense? By taking $\sqrt{(z-x)^2} = +(z-x)$ the above expression becomes
$$F_k = 2 G \frac{m M}{a^2} \frac{1}{t} \left. \Bigl( z - \sqrt{a^2 + (z - x)^2} \Bigr) \right|_{-t/2}^{t/2} =$$
$$= 2 G \frac{m M}{a^2} \Bigl( 1 + \frac{1}{t} \sqrt{a^2 + (t/2 + x)^2} - \frac{1}{t} \sqrt{a^2 + (t/2 - x)^2} \Bigr)$$
The difference of two square roots can be written as
$$\sqrt{a^2 + (t/2 + x)^2} - \sqrt{a^2 + (t/2 - x)^2} = \frac{(t/2 + x)^2 - (t/2 - x)^2}{\sqrt{a^2 + (t/2 + x)^2} + \sqrt{a^2 + (t/2 - x)^2}}$$
The expression for $F_k$ is now simplified into
$$F_k = 2 G \frac{m M}{a^2} \Bigl( 1 + \frac{2x}{\sqrt{a^2 + (t/2 + x)^2} + \sqrt{a^2 + (t/2 - x)^2}} \Bigr)$$
and for (very) small thickness $t \to 0$ it is further simplified into
$$F_k = 2 G \frac{m M}{a^2} \Bigl( \frac{x}{\sqrt{a^2 + x^2}} + 1 \Bigr)$$
This result is not equal to the result obtained in Approach I. The sign check also indicates the end result is not correct - for particle position $x > 0$ the force exerted by the disk on the particle is positive $F_k > 0$ which indicates that the disk pushes the particle away, which cannot be true. However, by taking $\sqrt{(z-x)^2} = -(z-x)$ the $+1$ term in the parenthesis becomes $-1$ and the end result equals the one I got in the Approach I. Is this really the root cause or am I missing something bigger here?
Thanks for reading!
 A: As I indicated in my comment, for any real or complex number $\ A\ $,$\ \sqrt{A}\ $ is uniquely defined to be the principal square root of $\ A\ $. For a non-negative real number $\ A\ $, the principal square root is the one which is non-negative.  Thus, in the integral
$$
\int_{-t/2}^{t/2} \left( \frac{(z-x)}{\sqrt{(z - x)^2}} - \frac{(z-x)}{\sqrt{a^2 + (z - x)^2}} \right) dz\ ,
$$
$\ \sqrt{(z - x)^2}=|z-x|=x-z\ $, because $\ x>\frac{t}{2}\ $, and therefore $\ x>z\ $ for all $\ z\in\left[-\frac{t}{2}, \frac{t}{2}\right]\ $.
However, while the integral is correct as written, there are two other factors that determine the sign of its integrand.
The first is the physics of the situation.  You're  computing the component of the gravitational force acting on the particle in the direction of the positive $z$ axis, which must be negative, because the actual direction of the force is from $\ (0,0,x)\ $ towards the origin. Thus, in the integral
$$
\int_{-t/2}^{t/2}\int_{0}^{a} \frac{\rho (z - x)}{( \rho^2 + (z - x)^2 )^{3/2}} d\rho dz\ ,
$$
the integrand is properly negative over the whole range of integration.
The second factor is the choice of substitution for evaluating the integral.  When you make the substitution $\ x=f(u)\ $ to evaluate the integral
$$
\int_{0}^{a} \frac{\rho (z - x)}{( \rho^2 + (z - x)^2 )^{3/2}} d\rho\ ,
$$
you're effectively appealing to the theorem that
$$
\int_{u_1}^{u_2}g(f(u))f'(u)du=\int_{f\left(u_1\right)}^{f\left(u_2\right)}g(\rho)d\rho\ .
$$
Here, $\ f\ $ can be any absolutely continuous function which makes the valuation of the integral on the left easier to see by inspection than that of the one on the right.  In your substitution, where $\ g(\rho)=\frac{\rho (z - x)}{( \rho^2 + (z - x)^2 )^{3/2}}\ $, you're effectively choosing $\ \rho=f(u)\ $ such that $\ u^2=\rho^2+(z-x)^2\ $.  Since $\ \rho\ge0\ $ over the whole range of the relevant integral, this means that $\ f(u)=\sqrt{u^2-(z-x)^2}\ $, $\ f'(u)=\frac{u}{\sqrt{u^2-(z-x)^2}}=\frac{u}{\rho}\ $, and the integral on the left becomes
$$
\int_{u_1}^{u_2}\frac{(z-x)u}{|u|^3}du\ .
$$
Here, there are two possible values you can choose for $\ u_1\ $ and $\ u_2\ $—namely $\ z-x\ $ or $\   x-z\ $ for $\ u_1\ $, and $\ \sqrt{a^2+(z-x)^2}\ $ or $\ {-}\sqrt{a^2+(z-x)^2}\ $ for $\ u_2\ $.  While it's simpler to choose the two positive values $\  x-z\ $ and $\ \sqrt{a^2+(z-x)^2}\ $ for  the limits, it actually doesn't matter which values you choose, as long as you're careful with the form of the integrand when $\ u<0\ $:
$$
\frac{u}{|u|^3}=\cases{\frac{1}{u^2}&if $\ u\ge0\ $,\\
{-}\frac{1}{u^2}&if $\ u<0\ $,}
$$
and use the identity $\ \int_\limits{a}^b=-\int_\limits{b}^a\ $ whenever $\ a>b\ $. The integral will always evaluate to the same quantity, no matter which combination of choices you make for the limits
