I know that if we fix $\alpha_1,\dots,\alpha_d\ge1$ and $S\ge1$ integers, and we define the $d$-dimensional process $(Z(n))_{n\ge0}$ as
follows: let $V$ be a Dirichlet random variable of parameter $\big(\frac{\alpha_1}S,\dots,\frac{\alpha_d}S\big)$ and let $(X_m)_{m\ge1}$ be a
sequence of random variables such that, for all $k\ge1$, for all $1\le i \le d$, $\mathbb{P}(X_k=e_i|V)=V_i$. Let $Z_i(0)=\alpha_i\ge1$ for all $1\le i\le d$ and, for all $n\ge1$, $Z(n)=Z(0)+S\sum^n_{k=1}X_k$. Then, $(Z(n))_{n\ge0}$ is a Pólya urn of initial composition $(\alpha_1,\dots,\alpha_d)$ and replacement matrix $S\cdot$ Id.
To understand the difference between $Z(n)$ and $Z_i(n)$, we have, for the canonical base in $\mathbb{R}^d$: $Z(n)=\sum^d_{i=1}Z_i(n)e_i$, or in words, for a $d$-colour Pólya urn $U(n)$, $U_i(n)$ denotes the number of balls of colour $i$ at time $n$.
Using this information, and with the help of the multidimensional CLT I want to prove the following:
Theorem. Assume that $(U(n))_{n\ge0}$ is a Pólya urn of initial composition
$(\alpha_1,\dots,\alpha_d)$ and replacement matrix $R=S\cdot Id$ for some $\alpha_1,\dots,\alpha_d$, $S\ge1$. Then, in distribution as $n\uparrow\infty$,$$\frac{U(n)-nSV}{S\sqrt{n}}\xrightarrow{d}\mathcal{N}(0,\Sigma^2)$$ where $\Sigma^2=\sum^d_{i=1}V_ie_i^te_i-V^tV$.
To find this, we take advantage of the fact that, in distribution, $U(n)=Z(n)$ $\forall n$
Conditionally on $V$, we have: $\mathbb{E}[X|V]=e_i\cdot V_i \ (\forall1\le i\le d)=V$
Now, we want to apply the multidimensional CLT to $Z(n)$, since it is defined by a summation of a constant $(Z(0))$ and a sum of random vectors. The $X_i$ are i.i.d. (identically distributed is clear by definiton and for independence note that $\mathbb{E}[X_iX_j]=\mathbb{E}[X_i]\mathbb{E}[X_j]$ for any pair $1\le i,j\le n$, $i\not=j$ using Law of total expectation).
So can we apply the CLT to $Z(n)$ or do we better apply it to $Z(n)-Z(0)$? Either way, by the CLT, conditionally on $V$ we get $$\frac{SX_1+\dots+SX_n-n\mathbb{E}[SX|V]}{\sqrt{n}}=\frac{Z(n)-Z(0)-nSV}{\sqrt{n}}\xrightarrow{d}\mathcal{N}(0,\Sigma^2)$$ where $\Sigma^2=\text{Cov}(X|V)=\sum^d_{i=1}V_ie_i^te_i-V^tV$, where $X$ is a copy of $X_1$.
Does this makes sense? Where does the $S^{-1}$ factor in the Theorem comes from?
Reference: PÒLYA URNS AND OTHER REINFORCEMENT PROCESSES
Best Answer
Notation: Fix $S \geq 1$ and $\alpha \in \mathbb{R}^d_+$. Define the process $\{U_t\}_{t \geq 0} \subset \mathbb{R}^d_+$ by $$ U_{t+1} = U_{t} + S e_{X_{t+1}}, \quad \mbox{for}~t=0, 1, 2, \ldots. $$ Above, we have $X_{t} | V \stackrel{\rm iid}{\sim} \mathsf{Categorical}(V)$ for all $t \geq 1$ and $V \sim \mathsf{Dirichlet}(\alpha)$.
Applying the CLT: Define $Y_t = U_t - U_0$, $$ Y_t = S \sum_{\tau = 1}^t e_{X_{\tau}}, \quad \mbox{and} \quad \mathbb{E}[Y_t \mid V] = t S V. $$ For any $t \geq 1$, and integers $1 \leq i,j\leq d$, $\mathrm{Cov}(e_{X_t} \mid V)_{ij} = \delta_{ij} V_i- V_i V_j$. Consequently, $$ \Sigma \equiv \mathrm{Cov}(e_{X_t} \mid V) = \sum_{i=1}^d V_i (e_i \otimes e_i) - V \otimes V. $$ Now we can apply the CLT conditionally. Define $$ W_t = \frac{Y_t - \mathbb{E}[Y_t \mid V]}{S\sqrt{t}} \equiv \frac{U_t - tSV}{S\sqrt{t}}, \quad t \geq 1. $$ Since $X_t$ are conditionally iid, the CLT yields conditionally $W_t \to \mathsf{N}(0, \Sigma)$ in distribution as $t \to \infty.$ Formally, this means for any $u \in \mathbb{R}^d$, $$ \lim_{t\to\infty}\mathbb{P}(W_t \leq u \mid V) \to \mathbb{P}(W_V \leq u \mid V) \quad \mbox{where} \quad W_V \sim \mathsf{N}\big(0, \Sigma(V)\big). $$ (Above $w \leq u $ means $w_i \leq u_i$ for $i = 1, \ldots, d$.)
Note that this also implies (e.g., by dominated convergence) $$ \lim_{t\to\infty} \mathbb{P}(W_t \leq u) = \mathbb{E}_V[\mathbb{P}(W_V \leq u)]. $$