Using Quantum Monte Carlo simulations, we study a series of models of fermions coupled to quantum Ising spins on a square lattice with $N$ flavors of fermions per site for $N=1,2$ and $3$. The models have an extensive number of conserved quantities but are not integrable, and have rather rich phase diagrams consisting of several exotic phases and phase transitions that lie beyond Landau-Ginzburg paradigm. In particular, one of the prominent phase for $N>1$ corresponds to $2N$ gapless Dirac fermions coupled to an emergent $mathbb{Z}_2$ gauge field in its deconfined phase. However, unlike a conventional $mathbb{Z}_2$ gauge theory, we do not impose the `Gausss Law by hand and instead, it emerges due to spontaneous symmetry breaking. Correspondingly, unlike a conventional $mathbb{Z}_2$ gauge theory in two spatial dimensions, our models have a finite temperature phase transition associated with the melting of the order parameter that dynamically imposes the Gausss law constraint at zero temperature. By tuning a parameter, the deconfined phase undergoes a transition into a short range entangled phase, which corresponds to Neel/Superconductor for $N=2$ and a Valence Bond Solid for $N=3$. Furthermore, for $N=3$, the Valence Bond Solid further undergoes a transition to a Neel phase consistent with the deconfined quantum critical phenomenon studied earlier in the context of quantum magnets.