OPE coefficients of descendants in the presence of null vectors Consider the OPEs of a 2D CFT, with nondegenerate primary fields labeled by some collection $\mathtt{Prim}$ of indices. The data specifying the OPE can be taken as a collection
$$\{C_{p,q}^{r;P}\}_{p,q,r\in \mathtt{Prim},P \in \mathrm{Part}}\subset \mathbb{C},$$ where $\mathrm{Part}$ is the collection of all multisets of positive integers (including the empty multiset). Hence, an element $P\in\operatorname{Part}$ can be uniquely identified with a partition of some nonnegative number into positive parts and vice versa. (To each such partition we can associate an element $L_{P}=L_{-n_1}\cdots L_{-n_K}$, where $n_1\geq\cdots\geq n_K\geq 1$ are the elements of $P$, with multiplicity.) The complex number $C^{r;P}_{p,q}$ is then the coefficient of $L_P$ applied to the $r$th primary field in the OPE between the $p$th and $q$th primary fields.
Given this setup, we can solve for $C^{r;P}_{p,q}$ in terms of $C^r_{p,q}=C^{r;\varnothing}_{p,q}$ using conformal invariance. See [1], Section 6.6.3. However, this algorithm seems to use non-degeneracy quite centrally, as solving for the coefficients $C^{r;P}_{p,q}$ at level $N$ requires inverting the Kac-Gram matrix at level $N$. Concretely, in order to solve Eq. 6.180 in [1], one must invert the matrix on the left-hand side, which is essentially the Kac-Gram matrix. Clearly, this fails when discussing the analogous problem in the presence of degenerate fields, e.g. when working with minimal models, since the Kac-Gram matrix is non-invertible. (The denominators in Eq. 6.181 can vanish, for instance.) This is discussed briefly in [1], Section 8.A, where the connection is made to the fusion rules of minimal models. However, I do not think that book contains a proof of something like the following, analogous to the result that applies when the Kac-Gram matrix is invertible at every level:
Letting $\{C^{r;P}_{p,q}\}_{p,q,r,P}$ denote a set of indended OPE coefficients for a theory with the primary field content specified by the minimal models, each $C^{r;P}_{p,q}$ is uniquely determined by the collection $\{C^{r;\varnothing}_{p,q}\}_{p,q,r}$ (up to the coefficients of null fields), and conversely given $\{C^{r;\varnothing}_{p,q}\}_{p,q,r} \subset \mathbb{R}$ (which may be subject to some specified finite set of constraints) there exists some collection $\{C^{r;P}_{p,q}\}_{p,q,r,P}$ of numbers such that the corresponding OPE is conformally invariant in the sense of Section 6.6 of [1].
Since this seems important to the construction of minimal models via the bootstrap, I assume it is true, but it does not appear straightforward to extract a proof from the content of [1]. Is there somewhere I can find a proof of this result? Requiring that the right-hand side of the $N$th level analogue of Eq. 6.180 is in the range of the matrix on the left-hand side for all $N$ naively imposes an infinite number of constraints on the conformal dimensions $h_p,h_q,h_r$ associated to $p,q,r$ when the $r$th primary field is degenerate, so from this perspective the proposition above would seem quite miraculous.
If one defines the term "fusion rules" to refer to the triples of $h_r,h_p,h_q$ such that the proposition holds, then the question becomes to show that the restrictions on $h_r,h_p,h_q$ described in Chapter 7 pf [1] are actually the fusion rules. (E.g., why does the discussion in Section 7.3
suffice?)
[1] Di Francesco, Mathieu, and Senechal, Conformal Field Theory
 A: Since you are concerned with conformal invariance rather than unitarity, there is no need to consider a full set of Virasoro primaries which close under fusion. We can just show that one primary coefficient $C^{r,\varnothing}_{p,q}$ determines all of the descendant coefficients $C^{r,P}_{p,q}$ where the $P$'s are taken to span the Verma module that has null states modded out.
To achieve this, the Kac-Gram matrix doesn't need to be invertible... it just needs to be symmetric. In terms of the "partition function" from number theory, the $N$-th block has a size of $p(N) \times p(N)$ so we can diagonalize it level-by-level. Since the Verma module for $r$ is degenerate, some of the eigenvalues will be zero but this is where we mod out. We keep only the eigenvectors with non-zero eigenvalues and from these we get a collection of operators
$$\mathcal{O}_{r,i} = \sum_{|P| = N} \alpha_{i,P} L_P \mathcal{O}_r$$
such that
$$\left < \mathcal{O}_{r,i}(z_1) \mathcal{O}_{r,j}(z_2) \right > = \frac{\delta_{ij}}{z_{12}^{2h_r + 2N}}.$$
Now we just need to know the 3pt coefficients which appear in $\left < \mathcal{O}_p \mathcal{O}_q \mathcal{O}_{r,i} \right >$. It is tedious, but for all $L_P$ generators, we can always "take them out" of $\left < \mathcal{O}_p \mathcal{O}_q L_P \mathcal{O}_r \right >$ by repeatedly using
$$L_{-n} \mapsto \mathcal{L}_{-n} = \sum_{i \neq 3} \frac{h_i(n - 1)}{(z_i - z_3)^n} - \frac{\partial_i}{(z_i - z_3)^{n-1}}$$
which the yellow book proves starting from the Virasoro Ward identity.
A: Here's a sketch of a proof, based partly on Connor's answer.
One first uses the symmetricness of the Kac-Gram matrix to find a basis $\{J_n\}_{n=0}^\infty \cup \{\Delta_n\}_{n=0}^\infty$ of the universal enveloping algebra of the negative part of the Virasoro algebra such that in the Verma module $V(c,h_r)$, the set $\{J_n|h_r\rangle\}_{n=0}^\infty$ is orthonormal and all $\Delta_n |h_r\rangle$ are null. Now, we can prove that the proposition in the original question will hold for $p,q,r \in \mathtt{Prim}$ if and only if $\mathcal{L}_{\Delta_n} \langle \mathcal{O}_p(w) \mathcal{O}_q(w')\mathcal{O}_r(z)\rangle=0$ holds for all $n$, where $\mathcal{L}_{\Delta_n}$ is the ordinary differential operator in $z$ associated to $\Delta_n$. (Here $\langle \mathcal{O}_p(w) \mathcal{O}_q(w')\mathcal{O}_r(z)\rangle$ denotes the 3-point function associated to $h_p,h_q,h_r$, with arbitrary nonzero coefficient.)
Now, we appeal to the claim that there exist two $\Delta,\Delta'$ (which I think we can take WLOG to be among the $\Delta_n$)  such that all of the $\Delta_n$ are in the left-submodule of the Virasoro enveloping algebra generated by $\Delta,\Delta'$. Hence, the ODE $\mathcal{L}_{\Delta_n} \langle \mathcal{O}_p(w) \mathcal{O}_q(w')\mathcal{O}_r(z)\rangle=0$ holds for all $n$ if and only if the two ODEs $\mathcal{L}_{\Delta} \langle \mathcal{O}_p(w) \mathcal{O}_q(w')\mathcal{O}_r(z)\rangle=0$, $\mathcal{L}_{\Delta'} \langle \mathcal{O}_p(w) \mathcal{O}_q(w')\mathcal{O}_r(z)\rangle=0$   both hold. Each of these two is equivalent to a constraint on $h_p,h_q,h_r$, and this is discussed adequately in Section 8.3 of [1].
