Let $p \ge 2$ be a positive integer, and let $Q \in \mathcal P(\mathbb Z_p)$ be a probability distribution on $\mathbb Z_p$.

>**Question.** What are necessary and sufficient conditions on $Q$ to ensure that the later admits a square-root w.r.t convolution, i.e such that there exists $D \in \mathcal P(\mathbb Z_p)$ verifying $D \star D = Q$ ?

I'm particularly, interested in the case where $Q$ is Zipf, i.e $Q(k) \propto (k+1)^{-\beta}$ for some $\beta \gt 1$.

An illustrative example
---
Consider the case where $p = 2$. Let $a := Q(0)$, $b:=Q(1)$, $x:=D(0)$, and $y := D(1)$. We are interested in the feasibility of the following system.
\begin{align}
a &= x^2 + y^2,\\
b &= xy + yx = 2xy,\\
1 &= x + y,\\
x &\ge 0,\\
y &\ge 0.
\end{align}
Substituting $y = x - 1$ gives $2x(1-x) = b$, i.e $2x^2 - 2x + b = 0$, which evaluates to
\begin{eqnarray}
x_\pm = \frac{2 \pm \sqrt{4 - 8b}}{4} = \frac{1 \pm \sqrt{1 - 2b}}{2}.
\end{eqnarray}
For this to be real, we require
$$
b \le 1/2.
$$
With this condition, note that $x_+ + x_- = 1$ and $x_+, x_- \ge 0$. Take $x=x_+$ and $y=x_-$.  Now, the first equation (the circle) reduces to the requirement

$$
a = x_+^2 + x_-^2 = 2\left(\frac{1}{4} + \frac{1-2b}{4}\right) = 1-b,
$$
i.e $a+b=1$, which is satisfied since $Q$ is a distribution. We conclude that in the case $p=2$, the answer to our question is affirmative if $b \le 1/2$. It's easy to check that this condition is also necessary (i.e $D$ doesn't exist when $b \gt 1/2$).


Edit: Exploring an idea based on Fourier analysis
----

This is based on some comments by user Paul Garrett.

----

So, $\mathbb Z_p$ is a compact (in fact finite!) group, and so we must have Fourier transform $\hat Q:\mathbb Z_p \to \mathbb Z_p$, given by
$$
\hat Q(k) = \sum_{x \in \mathbb Z_p} Q(x) e^{-2i\pi k / p}.
$$
The convolution theorem then gives
$$
\hat D^2(k) = \hat Q(k),\text{ for all }k \in \mathbb Z_p.
$$

For $p=2$, one has 
$$
\begin{split}
\hat Q(k) &= \sum_{x \in \mathbb Z_2} Q(x) e^{-i\pi kx} = Q(0) + Q(1) e^{-i\pi k} = a + be^{-i\pi k}.
\end{split}
$$
Thus, $\hat Q(0) = Q(0) + Q(1) = 1$ and $\hat Q(1) = Q(0) + Q(1) e^{-i\pi} = Q(0) - Q(1) = 1-b - b = 1 - 2b$. The convolution theorem above then translates to
$$
\hat D(0)^2 = 1,\quad \hat D(1)^2 = 1-2b.
$$
This is only solvable if $b \le 1/2$ (aha !), in which case we must have
$$
\hat D(0) = \pm 1,\quad \hat D(1) = \pm \sqrt{1-2b}.
$$

Inverting this hopefully recovers the result we previously obtained.