A partial answer:
Let $\psi \in \mathcal D(N)$; in particular, it follows that $\psi \in \mathcal D(N^{1/2})$. We want to show that these conditions are sufficient such that e.g. $a(f)\psi \in \mathcal D(N^{1/2})$. Then the product $a^\dagger(g) a(f) \psi$ would be well-defined.
To do so, I think we can use the fact that for all $\psi \in \mathcal D(N^{1/2})$ it holds that$^\ddagger$
$$||a(f)\psi||_F \leq ||f||_{\mathfrak h}\, ||N^{1/2}\psi||_F \quad ,\tag{1}$$
where $\mathfrak h$ denotes the underlying one-particle space and $f\in \mathfrak h$. We proceed by noting that for $M\in \mathbb N$ we have \begin{align} ||f||^2_\mathfrak h \sum\limits_{n=1}^{M+1} n^2 ||\psi_n||^2_n = ||f||^2_\mathfrak h \sum\limits_{n=1}^{M+1} n \underbrace{||N^{1/2}\psi_n||^2_n}_{=||N^{1/2}\psi_n ||_F^2} &\geq \sum\limits_{n=1}^{M+1}n\underbrace{||a(f)\psi_n||^2_{F}}_{=||(a(f)\psi)_{n-1}||_{n-1}^2} \\ &=\sum\limits_{n=0}^M n||(a(f)\psi)_n||^2_n + \sum\limits_{n=0}^M ||(a(f)\psi)_n||_n^2 \quad , \end{align} where in some intermediate steps we have identified $\psi_n$ with a vector in $F$ in an obvious way. We thus find
$$ ||f||^2_\mathfrak h \sum\limits_{n=0}^{M+1} n^2 ||\psi_n||^2_n \geq \sum\limits_{n=0}^M n||(a(f)\psi)_n||^2_n \tag{2} \quad .$$
For $\psi \in \mathcal D(N)$, the LHS of $(2)$ converges for $M\to \infty$ and therefore the RHS converges too, which shows that indeed $a(f)\psi \in \mathcal D(N^{1/2})$.
I further think that the domains $\mathcal D(N^{K/2})$ for all $K\in \mathbb N$ are dense subspaces, since the dense subspace $F_{\mathrm{fin}}\subset F$, consisting of all vectors $\psi=(\psi_n)_n$ with only finitely many non-zero $\psi_n$, is a subspace of these domains, i.e. we have $ F_{\mathrm{fin}} \subset \mathcal D(N^{K/2}) \subset F$ and the denseness of $F_{\mathrm{fin}}$ in turn implies the denseness of $\mathcal D(N^{K/2})$.
$^\ddagger$ It remains to show that $(1)$ holds in the said domain without using some result we want to prove. In Ref. 1. a proof is sketched, which seems plausible to me.