I've not seen this notation before. However, it does seem to make sense. The $M$-downsampler is defined by the equation:
$$
y_D[n] = x[Mn]
$$
Its $z$ transform is defined by the equation:
$$
\begin{align}
Y_D(z) &= \sum_{n=-\infty}^\infty y_D[n]z^{-n} \\
&= \sum_{n=-\infty}^\infty x[Mn] z^{-n}
\end{align}
$$
Apply a change of variable, letting $n' = Mn$. The ranges of the summation are unaffected by the change of variable since they extend to infinity.
$$
Y_D(z) = \sum_{n'=-\infty}^\infty x[n']z^{-n'/M}
$$
This looks similar to the $z$ transform of $x[n]$ itself. Recall that it is defined as:
$$
X(z) = \sum_{n=-\infty}^\infty x[n]z^{-n}
$$
By inspection, we can therefore conclude the following relationship between the $z$ transforms of $x[n]$ and $y_D[n]$:
$$
Y_D(z) = X\left(z^{1/M}\right)
$$
Therefore, the $z$ transform of the downsampler output is closely related to the $z$ transform of the input signal, which is to be expected. In the frequency domain, this results in an $M$-fold stretching of the signal's frequency content.
But how do you go from the above equation to the one you referenced in the paper? It gives a definition of $Y_D(z)$ in terms of $z$ only, while the expression we derived is a function of $z^{1/M}$. So for a particular value of $z$ that you would like to evaluate $Y_D(z)$ at, you would first calculate $z^{1/M}$ (i.e. take the $M$-th root of $z$) and then substitute that into $X(z)$. However, all nonzero $z \in \mathbb{C}$ have $M$ distinct $M$-th roots:
$$
\left\{ r_p,\ r_p e^{\frac{j2\pi}{M}},\ r_p e^{\frac{j2\pi2}{M}},\ \ldots\ ,\ r_p e^{\frac{j2\pi(M-1)}{M}} \right\}
$$
$$
= \left\{ r_p,\ r_p W,\ r_p W^{2},\ \ldots\ ,\ r_p W^{M-1} \right\}
$$
where $W_k$ is the DFT kernel value $e^{j2\pi k / M}$ referenced in your question, and $r_p$ is what I define to be the principal $M$-th root of the complex value $z$:
$$
r_p = \sqrt[M]{|z|} e^{\frac{j\angle{z}}{M}}
$$
That is, $z$'s principal $M$-th root $r_p$ is obtained by converting $z$ to polar form, taking the $M$-th root of $z$'s magnitude (which is a real number), and dividing $z$'s angle by $M$. The resulting values express $r_p$ in polar form.
Why go to all of this trouble? Because, as I noted before, the mapping from $Y_D(z)$'s domain to the domain of $X(z^{1/M})$ is not one-to-one. I'll now begin some handwaving. For any particular value of $z$ that you would like to evaluate $Y_D(z)$ for, there are $M$ corresponding points in $X(z^{1/M})$ that you could map to. Therefore, each of those $M$ points in $X(z^{1/M})$ contribute to the corresponding value of $Y_D(z)$. You then end up with a sum like that shown in the paper:
$$
Y_D(z) = \frac{1}{M} \sum_{k=0}^{M-1} X(r_p(z)W^k)
$$
where $r_p(z)$ refers to the principal $M$-th root calculation I showed earlier. In reality, you could pick any of $z$'s $M$-th roots as the principal one; I picked this definition because it's the most straightforward. If you were to properly and rigorously derive this relationship, I believe the factor of $\frac{1}{M}$ comes in because of a derivative of $z^{1/M}$.
In mathematician-speak, I believe this would be referred to as a composition of functions; $Y_D(z) = f(g(z))$, where $f(z) = X(z)$ and $g(z) = z^{1/M}$. In order to unroll the function composition and write $Y_D(z)$ as a function of $z$ only, you would chop the domain of $Y_D(z)$ into chunks that are one-to-one, invert the function over those intervals, and then sum the results with appropriate scaling factors. I've used this technique before to calculate the probability distribution function of a function of a random variable given the original random variable's pdf (e.g. to derive the pdf of $\sqrt{X}$ given $X$'s pdf), but the name of the technique escapes me.
