Assume that $\sqrt{n} (\boldsymbol{Z}_n - \boldsymbol{\mu}) \stackrel{\mathcal{D}}{\longrightarrow} \mathcal{N}(\boldsymbol{0},\Sigma)$, as $n\to \infty$, for some $\boldsymbol{\mu}\in \mathbb{R}^d$ and $\Sigma$ a symmetric positive definite matrix (here, $\stackrel{\mathcal{D}}{\longrightarrow}$ denotes the convergence in distribution). The delta method says that for a $\mathcal{C}^1$ function $h: \mathbb{R}^d \rightarrow \mathbb{R}$, we have
$$
\sqrt{n} (h(\boldsymbol{Z}_n) - h(\boldsymbol{\mu})) \stackrel{\mathcal{D}}{\longrightarrow} \mathcal{N}(\boldsymbol{0}, \nabla h (\boldsymbol{\mu})^{\top} \Sigma \, \nabla h (\boldsymbol{\mu})), \quad \text{as } n\to \infty.
$$

Now, my question is: Are there more precise results in the literature involving asymptotic expansion(s) for the DENSITY function of $h(\boldsymbol{Z}_n)$ (with appropriate conditions on $h$)?
We know that
$$
f_{h(\boldsymbol{Z}_n)}(t) = f_{\boldsymbol{Z}_n}(h^{-1}(t)) \, \left|\frac{d}{d t} h^{-1}(t)\right|, \quad t\in \mathbb{R},
$$
where the last term denotes the Jacobian of the transformation.
If we let $W\sim \mathcal{N}(\boldsymbol{0}, \nabla h (\boldsymbol{\mu})^{\top} \Sigma \, \nabla h (\boldsymbol{\mu}))$, would it be possible to obtain a result of the form:
$$
\frac{f_{h(\boldsymbol{Z}_n)}(t)}{f_{W}(t)} = 1 + \frac{\text{error}_1(t)}{\sqrt{n}} + \frac{\text{error}_2(t)}{n} + ~..., \quad \text{as } n\to \infty,
$$
with appropriate restrictions on $h$ ?
I found these papers:
- https://projecteuclid.org/journals/annals-of-statistics/volume-13/issue-1/A-Converse-to-Scheffes-Theorem/10.1214/aos/1176346604.full
- https://www.projecteuclid.org/journals/annals-of-statistics/volume-14/issue-3/On-a-Converse-to-Scheffes-Theorem/10.1214/aos/1176350065.full

 but they don't quite answer my question.