Is there a direct formula or a bound for the following Orlicz norm: $\|\xi \eta\|_{\psi_1}$, where $\xi, \eta \sim \mathcal{N}(0, 1)$ are standard normal random variables, $\psi_1(x) = e^x - 1$?
One way to bound this quantity would be $\|\xi \eta\|_{\psi_1} \leqslant \|\xi\|_{\psi_2} \|\eta\|_{\psi_2}$, and then compute the $\psi_2$-norm of a standard normal variable. Can this bound be improved? It gives the answer $8 / 3$.
A reminder: $$ \|\xi\|_{\psi_1} = \inf \{t > 0 | \mathbb{E}\exp\{|\xi| / t\} \leqslant 2\} $$
EDIT: It seems that the previous bound can indeed be improved. If one plots the approximate expectations, one can notice that the true value of the norm is below 2.
In order to calculate the sub-exponential norm, define $$ f(t) := \mathbb{E}[\exp(|\xi|/t)] $$ for $t > 0$. Then, we can simplify $f$ as follows by recalling the density of the Gaussian distribution: \begin{equation*} f(t) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty \exp(\frac{|z|}{t}) \exp(-\frac{z^2}{2}) dz \end{equation*} Now, I admit I used an integral calculator when I was first interested in this and plugged in the right side to get $$ f(t) = \exp(\frac{1}{2t^2}) \left( \text{erf}(\frac{1}{\sqrt{2}t}) + 1 \right) $$ where erf denotes the Gaussian error function. (If you are interested in the steps, I can try and do the calculation again by hand, let me know in the comments)
Since $f$ is decreasing on $\mathbb{R}_{>0}$, we can now simply set the right side equal to $2$ in order to get the sub-exponential norm, which I did using Geogebra to get the approximate value I mentioned above (1.3725).
Note that the same procedure can be done for the sub-gaussian norm, which gives you $\sqrt{8/3}$ for the standard normal distribution.