I cite the derivation on Wikipedia, here. The focus of that section of the article is to show that the $\chi_p^2$ statistic is asymptotically equivalent to the $\chi^2$ distribution.
Let $n$ be the number of observations, $m$ the number of cells and $p_i$ the probability of an observation (according to the null hypothesis of the supposed underlying distribution) falling in the $i$th cell, $1\le i\le m$. Let $\{k_i\}$ denote the configuration where $k_i$ occurences are observed in the $i$th cell. Let $\chi^2_p(\{k_i\},\{p_i\})$ denote the test statistic, and let $\chi^2_p(\{p_i\})$ denote its distribution.
For any arbitrary real $T$, assuming that observations are multinomially distributed by the supposed underlying distribution: $$P(\chi^2_p(\{p_i\})\gt T)=\sum_{\{k_i\}:\chi^2_p(\{k_i\},\{p_i\})\gt T}\frac{n!}{k_1!\cdots k_m!}\prod_{i=1}^mp_i^{k_i}$$ As $n\to\infty$, one may use Stirling's approximation for the factorial: $$P(\chi^2_p(\{p_i\})\gt T)\sim\sqrt{\frac{2\pi n}{\prod_{i=1}^m2\pi k_i}}\sum_{\{k_i\}:\chi^2_p(\{k_i\},\{p_i\})\gt T}\prod_{i=1}^m\left(\frac{n p_i}{k_i}\right)^{k_i}$$
With the substitution $x_i=\frac{k_i-np_i}{\sqrt{n}},\,1\le i\le m-1$ and noting that $k_m=np_m-\sqrt{n}\sum_{i=1}^{m-1}x_i$, one finds: $$\begin{align}P(\chi^2_p(\{p_i\})\gt T)&\sim\sqrt{\frac{2\pi n}{\prod_{i=1}^m2\pi k_i}}\sum_{\{x_i\}:\chi^2_p(\{x_i\sqrt{n}+np_i\},\{p_i\})\gt T}\\&\left\{\left(1-\frac{\sum_{i=1}^{m-1}x_i}{p_m\sqrt{n}}\right)^{-(np_m-\sqrt{n}\sum_{i=1}^mx_i)}\prod_{i=1}^{m-1}\left(1+\frac{x_i}{p_i\sqrt{n}}\right)^{-(np_i+x_i\sqrt{n})}\right\}\end{align}$$
Ugly as this is so far, I have no problems. However, as $n\to\infty$, Wikipedia pass from the above sum to an $(m-1)$ dimensional integral over the $x_i$, with everything unchanged save for an extra (mystifying) factor of: $$\prod_{i=1}^{m-1}\left(\sqrt{n}\,\mathrm{d}x_i\right)$$
$$\begin{align}P(\chi^2_p(\{p_i\})\gt T)\sim\sqrt{\frac{2\pi n}{\prod_{i=1}^m2\pi k_i}}\color{red}{\int_{\chi^2_p(\{x_i\sqrt{n}+np_i\},\{p_i\})\gt T}}\hspace{100pt}\\\left\{\left(1-\frac{\sum_{i=1}^{m-1}x_i}{p_m\sqrt{n}}\right)^{-(np_m-\sqrt{n}\sum_{i=1}^mx_i)}\prod_{i=1}^{m-1}\left(1+\frac{x_i}{p_i\sqrt{n}}\right)^{-(np_i+x_i\sqrt{n})}\right\}\color{red}{\left\{\prod_{i=1}^{m-1}\sqrt{n}\,\mathrm{d}x_i\right\}}\end{align}$$
This step I don't understand. That extra factor suggests that $\mathrm{d}x_i\sim\frac{1}{\sqrt{n}}$, but how can this be formalised? Or is it that, under the null hypothesis, $k_i-np_i\approx0$ as $n\to\infty$? In passing from a sum over discrete variables to a continuous integral, work must be done and the article's authors have not provided it. If anyone could help me see how the above sum transforms to an $m-1$ dimensional Riemann sum, with limiting difference $\sqrt{n}\mathrm{d}x_i$, that would be greatly appreciated. I understand that we can change our perspective and view the $x_i$ as a free ranging variable. I guess I just don't get the meaning of the symbol $\mathrm{d}x_i$ in this context.
Define $Y_{n,i}=(k_i-np_i)/\sqrt{np_i}$ and let $Y_n:=[Y_{n,1} \quad\cdots\quad Y_{n,m}]^{\top}$ such that $\chi_p^2(\{p_i\})=Y_n^{\top}Y_n$. By the CLT (assuming i.i.d. observations), $$ Y_{n}\xrightarrow{d}Y\sim N(0,\Sigma), $$ where $\Sigma=I_m-aa^{\top}$ with $a=[\sqrt{p_1} \quad\cdots\quad \sqrt{p_m}]^{\top}$. Let $Q\Lambda Q^{\top}$ be the eigendecomposition of $\Sigma$, i.e., $\Sigma=Q\Lambda Q^{\top}$ and $QQ^{\top}=I_m$. Then $$ Y^{\top}Y= Y^{\top}QQ^{\top}Y, $$ with $Q^{\top}Y\sim N(0,\Lambda)$. It is not hard to check that $\Lambda=\operatorname{diag}(0,1,\ldots,1)$. Thus, $Y^{\top}Y$ is the sum of squares of $m-1$ independent $N(0,1)$ random variables. Finally, using the continuous mapping theorem, $$ \chi_p^2(\{p_i\})\xrightarrow{d}Y^{\top}Y\sim \chi_{m-1}^2. $$