Let $\mathcal{F}$ be an $\mathbb{R}$-algebra. We write $|\mathcal{F}|$ for the set of surjective algebra homomorphisms $\mathcal{F}\rightarrow\mathbb{R}$, and we say $\mathcal{F}$ is geometric when $\mathcal{I}(\mathcal{F}):=\bigcap\limits_{p\in|\mathcal{F}|}\ker p=\{0\}$.
Let $f\in\mathbb{R}[x_1,\dots,x_n]$. When is $\mathbb{R}[x_1,\dots,x_n]/f\mathbb{R}[x_1,\dots,x_n]$ geometric?
Here's what I know so far: any algebra homomorphism $\mathbb{R}[x_1,\dots,x_n]/f\mathbb{R}[x_1,\dots,x_n]$ is a homomorphism out of $\mathbb{R}[x_1,\dots,x_n]$ which sends $f$ to $0$ - in other words, it is evaluation at a root of $f$. Thus $[g]\in\mathcal{I}(\mathbb{R}[x_1,\dots,x_n]/f\mathbb{R}[x_1,\dots,x_n])$ when $g$ is $0$ on all roots of $f$. Hence $\mathbb{R}[x_1,\dots,x_n]/f\mathbb{R}[x_1,\dots,x_n]$ is geometric when a polynomial being $0$ on the roots of $f$ implies it is a multiple of $f$.
Since $\mathbb{R}[x_1,\dots,x_n]$ is a UFD, writing $f$ as the product $P_1^{\alpha_1}\dots P_k^{\alpha_k}$ of powers of pairwise distinct non-constant irreducible polynomials $P_i$, we note $f$ is $0$ exactly when $P_1\dots P_n$ is. Hence a necessary condition is that $f\mid P_1\dots P_n$, so $\alpha_i=1$ for all $i$.
We also notice if the roots of $P_i$ are a subset of the roots of $P_j$, $f$ is $0$ exactly when $f/P_i$. If the quotient is geometric then $f\mid f/P_i$, so another necessary condition is that this never occurs. I'm not sure where to go from here.
This is a good question and I'm here to tell you that it's not so easy to answer in general. Your work so far is good, and your deduction that $f$ must be squarefree is correct, and there's not a whole lot more in terms of easy conditions one can give on $f$.
$\bigcap_{p\in|\mathcal{F}|} \ker p=0$ is the case exactly when $(f)=I(V_{\Bbb R}(f))$ - that is, when $\mathcal{F}$ is the coordinate algebra of a real variety in some $\Bbb R^n$. The kernel you care about is the quotient of $I(V_{\Bbb R}(f))$ (all the functions in $\Bbb R[x_1,\cdots,x_n]$ vanishing where $f$ does) by $(f)$. By the real Nullstellensatz, $I(V_{\Bbb R}(f))$ is the real radical of $(f)$, also denoted $\sqrt[\Bbb R]{(f)}$, and it may be described as the ideal $$\{ p\in\Bbb R[x_1,\cdots,x_n] \mid \exists m\in\Bbb Z_{>0}, \exists q_j\in\Bbb R[x_1,\cdots,x_n] \text{ such that } p^{2m}+\sum_j q_j^2\in I\}.$$
This is not so easy to compute in general. One source of difficulty is that if some element of $(f)$ can be written as a sum of squares $\sum_i f_i^2$, then each $f_i$ belongs to the real radical - so you have to be able to understand all the ways in which multiples of $f$ decompose in to sums of squares in order to compute this. For instance, if $f=x^2+y^2$, a perfectly respectable irreducible polynomial, then $\sqrt[\Bbb R]{(x^2+y^2)}=(x,y)$.
In general, there are algorithms available to compute $\sqrt[\Bbb R]{I}$, but I don't know the current state of the field very well - I'm under the impression that it's still an active area of research, with some new techniques showing up even this year. A cursory search of the literature gives "Computations of real radicals of polynomial ideals" by Becker and Neuhaus (doi link) and part II, by Neuhaus alone (doi link, includes full article) as the start of published work in the area, and "Computing real radicals by moment optimization" by Baldi and Mourrain (arXiv link) from earlier this year as the most modern thing I can find.