The case $p=2$ is treated briefly in the final section of the cited
Schmuland paper, which gives a picture of the distribution. The observations
in the paper's first few sections adapt to arbitrary $p>1$: the distribution
is not expected to have a simple formula, but the moment-generating function
has a product formula (exhibited below),
which we can use to compute each power moment ${\mathbb E}(X_p^{2k})$
of this random variable $X_p = \sum_n \epsilon_n/n^p$
as a polynomial in $\zeta(2p), \zeta(4p), \zeta(6p), \ldots, \zeta(2kp)$.
[The odd-order moments vanish by symmetry.]
This is because the distribution of $X_p$ is the convolution of
an infinite series of distributions, the $n$-th of which is supported on
$\pm 1 /n^p$ each with probability $1/2$; therefore
$$
{\mathbb E}(\exp(tX_p))
= \prod_{n=1}^\infty \left(\frac12 e^{t/n^p} + \frac12 e^{-t/n^p} \right)
= \prod_{n=1}^\infty \cosh(t/n^p).
$$
We can recover the power moments from the expansion of
${\mathbb E}(\exp(tX_p))$ in a Taylor series about $t=0$,
writing
$$
\begin{eqnarray}
\log {\mathbb E}(\exp(tX_p))
&=& \sum_{n=1}^\infty \log(\cosh(t/n^p)) \cr
&=& \sum_{n=1}^\infty
\frac12 \left(\frac{t}{n^p_{\phantom1}}\right)^2
- \frac1{12} \left(\frac{t}{n^p_{\phantom1}}\right)^4
+ \frac1{45} \left(\frac{t}{n^p_{\phantom1}}\right)^6
- \frac{17}{2520} \left(\frac{t}{n^p_{\phantom1}}\right)^8
+ - \cdots \cr
&=&
\frac{\zeta(2p)}{2} t^2
- \frac{\zeta(4p)}{12} t^4
+ \frac{\zeta(6p)}{45} t^6
- \frac{17\zeta(6p)}{2520} t^8
+ - \cdots
\end{eqnarray}
$$
and exponentiating.
While the random variables $X_p$ may be no more than objects of curiosity,
similar constructions arise naturally in analytic number theory, such as
the value at a fixed $s>1$ of the Dirichlet $L$-series associated to a random
real character $\chi$. In this example, the terms $\chi(n)/n^s$ in the
Dirichlet series are correlated, but the Euler product
$L(s,\chi) = \prod_l (1 - \chi(l)/l^s)^{-1}$ has independent factors, so
$\log L(s,\chi)$ is an infinite convolution of the same kind. Likewise
if the $\epsilon_n$ were random complex numbers of unit length:
the analogue could be $\log L(s,\chi)$ for a random Dirichlet charater $\chi$
that need not be real, or $\log \zeta(\sigma + it)$ for fixed $\sigma>1$
and random real $t$. For these sums of complex-valued random variables,
the factors in the moment generating function get more complicated than
hyperbolic cosines, but are still tractable.
