Derivation of conformal OPE in D>2 Operator product expansion says that, the product of two primary fields(of same dimension in this case) can be expanded as sum of primaries and their descendants
$$\phi_1(x)\phi_2(0) = {\Large \Sigma_\mathcal{O}}\lambda_\mathcal{O}C_\mathcal{O}(x,\partial_y)\mathcal{O}(y)|_{y=0}
$$
where the summation $\Sigma_\mathcal{O}$ is over primaries. Descents appear when acted upon by the derivatives in $C_\mathcal{O}(x,\partial_y)$. Considering the three point function and since we know two point functions are diagonal we get,
\begin{equation}\langle\phi_1(x)\phi_2(0)\Phi(z)\rangle = \lambda_\Phi C_\Phi(x,\partial_y)\langle\Phi(y)|_{y=0}\Phi(z)\rangle \hspace{0.2cm} (1)
\end{equation}
Now using known forms of two and three point functions below 
$$\langle\phi_1(x)\phi_2(x_2)\Phi(x_3)\rangle = \frac{\lambda_\Phi }{|x_{12}|^{\Delta_1+\Delta_2-\Delta_3}|x_{23}|^{\Delta_2+\Delta_3-\Delta_1}|x_{13}|^{\Delta_1+\Delta_3-\Delta_2}}
$$
$$ \langle \Phi(y)\Phi(z)\rangle = \frac{1}{|y-z|^{2\Delta_\Phi}}
$$
one is supposed to fix the the constants $\alpha, \beta$ in $C_\Phi(x,\partial_y)$ by assuming a form
$$
C_\Phi(x,\partial_y) = \frac{1}{|x|^{2\Delta -\Delta_\Phi}}\Big[1+ \frac{1}{2}x^\mu\partial_\mu + \alpha x^\mu x^\nu \partial_\mu\partial_\nu + \beta x^2\partial^2 + ...\Big]
$$ 
Here, the dimensions of $\phi_1$ and $\phi_2$ are each $\Delta$ and that of $\Phi$ is $\Delta_\Phi$.
Now I can see using the three point function with insertions at $x,0,z$ on the LHS of (1) is 
$$\langle\phi_1(x)\phi_2(0)\Phi(z)\rangle = \frac{\lambda_\Phi }{|x|^{2\Delta -\Delta_\phi}|z|^{\Delta_\Phi}|z-x|^{\Delta_\Phi}}
$$
the leading term when expanding about $x$, $\frac{\lambda_\Phi }{|x|^{2\Delta -\Delta_\phi}|z|^{2\Delta_{\Phi}}} $ matches with the RHS of eq. (1) but can't figure out how to find the the coefficients of higher order terms. I find myself trying to evaluate a binomial expansion of $|z-x|^{\Delta_\Phi}$ where now the points are in $D$ dimensional space and that I am not able to do. I am only trying to get till two orders. Any help is appreciated.
 A: This is how you evaluate the series for
$$
\frac{1}{|z-x|^{\Delta_\phi}}
$$
for small $x$. First, you pull out $|z|$,
$$
\frac{1}{|z|^{\Delta_\phi}|e-\xi|^{\Delta_\phi}},
$$
where $e=\frac{z}{|z|}$, write $\xi=\frac{x}{|z|}$ and now work with
$$
\frac{1}{|e-\xi|^{\Delta_\phi}}=\left[\frac{1}{(e-\xi)^2}\right]^{\Delta_\phi/2},
$$
Using $(e-\xi)^2=1-2(e\cdot\xi)+\xi^2$ and replacing $\xi\to\epsilon\xi$ to track the order you get just
$$
\left[\frac{1}{1-2\epsilon(e\cdot\xi)+\epsilon^2\xi^2}\right]^{\Delta_\phi/2},
$$
which is now just a usual function of scalar argument $\epsilon$, which you can expand out in powers of $\epsilon$ by hand or using Mathematica. To get the general answer, set $\epsilon=t/|\xi|$ and get
$$
\left[\frac{1}{1-2t\frac{(e\cdot\xi)}{|\xi|}+t^2}\right]^{\Delta_\phi/2}=\sum_{j=0}^\infty C^{(\Delta_\phi/2)}_j\left(\frac{e\cdot\xi}{|\xi|}\right)t^j=\sum_{j=0}^\infty C^{(\Delta_\phi/2)}_j\left(\frac{e\cdot\xi}{|\xi|}\right)|\xi|^j\epsilon^j,
$$
by definition of the Gegenbauer polynomials, as noted by Abdelmalek Abdesselam.
