EDIT. Your function $\mathfrak{h}$ is not modular of any level. Indeed $\mathfrak{h}$ has Fourier coefficients in $\mathbf{Q}$. A standard argument shows that any modular function $h$ with Fourier coefficients in some number field $K$ must have bounded denominators: take $n$ large enough such that $h \Delta^n$ is a cusp form. Shimura has proved that the space of cusp forms is generated by modular forms with integral coefficients. This implies that $h \Delta^n$ has bounded denominators. Since $\Delta$ is invertible in $\mathbf{Z}((q))$, the function $h$ also has bounded denominators. So if $\mathfrak{h}$ were modular then it would have bounded denominators. But $\mathfrak{h}^2$ has integral coefficients, so by Gauss's Lemma $\mathfrak{h}$ would also have integral coefficients. But the Fourier expansion of $\mathfrak{h}$ begins with $q^{-1/96} (1+\frac12 q^{3/2}+\cdots)$, a contradiction.

More generally, it is not always true that the $d$-th root of a modular function of level $N$ is modular of level $dN$.

One explanation is as follows: let $\Gamma=\Gamma(N)$ be the principal congruence subgroup of $\mathrm{SL}_2(\mathbf{Z})$ with $N \geq 3$. Then $Y(N)=\mathcal{H}/\Gamma(N)$ is a Riemann surface of some genus $g_N$ with $t_N$ points removed (the cusps). The group $\Gamma(N)$ is torsion-free, hence is isomorphic to the fundamental group of $Y(N)$. But it is known that this fundamental group is free of rank $2g_N+t_N-1$. In particular its rank is at least $2g_N$, and for any prime $p$, the maximal abelian $p$-torsion quotient of $\Gamma(N)$ is a $\mathbf{F}_p$-vector space of dimension at least $2g_N$, which grows roughly like $c \cdot N^3$. On the other hand, if $p$ divides $N$ then it is easy to see that $\Gamma(N)/\Gamma(pN)$ is isomorphic to $(\mathbf{Z}/p\mathbf{Z})^3$, which is much smaller. Now if $u$ is a modular unit of level $N$, then $u^{1/p}$ is a non-vanishing function on $\mathcal{H}$. In general the stabilizer of $u^{1/p}$ will be a normal subgroup of index $p$ in $\Gamma(N)$, but by the above analysis I see no reason this stabilizer will contain $\Gamma(pN)$.

You can construct infinitely many examples as follows (there may be simpler and more explicit ways to do that). Let $U$ be the group of all modular units of all levels, and let $S$ be the subgroup of $U$ generated by the Siegel units ($S$ is known as the Siegel group). Kubert has shown that $U/S$ is an abelian group of exponent $2$ and infinite rank (see The square root of the Siegel group, Introduction and Thm 4.19). Now take any modular unit $u \in U \backslash S$. If $v=u^{1/2}$ were modular then $u=v^2$ would belong to $U^2$, hence to $S$, contradiction. This shows that in a suitable sense, the square roots of most modular units are not modular of any level.

EDIT. Your function $\mathfrak{h}$ is not modular of any level. Indeed $\mathfrak{h}$ has Fourier coefficients in $\mathbf{Q}$. A standard argument shows that any modular unit $u$ with Fourier coefficients in some number field $K$ must have bounded denominators: take $n$ large enough such that $u \Delta^n$ is a cusp form. Shimura has proved that the space of cusp forms is generated by modular forms with integral coefficients. This implies that $u \Delta^n$ has bounded denominators. Since $\Delta$ is invertible in $\mathbf{Z}((q))$, the function $u$ also has bounded denominators. So if $\mathfrak{h}$ were modular then it would have bounded denominators. Since $\mathfrak{h}^2 = \mathfrak{f}$ has integral coefficients, Gauss's Lemma implies that $\mathfrak{h}$ also has integral coefficients. But the Fourier expansion of $\mathfrak{h}$ begins with $q^{-1/96} (1+\frac12 q^{3/2}+\cdots)$, a contradiction.

More generally, it is not always true that the $d$-th root of a modular function of level $N$ is modular of level $dN$.

One explanation is as follows: let $\Gamma=\Gamma(N)$ be the principal congruence subgroup of $\mathrm{SL}_2(\mathbf{Z})$ with $N \geq 3$. Then $Y(N)=\mathcal{H}/\Gamma(N)$ is a Riemann surface of some genus $g_N$ with $t_N$ points removed (the cusps). The group $\Gamma(N)$ is torsion-free and is isomorphic to the fundamental group of $Y(N)$. But it is known that this fundamental group is a free group of rank $2g_N+t_N-1 \geq 2g_N$. In particular for any prime $p$, the maximal abelian $p$-torsion quotient of $\Gamma(N)$ is a $\mathbf{F}_p$-vector space of dimension at least $2g_N$, which grows roughly like $c \cdot N^3$. On the other hand, if $p$ divides $N$ then it is easy to see that $\Gamma(N)/\Gamma(pN)$ is isomorphic to $(\mathbf{Z}/p\mathbf{Z})^3$, which is much smaller. Now if $u$ is a modular unit of level $N$, then $u^{1/p}$ is a non-vanishing function on $\mathcal{H}$. In general the stabilizer of $u^{1/p}$ will be a normal subgroup of index $p$ in $\Gamma(N)$, but by the above analysis I see no reason this stabilizer will contain $\Gamma(pN)$.

