Let $f(y) = e^{y-e^y}$. I want to find $$ f_{(n)} = \underbrace{f * f * \dots * f}_{n} $$ where $*$ denotes convolution. I'm sure that Fourier or Laplace transforms are the key but I don't have a lot of experience with those.
So far I worked out $$ M_f(t) = \int_{\mathbb R}f(y)e^{ty}\,\text dy = \int_{\mathbb R}e^{-e^y}(e^y)^te^y\,\text dy = \int_0^\infty e^{-x}x^{t}\,\text dx = \Gamma(t+1). $$
(I'm coming at this from probability so I'm using moment generating function notation).
Does this mean that now I have $$ M_{f_{(n)}}(t) = \prod_{j=1}^n M_f(t) = \Gamma(t+1)^n? $$
And furthermore, how do I invert this to get $f_{(n)}$? I read through the wikipedia articles on inverse Fourier and Laplace transforms but they involve complex analysis which I don't know much of, so I'm not sure how to deal with the imaginary units, but also I don't think I need those for my particular real-valued case, do I? Thanks a lot for any help or references on the topic.
You are correct so far. Now you need to exploit the relation between the bilateral Laplace transform $\mathfrak{B}$ and the Mellin transform $\mathfrak{M}$, in particular you need a version of the Cahen-Mellin integral $$ e^{-y} = \frac{1}{2\pi i}\int_{c-i\infty}^{c+i\infty}\,\Gamma(s)\,y^{-s}\,ds\qquad (c>0) $$ where $\Gamma(s)$ is replaced by $\Gamma(s+1)^n$, i.e. you need to compute the residues of $\Gamma(s+1)^n y^{-s}$ at $s=-1,-2,-3,\ldots$. At the end of the process, by replacing $y$ with $e^x$ you will get $f_{(n)}(x)$.
Every negative integer is a pole of order $n$ for the function $\Gamma(s+1)^n$, so the computation of the needed residues is doable but becomes more and more involved as $n$ increases.
The $n=2$ case is yet pretty non-elementary
$$\begin{eqnarray*} \operatorname*{Res}_{s=-m}\Gamma(s+1)^2 y^{-s}&=&\lim_{s\to -m}\frac{d}{ds}(s+m)^2\Gamma(s+1)^2 y^{-s}\\&=&y^m\lim_{s\to -m}2(m+s)\Gamma(1+s)^2-(m+s)^2\Gamma(1+s)^2\log(y)+2(m+s)^2\Gamma(s+1)^2\psi(s+1)\\&=&\frac{y^m}{(m-1)!^2}\left(2H_{m-1}-2\gamma-\log y\right)\end{eqnarray*}$$ and it leads to $$f_{(2)}(x) = 2e^x \cdot K_0(2e^{x/2}) $$ where the Bessel $K_0$ function also appears in the normal product distribution.
The central limit theorem seems to be a simpler way for dealing with $f_{(n)}(x)$ as $n\to +\infty$, for sure: $f_{(n)}(x)$ converges to the distribution of a normal variable with mean value $-n\gamma$ and variance $\frac{\pi^2}{6} n$.