Possible next step but not yet the final answer:
Working out first few terms in Mathematica and searching the OEIS, I found that my polynomials can be expressed in terms of "B-analogues of Stirling numbers of the second kind" $T(n,k)$, see sequence A039755 for definition and references.
Using the generating function from OEIS, I get
$p_k(y) = \sum_{m=0}^{k} T(k, m) (-1)^m y^{2 m} = \frac{\partial^k}{\partial z^k} g(z,y)$
where $g(z) = e^{z-(1/2)y^2 e^{2 z}} e^{y^2/2}$.
Thus $F(y) = \sum_{l=0}^{\infty} \frac{1}{l!} (i \beta)^l y e^{-y^2/2}\frac{\partial^{2 l}}{\partial z^{2 l}} g(z,y) $ .
Unfortunatelly, I have the square in the exponent and $l!$ does not match the dergee $2l$ in the expansion, and I don't see how to proceed.
UPDATE-1
Using "Dobinski-type identity" for the "row generating polynomials" (see OEIS entry), which in my notation is $p_k(y)=e^{y^2/2} \sum_{q=0}^{\infty} (q!)^{-1} (-2)^{-q} (2 q+1)^k y^{2 q}$, results in a series repsentation:
$F(y,\beta)=\sum_{q=0}^{\infty} (q!)^{-1} (-2)^{-q} e^{i \beta (2 q+1)^2} y^{2 q} $
This gives a useful asymptotic expansion of $F$ in terms of $\beta$ which, however does not capture correctly the asymptotics at large $y$.
Will try to manupulate this series, any advice is appreciated.