The distribution of the product of a multivariate normal and a lognormal distribution

I will start looking at some simpler case. Let $X, Z$ be independent standard normal random variables, and let $Y=e^Z$. Then $Y$ is lognormal, and we are interested in the distribution of the product $XY= Y mid X mid S$ where $S=text{sign}(X)$ is independent of $X$. I will study the distribution of $mid XYmid = mid X mid Y$. Then finally multiplying by $S$ simply represents a random reflection in the vertical axis. Using the Mellin transform simplifies working with the product of independent (and positive) random variables. The Mellin transform of (almost surely) positive $X$ is $$ DeclareMathOperator{E}{mathbb{E}} M_X(s) = E X^s = E e^{log(X^s)}=E e^{s log X} =K_{log X}(s) $$ where $K$ represents the moment generating function (mgf). So the Mellin transform of $X$ is simply the mgf of $log X$.

For our random variables we find (I used maple here) $$ M_{mid X mid}(s) = frac{2^{s/2}Gamma((s+1)/2)}{sqrtpi} \ M_Y(s) = e^{s^2/2} $$ Then for the Mellin transform of the product: $$ M_{mid Xmid Y}(s) =E (mid X mid Y)^s=E (mid Xmid^s Y^s)=\ E mid X mid^s E Y^s qquad text{(by independence)} \ = frac{e^{(slog 2 + s^2)/2}Gamma((s+1)/2)}{sqrtpi} $$ Then we need the inverse transform of this, or we can use it as a basis for approximation, maybe saddlepoint approximation.

For the multivariate case, if the multinormal vector has as components two standard normal variables, the above will give both marginals. That give the solution in the case of independence. I will try to come back for the dependent case.