There have been many studies on “ \({\mathbb {Z}}_2\) ” harmonic functions, differential forms, or spinors recently, for example Donaldson (Q J Math 72:1–2, 2021), He (arXiv preprint, 2022), Taubes and Wu (Proceedings of Gökova geometry-topology conference 2018–2019, 2021; J Differ Geom 128(1), 379–462, 2024). It is then natural to ask what the most simple case is like: \({\mathbb {Z}}_2\) harmonic functions on \({\mathbb {R}}^2\) with point singularities? This paper focus on this issue. Unlike the more ambitious ongoing program of studying \({\mathbb {Z}}_2\) harmonic functions f such that \(|df| = 0 \) on the singular locus K (see, for example, Donaldson (Q J Math 72:1–2, 2021) and He (arXiv preprint, 2022), this paper only requires |df| to have locally finite \(L^2\) norms near K, which means, |df| is allowed to be unbounded near K. To emphasize the difference and to reduce the confusion, we call them singular \({\mathbb {Z}}_2\) harmonic functions. This paper studies singular \({\mathbb {Z}}_2\) harmonic functions and finds its unexpected relationship with the polynomial Pell equation.