Heavy-tail Lambert W x F distributions seem to be what you are looking for (disclaimer: I am the author). They arise from a parametric, non-linear transformation of a random variable (RV) $X \sim F$, to a heavy-tailed version $Y \sim \text{Lambert W} \times F$. For $F$ being Gaussian they reduce to Tukey's $h$ distribution.
They have one parameter $\delta \geq 0$ that regulates the degree of tail heaviness (you can also choose different left and right heavy tails). In its most basic form it transforms a standard Normal $U \sim \mathcal{N}(0,1)$ to a Lambert W $\times$ Gaussian $Z$ by
$$
Z = U \exp\left(\frac{\delta}{2} U^2\right)
$$
If $\delta > 0$ $Z$ has heavier tails than $U$; for $\delta = 0$ $Z$ falls back to your initial RV $U$.
If you don't want to use the Gaussian distribution as your baseline, you can just create other Lambert W versions of your favorite distributions, e.g. t, uniform, gamma, exponential, beta, ... However, for the purpose of heavy-tails the Gaussian seems to be the obvious baseline reference.
Since this heavy-tail generation is based on transformations of RVs/data rather than a manipulation of cdfs or pdfs, it is very convenient to implement: just use your code to simulate any standard RV, add one line that transforms it, and then you have a random sample from the heavy-tail version of your initial RV.
In R this becomes (using the LambertW package)
library(LambertW)
set.seed(1)
zz = rLambertW(n=1000, distname = "normal", beta = c(0,1), delta = 0.5)
normfit(zz)
You can also fit the best model to the data using a maximum likelihood estimator (MLE)
model = MLE_LambertW(zz, distname = "normal", type = "h")
summary(model)
plot(model)