You can construct infinitely many examples as follows (there may be simpler and more explicit ways to do that). Let $U$ be the group of all modular units of all levels, and let $S$ be the subgroup of $U$ generated by the Siegel units ($S$ is known as the Siegel group). Kubert has shown that $U/S$ is an abelian group of exponent $2$ and infinite rank (see The square root of the Siegel group, Introduction and Thm 4.19). Now take any modular unit $u \in U \backslash S$. If $v=u^{1/2}$ were modular then $u=v^2$ would belong to $U^2$, hence to $S$, contradiction. This shows that in a suitable sense, the square roots of most modular units are not modular of any level.

I don't know the answer for your specific modular function, but in general it is not always true that the $d$-th root of a modular function of level $N$ is modular of level $dN$.

One explanation is as follows: let $\Gamma=\Gamma(N)$ be the principal congruence subgroup of $\mathrm{SL}_2(\mathbf{Z})$ with $N \geq 3$. Then $Y(N)=\mathcal{H}/\Gamma(N)$ is a Riemann surface of some genus $g_N$ with $t_N$ points removed (the cusps). The group $\Gamma(N)$ is torsion-free, hence is isomorphic to the fundamental group of $Y(N)$. But it is known that this fundamental group is free of rank $2g_N+t_N-1$. In particular its rank is at least $2g_N$, and for any prime $p$, the maximal abelian $p$-torsion quotient of $\Gamma(N)$ is a $\mathbf{F}_p$-vector space of dimension at least $2g_N$, which grows roughly like $c \cdot N^3$. On the other hand, if $p$ divides $N$ then it is easy to see that $\Gamma(N)/\Gamma(pN)$ is isomorphic to $(\mathbf{Z}/p\mathbf{Z})^3$, which is much smaller. Now if $u$ is a modular unit of level $N$, then $u^{1/p}$ is a non-vanishing function on $\mathcal{H}$. In general the stabilizer of $u^{1/p}$ will be a normal subgroup of index $p$ in $\Gamma(N)$, but by the above analysis there is no reason this stabilizer will contain $\Gamma(pN)$.

If you want an explicit example, you can proceed as follows (there may be simpler and more explicit ways to do that). Let $U$ be the group of all modular units of all levels, and let $S$ be the subgroup of $U$ generated by the Siegel units ($S$ is known as the Siegel group). Kubert has shown that $U/S$ is an abelian group of exponent $2$ and infinite rank (see The square root of the Siegel group, Introduction and Thm 4.19). Now take any modular unit $u \in U \backslash S$. If $v=u^{1/2}$ were modular then $u=v^2$ would belong to $U^2$, hence to $S$, contradiction.

EDIT. Your function $\mathfrak{h}$ is not modular of any level. Indeed $\mathfrak{h}$ has Fourier coefficients in $\mathbf{Q}$. A standard argument shows that any modular function $h$ with Fourier coefficients in some number field $K$ must have bounded denominators: take $n$ large enough such that $h \Delta^n$ is a cusp form. Shimura has proved that the space of cusp forms is generated by modular forms with integral coefficients. This implies that $h \Delta^n$ has bounded denominators. Since $\Delta$ is invertible in $\mathbf{Z}((q))$, the function $h$ also has bounded denominators. So if $\mathfrak{h}$ were modular then it would have bounded denominators. But $\mathfrak{h}^2$ has integral coefficients, so by Gauss's Lemma $\mathfrak{h}$ would also have integral coefficients. But the Fourier expansion of $\mathfrak{h}$ begins with $q^{-1/96} (1+\frac12 q^{3/2}+\cdots)$, a contradiction.

More generally, it is not always true that the $d$-th root of a modular function of level $N$ is modular of level $dN$.

One explanation is as follows: let $\Gamma=\Gamma(N)$ be the principal congruence subgroup of $\mathrm{SL}_2(\mathbf{Z})$ with $N \geq 3$. Then $Y(N)=\mathcal{H}/\Gamma(N)$ is a Riemann surface of some genus $g_N$ with $t_N$ points removed (the cusps). The group $\Gamma(N)$ is torsion-free, hence is isomorphic to the fundamental group of $Y(N)$. But it is known that this fundamental group is free of rank $2g_N+t_N-1$. In particular its rank is at least $2g_N$, and for any prime $p$, the maximal abelian $p$-torsion quotient of $\Gamma(N)$ is a $\mathbf{F}_p$-vector space of dimension at least $2g_N$, which grows roughly like $c \cdot N^3$. On the other hand, if $p$ divides $N$ then it is easy to see that $\Gamma(N)/\Gamma(pN)$ is isomorphic to $(\mathbf{Z}/p\mathbf{Z})^3$, which is much smaller. Now if $u$ is a modular unit of level $N$, then $u^{1/p}$ is a non-vanishing function on $\mathcal{H}$. In general the stabilizer of $u^{1/p}$ will be a normal subgroup of index $p$ in $\Gamma(N)$, but by the above analysis I see no reason this stabilizer will contain $\Gamma(pN)$.

You can construct infinitely many examples as follows (there may be simpler and more explicit ways to do that). Let $U$ be the group of all modular units of all levels, and let $S$ be the subgroup of $U$ generated by the Siegel units ($S$ is known as the Siegel group). Kubert has shown that $U/S$ is an abelian group of exponent $2$ and infinite rank (see The square root of the Siegel group, Introduction and Thm 4.19). Now take any modular unit $u \in U \backslash S$. If $v=u^{1/2}$ were modular then $u=v^2$ would belong to $U^2$, hence to $S$, contradiction. This shows that in a suitable sense, the square roots of most modular units are not modular of any level.