While $1+1$ dimensional growth models in the Kardar-Parisi-Zhang universality class have witnessed an explosion of activity, higher dimensional models remain much less explored. The special case of $2+1$ dimensions is particularly interesting as it is, in physics parlance, neither ultraviolet nor infrared super-renormalizable. Canonical examples include the stochastic heat equation (SHE) with multiplicative noise and directed polymers. The models exhibit a weak to strong disorder transition as the inverse temperature, up to a logarithmic (in the system size) scaling, crosses a critical value. While the sub-critical picture has been established in detail, very recently [CSZ '23] constructed a scaling limit of the critical $2+1$ dimensional directed polymer partition function, termed as the critical $2d$ Stochastic Heat Flow (SHF), a random measure on $\mathbb{R}^2.$ The SHF is expected to exhibit a rich intermittent behavior and consequently a rapid growth of its moments. The $h^{th}$ moment was known to grow at least as $\exp(Ω(h^{2}))$ (a consequence of the Gaussian correlation inequality) and at most as $\exp(\exp (O(h^2)))$. The true growth rate, however, was predicted to be $\exp(\exp (Θ(h)))$ in the late nineties [R '99]. In this paper we prove a lower bound of the $h^{th}$ moment which matches the predicted value, thereby exponentially improving the previous lower bound. We also obtain rather sharp bounds on its upper tail. The key ingredient in the proof involves establishing a new connection of the SHF and moments thereof to the Gaussian Free Field (GFF) on related Feynman diagrams. This connection opens the door to the rich algebraic structure of the GFF to study the SHF. Along the way we also prove a new monotonicity property of the correlation kernel for the SHF as a consequence of the domain Markov property of the GFF.