The Ising model is widely regarded as the most studied model of spin-systems in statistical physics. The focus of this paper is its dynamic (stochastic) version, the Glauber dynamics, introduced in 1963 and by now the most popular means of sampling the Ising measure. Intensive study throughout the last three decades has yielded a rigorous understanding of the spectral-gap of the dynamics on $\Z^2$ everywhere except at criticality. While the critical behavior of the Ising model has long been the focus for physicists, mathematicians have only recently developed an understanding of its critical geometry with the advent of SLE, CLE and new tools to study conformally invariant systems. A rich interplay exists between the static and dynamic models. At the static phase-transition for Ising, the dynamics is conjectured to undergo a critical slowdown: At high temperature the inverse-gap is O(1), at the critical $\beta_c$ it is polynomial in the side-length and at low temperature it is exponential in it. A seminal series of papers verified this on $\Z^2$ except at $\beta=\beta_c$ where the behavior remained a challenging open problem. Here we establish the first rigorous polynomial upper bound for the critical mixing, thus confirming the critical slowdown for the Ising model in $\Z^2$. Namely, we show that on a finite box with arbitrary (e.g. fixed, free, periodic) boundary conditions, the inverse-gap at $\beta=\beta_c$ is polynomial in the side-length. The proof harnesses recent understanding of the scaling limit of critical Fortuin-Kasteleyn representation of the Ising model together with classical tools from the analysis of Markov chains.