(Cross-posted from mathoverflow: No answers yet; bounty there expires in less than 24 hours)
I'm looking for a function $f_p(x)$ with real parameter $p>0$ satisfying $$ \int_0^\infty f_p(x)x^{s-1}dx=e^{-p\psi(s)} $$ where $\psi(s)$ is the usual digamma function. The inverse Mellin formula is $$ f_p(x)=\frac{1}{2\pi i}\int_{c-i\infty}^{c+i\infty}x^{-s}e^{-p\psi(s)}ds $$ Normally we would try the residue theorem; this is possible but there are essential singularities and it is not obvious (to me) how to find the residues. The preprint www.arxiv.org/abs/1312.1604 gives asymptotic expansions for the exponential of the digamma function that may be useful.
Edit: I think $f_p(x)$ is supported on the unit interval; note that the first integral is positive and does not grow exponentially with $s$. Expanding in the relevant Chebyshev polynomials, ie
$ f_p(x)=\sum_{n=0}^\infty F_{n,p}\frac{T_{2n}(\sqrt{x})}{\sqrt{x(1-x)}} $
and noting that the digamma at integer values is just the harmonic numbers, $\psi(n)=H_{n-1}-\gamma$, we have
$ e^{-p(H_n-\gamma)}=\sum_{m=0}^\infty F_{m,p}\int\frac{T_{2m}(\sqrt{x})}{\sqrt{x(1-x)}}x^ndx=\sum_{m=0}^n\frac{\pi F_{m,p}}{2^{2n-1}} \frac{n (2n-1)!}{(n+m)!(n-m)!} $
where $n(2n-1)!$ is taken to be unity when $n=0$. The integral was done on mathematica by choosing values of $n$ and observing a pattern (would be good to have a reference/proof). This gives an explicit set of equations to find the $F_{m,p}$. These can then be used to reconstitute $f_p(x)$ numerically, using 200 terms in the expansion. The plot shows $f_p(x)$ for $p\in\{1/3,2/3,\ldots,4\}$ from top to bottom near $x=1$, ie behaviour roughly as $c(1-x)^{p-1}$. Convergence starts to break down at both ends of this range of $p$ (the observed oscillations).
But the question remains as to whether $f_p(x)$ can be expressed explicitly in terms of known functions.