This question is in some sense a follow-up to [1]: is it known how to construct the primary field operators of the unitary minimal models $\mathcal{M}(m+1,m)$ in the Coulomb gas formalism? (This would imply an answer to my original question, I think.) As far as I can tell, the original work of Dotsenko and Fateev [2] is limited to the 4-point functions and does not give a direct definition of the field operators (as honest operators, on some not-necessarily "minimal" Hilbert space(s)) themselves (though, at a physicist's level of rigor, it should fix the N-point functions via the conformal bootstrap, which ought then determine the field operators via a reconstruction theorem). The (physics?) paper [3] of Felder seems to me to strongly suggest that such a direct definition is possible. But, as far as I can tell, it itself falls a bit short of a rigorous definition.
Here is my incomplete attempt at making sense of Felder's definition in the simplest interesting case, the spin field in the critical Ising model $\mathcal{M}(4,3)$. Hopefully this will indicate what I have in mind, but unfortunately it requires a bit of explanation. First, some generalities about vertex operators. I will use the notation of [3][4]. In [3], we have defined some Verma models $\mathcal{F}_{\alpha,\alpha_0}$ of a central extension of the Heisenberg algebra that also come equipped with representations of the Virasoro algebra. Here, $\alpha \in \mathbb{R}$ is arbitrary and $\alpha_0$ is some fixed number. For each $z\in \mathbb{C}$ and $\beta \in \mathbb{R}$, we have a formal vertex operator $V_{\beta}(z)$ which is defined as a binary form $$\mathcal{F}_{\alpha+\beta,\alpha_0}\times\mathcal{F}_{\alpha,\alpha_0} \to \mathbb{C}$$ (see pg. 14 of [5]). In [5] is constructed a nested 1-parameter family $\{\mathcal{H}_{\alpha,\alpha_0,r}\}_{r>0}$ of completions of $\mathcal{F}_{\alpha,\alpha_0}$, and then it can be shown that, whenever $r<|z|<R$, $V_\beta(z)$ defines a Hilbert-Schmidt (in particular, bounded) map $\mathcal{H}_{\alpha,\alpha_0,r}\to \mathcal{H}_{\alpha+\beta,\alpha_0,R}$ depending analytically on $z$. Then, for $z_1,\ldots,z_N$ with $r<|z_N|<\cdots<|z_1|<R$ we have a composition $V_{\beta_1}(z_1)\cdots V_{\beta_N}(z_N)$ which is well-defined in the sense that it does not depend on which intermediate $\mathcal{H}_{\bullet,\alpha_0,\rho}$'s we use. One can then show that the product analytically continues to the universal cover of the moduli space $\mathcal{M}_{N,r,R}$ of $N$ distinct points in the annulus $\{z: r<|z|<R\}$. In fact, the universal cover is overkill, and it continues to the cover corresponding to the commutator subgroup of $\pi_1(\mathcal{M}_{N,r,R})$. Now here is the bit specific to the Ising model: let $$ \mathcal{H}_r = \overline{\bigoplus}_{s\in \mathbb{Z}} \mathcal{H}_{s\sqrt{3}/4,\alpha_0,r}\otimes \mathcal{H}_{s\sqrt{3}/4,\alpha_0,r},$$ where $\alpha_0 = 1/(4\sqrt{3})$, and, for $C_1,C_2>0$, $$ \sigma(z) = C_1 V_{\sqrt{3}/4}(z)V_{\sqrt{3}/4}(z^*) +C_2 \Big( \int_{P(z)}V_{-\sqrt{3}/2}(w)V_{\sqrt{3}/4}(z) \,\mathrm{d} w \Big)\Big( \int_{P(z^*)}V_{-\sqrt{3}/2}(w^*)V_{\sqrt{3}/4}(z^*) \,\mathrm{d} w^* \Big),$$ which we consider as a bounded operator $\mathcal{H}_r\to \mathcal{H}_R$ for $r<|z|<R$. Here $P(z) \in [\pi_1(\mathcal{M}_{N,r,R}),\pi_1(\mathcal{M}_{N,r,R})]$ is a Pochhammer contour with respect to $z$ and the origin. (The $V(z)$'s are acting on the left factors of the tensor products, and the $V(z^*)$'s are acting on the right factors.) This should be a well-defined $(1,2)$-primary field, and I suspect that it is (in some appropriate sense) the critical Ising spin field. So, I believe that for an appropriate choice of $C_1C_2$ (only the product shows up in correlation functions with $|0,\alpha_0\rangle \in \mathcal{F}_{0,\alpha_0}$), we have $$ \langle \alpha_0,0| \sigma(z_1)\cdots \sigma(z_N)|0,\alpha_0\rangle = \Big[ \sum_{\substack{\varsigma_1,\ldots,\varsigma_N \in \{-1,+1\} \\ \varsigma_1+\cdots+\varsigma_N = 0}} \prod_{1\leq m < n \leq N} |z_m-z_n|^{\varsigma_m\varsigma_n/2} \Big]^{1/2}. $$ This is automatically true for $N$ odd, and for $N=2$ we have, via a surprisingly nontrivial computation, $$\langle \alpha_0,0| \sigma(z_1)\sigma(z_2)|0,\alpha_0\rangle \propto C_1C_2 |z_1-z_2|^{-1/4},$$ where the constant of proportionality is explicitly computable (and nonzero), so the question is whether or not this continues to hold for $N=4,6,8,\cdots$. For any $N\in 2\mathbb{N}^+$, the left-hand side is some explicit iterated contour integral, but I do not know how to compute it. I attempted to test this for $N=4$ using some sketchy numerics and got agreement within ~3%, but Mathematica's NIntegrate was not very happy with the $\smash{(z-w)^{-3/4}}$'s that show up in the integrals...
Update: I did the numerics for the 4-point function properly, and there is complete agreement in this case.
References:
[1] Wightman QFTs corresponding to minimal models?
[2] Dotsenko and Fateev, Conformal algebra and multipoint correlation functions in 2D statistical models
[3] Felder, BRST Approach to Minimal Models
[4] Francesco, Mathieu, and Senechal, Conformal Field Theory, Section 9.B
[5] Bovier, A rigorous treatment of conformal blocks in a model of bosonic conformal field theory
