  We study functions related to the experimentally observed Havriliak-Negami
dielectric relaxation pattern in the frequency domain
$\sim[1+(i\omega\tau_{0})^{\alpha}]^{-\beta}$ with $\tau_{0}$ being some
characteristic time. For $\alpha = l/k< 1$ ($l$ and $k$ positive integers) and
$\beta > 0$ we furnish exact and explicit expressions for response and
relaxation functions in the time domain and suitable probability densities in
their "dual" domain. All these functions are expressed as finite sums of
generalized hypergeometric functions, convenient to handle analytically and
numerically. Introducing a reparameterization $\beta = (2-q)/(q-1)$ and
$\tau_{0} = (q-1)^{1/\alpha}$ $(1 < q < 2)$ we show that for $0 < \alpha < 1$
the response functions $f_{\alpha, \beta}(t/\tau_{0})$ go to the one-sided
L\'{e}vy stable distributions when $q$ tends to one. Moreover, applying the
self-similarity property of the probability densities $g_{\alpha, \beta}(u)$,
we introduce two-variable densities and show that they satisfy the integral
form of the evolution equation.
