Given a non-zero polynomial $P(x)$, we study Fuchsian differential operators of the form $L=\partial_x^2-u(x)$ such that for all $\lambda\in\mathbb{C}$ the operator $L+\lambda P(x)$ is monodromy free. We prove that all such operators are obtained from populations of critical points of ${\widehat{\mathfrak{sl}}_2}$ master functions. Moreover, we show that the reproduction procedure of critical points corresponds to a Darboux transformation of operator $P^{-1}(x)L$. As a result, we obtain a classification of all operators $L$ with such properties in the case of $P(x)=x^k$.