Assume that $\sqrt{n} (\boldsymbol{Z} - \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}) - 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})$ ?
We know that
$$
f_{h(\boldsymbol{Z})}(\boldsymbol{z}) = f_{\boldsymbol{Z}}(h^{-1}(\boldsymbol{z})) \, \left|\frac{d}{d \boldsymbol{z}} h^{-1}(\boldsymbol{z})\right|,
$$
where the last term denotes the Jacobian of the transformation.
If we let $\boldsymbol{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})}(\boldsymbol{z})}{f_{\boldsymbol{W}}(\boldsymbol{z})} = 1 + \frac{\text{error 1}}{\sqrt{n}} + \frac{\text{error 2}}{n} + ~..., \quad \text{as } n\to \infty,
$$
with appropriate restrictions on $h$ of course ?
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.