Gaudin Models and Multipoint Conformal Blocks: General Theory

The construction of conformal blocks for the analysis of multipoint correlation functions with $N>4$ local field insertions is an important open problem in higher dimensional conformal field theory. This is the first in a series of papers in which we address this challenge, following and extending our short announcement in [Phys. Rev. Lett. 126, 021602]. According to Dolan and Osborn, conformal blocks can be determined from the set of differential eigenvalue equations that they satisfy. We construct a complete set of commuting differential operators that characterize multipoint conformal blocks for any number $N$ of points in any dimension and for any choice of OPE channel through the relation with Gaudin integrable models we uncovered in [Phys. Rev. Lett. 126, 021602]. For 5-point conformal blocks, there exist five such operators which are worked out smoothly in the dimension $d$.


Introduction and Summary of Results
In the last few decades conformal field theories have steadily gained importance, first in d = 2 and then for d > 2. They provide a powerful source of new paradigms that advance our understanding of quantum field theory deep in the quantum regime which is inaccessible to perturbation theory.
Through the AdS/CFT correspondence conformal field theories have even begun to teach us about some of the deepest mysteries of quantum gravity.
In a conformal field theory, Poincare symmetry is enhanced to conformal symmetry. This symmetry enhancement turns out to be highly constraining, even in d > 2, where the conformal group is finite dimensional and only contains d + 1 generators outside of the Poincare group. These few additional generators provide enormous control over the operator product expansion (OPE) of local fields and in fact fix these products up to some numerical OPE coefficients. The latter determine all higher correlations, at least in principle. In practice, however, higher correlations with more than N = 3 field insertions possess an intricate dependence on the insertion points of the local fields and their simplicity only becomes visible after expanding correlators in a basis of conformal blocks. The latter are very much like plane waves in ordinary Fourier analysis, i.e. conformal blocks (or rather the closely related conformal partial waves) are to conformal symmetry what plane waves are to the translation group.
Once the basis of conformal blocks is known, one can expand the correlation function. Conformal symmetry then implies that the coefficients are simply products of the coefficients that appear in the OPE.
All this is of course well known since the early days of conformal field theory, see [2][3][4]. But the pioneers of conformal field theory were only able to write down integral formulas for conformal blocks [5,6]. These so-called shadow integrals are very much like Feynman integrals. In particular, it is a considerable challenge to relate shadow integrals to some known special functions, to read off their properties and to find relations between them, at least for d > 2 where the integrations cannot be performed directly. There has been little progress on this problem until Dolan and Osborn shifted attention from the integrals to the differential equations that these integrals satisfy [7,8]. Once they had constructed their so-called Casimir differential operators they were able to extract a wealth of information on 4-point blocks, see also [9][10][11]. This work has been the decisive input for the modern conformal bootstrap program (see [12] for a review and many references to the original literature) and its spectacular applications to the d = 3 Ising model, in particular [13][14][15][16]. More recently, it was observed that the Casimir differential operators studied by Dolan and Osborn could be identified with the Hamiltonians of an integrable 2-particle quantum mechanics system of Calogero-Sutherland type [17]. The observation was later explained in the context of harmonic analysis for the conformal group [18][19][20][21]. Calogero-Sutherland Hamiltonians and their eigenfunctions have actually been studied in mathematics for several decades where it has become instrumental in developing the modern theory of multivariate hypergeometric functions. After the relation with conformal blocks was noticed it became clear that most of the results on 4-point blocks had been known long before the work in conformal field theory, mostly from the early work on Calogero-Sutherland models by Heckman and Opdam [22], see [23].
In principle it is possible to probe a conformal field theory by studying the entire set of 4-point functions for the infinite number of (spinning) conformal primary fields. And even though the theory of spinning conformal blocks is reasonably well developed by now, see e.g. [24][25][26][27], working with an infinite number of correlation functions seems impossible at least in the absence of any additional algebraic structure. In practice, it may appear more promising to study higher correlations of a small set of fundamental fields of the theory. The 8-point function of a fundamental scalar field, for example, already contains as much dynamical information as an infinite set of spinning 4-point functions. Given that it took more than 30 years to develop a theory of 4-point blocks, it may seem like a daunting task to actually extend this theory to multi-point blocks, though there has been significant activity in this direction lately, see e.g. [28][29][30][31][32][33][34][35][36][37][38][39][40][41]. The embedding of 4-point blocks into harmonic analysis of the conformal group and the profound relations with integrable quantum mechanical models may be viewed as a strong hint towards the appropriate mathematical tools for the study of multi-point blocks. And in fact, as we have announced in [1], the relation to integrable models is not restricted to the case of 4-point blocks. As we have outlined, multi-point blocks for any number N of local fields can be characterized as joint eigenfunctions of a complete set of commuting differential operators. The latter were argued to arise as Hamiltonians of Gaudin integrable models.
Our goal in this paper and its sequels is to substantiate our claims and to develop them into a full theory of multi-point blocks. In this paper we perform the first step of this program and explain in detail the relation between multi-point conformal blocks and Gaudin integrable systems. We will actually go further than [1] which only sketched the construction of differential operators for the socalled comb channel conformal blocks [28]. The results described in this work apply to any channel, i.e. they provide a complete set of commuting differential operators for whatever channel one is interested in, including e.g. the snowflake channel that has also received some attention lately, [34,36]. As we shall explain below, in any given channel, the differential operators split into two groups: the Dolan-Osborn-like Casimir differential operators that measure the weight and spin of intermediate fields, and the novel vertex differential operators that measure a choice of tensor structure in the operator product of intermediate (spinning) fields. In a second paper of this series we will address such internal vertices for the comb channel of 3-and 4-dimensional conformal field theories. These vertices give rise to a single vertex differential operator, see [1] for an example, that we will identify as the Hamiltonian of an elliptic Calogero-Moser system for some complex crystallographic reflection group, first obtained in [42]. Extensions to vertices with more than one degree of freedom, which are relevant e.g. for the snowflake channel in 4D, as well as a theory of solutions for single and multi-variable vertices will be addressed in the future.
The main results of this work are the following. We construct the Casimir and vertex differential operators that are simultaneously diagonalized by conformal blocks. Commutativity of this set of operators is established by realizing it as a limit of the Gaudin integrable model. The connection to the Gaudin model seems to be necessary for the full proof, although some parts of the proof can be established by elementary arguments. Further, we provide evidence that independent Casimir and vertex operators are equal in number to independent cross ratios and can thus be used to completely characterize conformal blocks. This is done by exhibiting a number of relations between the operators, leaving us with an explicit commuting system with no apparent further dependencies. Finally, we construct the five commuting Casimir and vertex operators for 5-point functions explicitly as differential operators in the five cross ratios.
For the remainder of this introduction we will give a more precise summary of the main results and place them in the context of conformal field theory.

Summary of Results
In this work we consider correlation functions of N local scalar primary fields φ i with conformal weights denoted by ∆ i , (1.1) The fields are inserted at points x i , i = 1, . . . , N , of d-dimensional Euclidean space R d . Scalar primary fields are characterized by their commutation relations with the generators of the conformal algebra which take the form where T α are the usual first order differential operators and α runs through the set of conformal generators.
Prerequisites and notation. Conformal correlation functions are symmetric with respect to the exchange of any two fields. On the other hand, the way we compute them is not. In fact, N -point correlation functions are evaluated by a repeated use of operator product expansions (OPEs) and the precise sequence of performing these expansions determines what is known as an OPE channel. These channels can be represented as (plane) tree diagrams C N OPE with enumerated leaves. Figure 1 shows one such example for the 10-point function. Given a fixed enumeration of the external scalar fields, the number of OPE channels is (2N − 5)!!. Each of these channels is associated with a system of conformal blocks that provide a Fourier-like expansion of the correlation function. After stripping off some appropriate prefactors Ω(x i |∆ i ), the conformal blocks depend only on conformally invariant cross ratios. It is well known that one can form of such independent cross ratios from a set of N points x i ∈ R d . Note that this number is independent of the OPE channel. The conformal blocks we want to expand our correlation functions in must depend on the same number of labels. For example, in the case of a 4-point function in d > 1 one has two cross ratios. The associated blocks are famously labeled by the conformal weight ∆ and the spin l of the field that is exchanged in the intermediate channel. Since the latter transforms in a symmetric traceless tensor representation of the rotation subgroup, a single number l is sufficient to characterize the spin.
When dealing with higher correlation functions, it is easy to see that the quantum numbers of fields exchanged in the intermediate channels are not sufficient. The precise number of such intermediate field labels does depend on the channel topology, at least for N > 5, but it is always strictly smaller than the number n cr of cross ratios. As an example let us consider N = 5. In this case there exist 15 OPE channels but they all possess the same topology. In Figure 2 we have displayed one of the fifteen OPE channels. All other 14 channels are related to this one by a permutation of the leaves, modulo symmetries of the bare OPE diagram that has been stripped of its leaves. As one can readily see, and are referred to as Casimir differential operators, since they are straightforward generalizations of the Casimir differential operators constructed for N = 4 by Dolan and Osborn. Second, differential operators that measure choices of tensor structure, the first example of which was introduced recently in [1], are referred to as vertex differential operators. Let us note that for scalar blocks, the choice of tensor structures and hence the vertex differential operators are relevant as soon as multiple non-scalar exchanges are involved. These types of blocks have only been considered very recently in [41] and [38, Appendix E] for five-point blocks, and in a certain limit in [43] for five or six scalar legs.
Before we can describe our results on both types of differential operators we need to set up some notation. Given an OPE channel we enumerate internal lines by Latin indices r = 1, . . . , N − 3 and vertices by Greek indices ρ = 1, . . . , N − 2. There are no rules on how to enumerate these two sets of objects. OPE diagrams are (plane) trees and hence by cutting any internal line with label r we separate the diagram into two disconnected pieces. Hence, r is associated with a partition of the external fields into two disjoint sets, Similarly, any vertex ρ gives rise to a partition of N into three disjoint sets Given any subset I ⊂ N we can define the following set of first order differential operators in the insertion points x i , Let us note that for two disjoint sets I 1 , I 2 ⊂ N we have Casimir differential operators. With this notation it is very easy to construct the differential operators that measure the quantum numbers of the intermediate fields, The number r d = [(d + 2)/2] denotes the rank of the conformal Lie algebra. In even dimensions d, the symmetric invariant tensor κ p of order p = 2r d = d + 2 actually possesses a square root of order p = r d that also commutes with all generators of the conformal algebra. This so-called Pfaffian differential operator has the same form as in (1.8), but with a symmetric invariant tensor κ p of order p = d/2 + 1, (1.10) When d = 4k + 2, the symmetric invariants of order p = d/2 + 1 are twofold degenerate and we should use two different symbols for these two invariants of order d/2 + 1. In order not to clutter notation too much, we decided to ignore this distinction. In other words, we will consider κ d/2+1 as a pair of symmetric invariants when d = 4k + 2.
In our formulas for the differential operators we have placed a subscript |G to stress that they are defined as operators acting on correlations functions, i.e. on functions G that satisfy the conformal Ward identities In our construction of the differential operators we have favored the set I r,1 over I r,2 . But from the conformal Ward identities we can conclude that Though some caution is needed when we apply this relation to the evaluation of the Casimir differential operator, see Subsection 2.1 for details, it is not difficult to see that all differential operators of even order come out the same if we pick I r,2 rather than I r,1 . There is only one family for which the set matters, namely for the Pfaffian operators when d is a multiple of four. In that case the operator flips sign when we change the set. Of course, overall factors are a matter of convention and hence of no concern. Therefore, we shall drop the reference to the set we use in the construction of Casimir differential operators, writing D p r instead of D p r,1 . An important point to note is that the Casimir differential operators need not be independent. To illustrate this, consider the case N = 4 for d > 2. Since all external fields are assumed to be scalar, the single intermediate field is a symmetric traceless tensor and is hence characterized by two numbers only, its weight ∆ and spin l. These can be measured by the Casimir differential operators D 2 and D 4 .
But starting from d = 4, the conformal algebra possesses Casimir elements of higher order which are independent in general, but become dependent on the lower order ones when evaluated on symmetric traceless tensors. More generally, the number of independent Casimir differential operators at a given internal line r is given by 12) and |I| denotes the order of the set I. 1 We can determine the total number of independent Casimir differential operators to be Let us note that the total number of Casimir differential operators does depend on the topology of the OPE channel, not just on the number N of points. In the case of N = 6 and d ≥ 4, for example, there are n cdo (C N =6 comb , d) = 7 Casimir differential operators in the comb channel, while the snowflake channel admits only n cdo (C N =6 snowflake , d) = 6 of such operators.
Vertex differential operators. What we have described so far is nothing new, and can be established by elementary means. But as we have explained, starting from N = 5 the Casimir differential operators do not suffice to resolve all quantum numbers of the conformal blocks, i.e. n cdo is strictly smaller that n cr for all OPE channels. Our main task is to construct additional differential operators that can measure the choice of tensor structures at the vertices independently of the weights and spins of the intermediate fields, i.e. we need to find a complete set of vertex differential operators that commute with the Casimir differential operators and among themselves. In this work we describe how to accomplish this task, for any number N of external scalar fields and any OPE topology. One central claim is that these vertex differential operators take the form D p,ν ρ,12 = κ α1,...,αν ,αν+1,...,αp p T (Iρ,1) α1 · · · T (Iρ,1) αν T (Iρ,2) αν+1 · · · T (Iρ,2) αp |G (1.14) where ν = 1, . . . , p−1 and p = 2, 4, . . . , 2r d = d+1 when d is odd. For even d, we let p run through even integers until we reach d and add a set of Pfaffian vertex operators Pf ν ρ,12 , ν = 1, 2, . . . d/2 which are constructed with a symmetric invariant tensor κ p of order p = d/2 + 1. Let us note that the definition of all these vertex differential operators also makes sense for ν = 0 and ν = p. The corresponding objects coincide with Casimir differential operators for the links that enter the first and second leg of the vertex. This is why we have excluded them from our list. The remaining operators still allow us to reconstruct the Casimir operators for the link that enters the third leg. Therefore, there is one linear relation for each value that p can assume, i.e. we have r d linear relations in total. One may use these relations to eliminate e.g. the operator with ν = p/2.
Let us note that the definition of the vertex operators D p,ν ρ,ij depends on the choice of labeling of the subsets I ρ,j forming the partition N = I ρ,1 ∪ · I ρ,2 ∪ · I ρ,3 associated with the vertex ρ, which is arbitrary. However, the algebra generated by the vertex operators D p,ν ρ is in fact independent of this choice: more precisely, the vertex operators constructed from another choice of labeling of the I ρ,j 's are linear combinations of the operators D p,ν ρ,12 , modulo the use of the conformal Ward identities (1.11), see Section 2.2.
The number of vertex differential operators at a given vertex is now easy to count. Taking into account that one additional linear relation among the operators listed in eq. (1.14), one finds The first key result of this work is that these vertex differential operators commute among themselves and with the Casimir differential operators. Commutation between Casimir and vertex differential operators is obvious. Similarly, it is easy to show that two vertex operators commute if they are associated with different vertices ρ = ρ . The deepest part of our claim concerns the fact that also vertex operators associated with the same vertex commute. It does not seem straightforward to prove this statement by elementary manipulations. Below we shall use an indirect strategy in which we identify these vertex differential operators with Hamiltonians of some Gaudin integrable system defined on a 3-punctured sphere. For the latter, commutativity has already been established.
Of course the vertex differential operators we listed may not all be independent, as for the Casimir operators, see discussion above. In order to count the number of independent vertex differential operators, we shall employ the depth function d(I, d) we introduced in eq. (1.12). For a given vertex ρ inside an OPE channel C N OPE , the number of independent vertex differential operators is expected to be equal to the degrees of freedom associated to this vertex where d ρ,i = d(I ρ,i , d) with i = 1, 2, 3. The inequality is saturated for vertices ρ with d ρ,i = r d . For the special vertices that can appear in the comb channel and in which one of the legs is scalar, the formula becomes Here ρ m is a vertex with d ρm,1 = m, d ρm,2 = 1 and d ρm,3 = m + 1, see Figure 3, and ρ d is the maximal comb channel vertex with d ρr d ,1 = r d = d ρr d ,3 . The total number n vdo (C N OPE , d) of vertex differential operators is obtained by summing over all N − 2 vertices, i.e.
At least for the comb channel, it is easy to verify that the number of independent Casimir and vertex differential operators coincides with the number of cross ratios, The formula holds of course for all OPE channels. Below we shall exhibit the relations among vertex differential operators that are responsible for the reduction from the n v (d) operators in our list (1.14) (with ν = p/2 removed) to the n vdo,ρ (C N OPE , d) independent vertex differential operators that are needed to characterize the vertex ρ. This is the second key result of this work. It will allow us in particular to determine the precise order of each independent vertex differential operator. While Gaudin models for the 3-punctured sphere only enter the discussion as a convenient tool to construct commuting vertex differential operators at the individual vertices, the relation between conformal blocks and Gaudin models turns out to reach much further. In fact, is is possible to embed the whole set of Casimir and vertex differential operators for arbitrary scalar N -point functions into Gaudin models on the N -punctured sphere. The latter contains N additional complex parameters that are not present in correlation functions. In the Gaudin integrable model these parameters correspond to the poles of the Lax matrix and they enter all Gaudin Hamiltonians. By considering different limiting configurations of these parameters is is possible to recover the full set of Casimir and vertex differential operators for scalar N -point functions in all the OPE channels. This construction not only embeds our differential operators into a unique Gaudin integrable model, but also shows that operators in different channels are related by a smooth deformation.
Let us finally outline the content of the following sections. Section 2 is mostly devoted to the study of the individual vertices. After a brief discussion of commutativity for Casimir differential operators and also vertex differential operators assigned to different vertices, we shall zoom into the individual vertices for most of Section 2. In Section 2.2 we construct the vertex differential operators in terms of the commuting Hamiltonians of a 3-site Gaudin integrable system. Section 2.3 addresses the relations between these operators for restricted vertices. The main purpose of Section 3 is to embed the whole set of Casimir and vertex differential operators for arbitrary scalar N -point functions into Gaudin models on the N -punctured sphere. In Section 4 we discuss one concrete example, namely we construct all five differential operators that characterize the blocks of a scalar 5-point function in any d ≥ 3. Two of these operators are of order two while the other three are of fourth order. The paper concludes with a summary, outlook to further results and a list of interesting open problems to be addressed.

The Vertex Integrable System
The aim of this section is to address the key new element in the construction of multi-point conformal blocks for d ≥ 3: the vertices themselves. In the first subsection we shall show that the construction of commuting differential operators for scalar N -point blocks can be reduced by rather elementary arguments to the construction of commuting differential operators for 3-point functions of spinning fields. Recall that the dependence on spin degrees of freedom can be encoded in auxiliary variables from which one is often able to construct non-trivial cross ratios, 2 even in the case of a 3-point function.
Constructing sufficiently many commuting vertex differential operators that act on such cross ratios of the 3-point function requires more powerful technology from integrability which we shall turn to in the second subsection. There we construct commuting differential operators for vertices from the Gaudin integrable model for the 3-punctured sphere. The basic construction provides n v (d) of such commuting operators and hence sufficiently many even for the most generic vertices. For special vertices, such as those appearing in the comb channel with external scalars, there exist linear relations between these operators. These are the subject of the third subsection.

Reduction to the vertex systems
Our goal here is to prove that Casimir and vertex operators constructed around different vertices of an OPE diagram commute. More precisely we shall show that for every pair (r, p), (r , q) of Casimir operators and any choice (ρ, q, ν) of a vertex operator, including the Pfaffian Casimir and vertex operators that appear at order p, q = d/2 + 1 when the dimension d is even. Note that the individual vertex differential operators also depend on the choice a ∈ {12, 23, 13} of a pair of legs. In addition, we shall also establish that vertex operators associated with different vertices commute, and all triples (p, ν, a) and (q, µ, b). This leaves only the commutativity of operators attached to the same vertex which is deferred to the next subsection.
The properties (2.1) and (2.2) are in fact elementary. They require global conformal invariance, the tree structure of OPE diagrams and the commutativity property (1.7). To begin with, let us recall that we associated two disjoint sets I r,1 and I r,2 to every link. As we pointed out before, the Casimir differential operators (1.8) do not depend on whether we used the generators T (Ir,1) α or T (Ir,2) α to construct them. Let us briefly discuss the details of the proof. Note that we think of D as an operator acting on correlation functions G. This is signaled by the subscript |G in the definition of the Casimir differential operators. In evaluating products of first order differential operators, one can only apply the Ward identity (1.11) to the rightmost operator, which acts directly on the correlation function, and not on some derivative thereof. But once we have converted the rightmost operators T (Ir,1) into −T (Ir,2) , they will commute with all operators to their left, such that we can freely move them all the way to the left and proceed to apply the Ward identity to the next set of first order operators, and so on. If we finally take into account that the invariants κ of the conformal Lie algebra are symmetric, we arrive at an expression for the Casimir differential operators in terms of T (Ir,2) .
A similar analysis can be carried out for vertex operators, see also Subsection 2.2. Without loss of generality we can assume that we have constructed our vertex operators in terms of the generators T (Iρ,1) and T (Iρ,2) and want to switch to constructing them from T (Iρ,1) and T (Iρ, 3) instead. To do so we make use of the invariance condition that follows from relations (1.5) and (1.11). After we apply this to the rightmost operator in the vertex differential operator, we use the commutativity property (1.7) to move the generators T (Iρ, 1) and T (Iρ, 3) to the left of T (Iρ,2) . We continue this replacement process until all the generators T (Iρ, 2) are removed and using symmetry of the tensor κ we find along with a similar relation for the Pfaffian vertex operators for p = d/2 + 1 and d even. Since the last term in this sum with µ = p − ν is just a Casimir operator, we have managed to express all vertex differential operators that are constructed from the generators associated with I ρ,1 and I ρ,2 as a linear combination of the vertex operators associated with the pair I ρ,1 and I ρ,3 and a Casimir operator.
This is the main input in proving the commutativity statements (2.1) and (2.2).
Let us start with a pair of links r, r . Each of these links divides the set of external points into the two disjoint sets I r,1 , I r,2 and I r ,1 , I r ,2 respectively. Since the OPE diagram is a tree, it is always possible to find a pair i, j = 1, 2 such that I r,i ∩ I r ,j = ∅. For this choice because of the commutativity property (1.7). This proves our first claim. Note that the same arguments also apply to the case in which r = r and also if one or both operators are Pfaffian, i.e. if d is even and p, q = d/2 + 1.
ρ Let us now extend this argument to include vertex differential operators. In order to prove that the Casimir operators associated to a link r commute with the vertex operators associated to any vertex ρ we recall that any choice of a vertex ρ on an OPE diagram divides the diagram into the three distinct branches that are glued to the vertex, and we denote these by I ρ,j as in Figure 4. A quick glance at Figure 4 suffices to conclude that given r and ρ it is possible to find a pair i ∈ 1, 2 and j ∈ 1, 2, 3 such that I r,i ⊂ I ρ,j , since the link r must be in one of the three branches. It follows that Commutativity of Casimir and vertex differential operators then follows since we can construct the Casimir differential operators in terms of the generators for I r,i while using the generators for I ρ,j1 and I ρ,j2 for the vertex differential operators.
Let us finally consider any two distinct vertices ρ and ρ on an OPE diagram. As we highlight in Figure 5, any configuration of two vertices divides the diagram into five parts: four external branches Figure 5. Schematic representation of a generic OPE diagram with focus on two internal vertices. Operators supported around distinct vertices trivially commute, as they can be written in terms of generators that belong to different branches.
I ρ,j , I ρ ,j , j = 1, 2, attached to only one vertex, and the central part C of the diagram that is attached to both vertices ρ and ρ . Following what we just claimed with focus on one vertex, we can use diagonal conformal symmetry to rewrite the operators around ρ and ρ to depend on disjoint sets of legs I ρ,j , I ρ ,j with j = 1, 2. Since generators associated to these sets commute (1.7), it follows automatically that operators constructed around different vertices must commute as well.
This implies that to prove commutativity of our set of operators, we can just focus on operators that live around one single vertex. To prove the commutativity of these vertex operators we will now make use of the integrability technology that is provided by Gaudin models.

The vertex system and Gaudin models
In this section, we will explain how the operators (1.14) associated with a vertex ρ in the OPE diagram naturally arise from a specific Gaudin model, which in particular will provide us with a proof of their commutativity. Let us start by reviewing briefly how Gaudin models are defined [44,45].
They are integrable systems naturally constructed from a choice of a simple Lie algebra g. Having in mind applications of these systems to conformal field theories, we will choose g to be the conformal Lie algebra so(d + 1, 1) of the Euclidean space R d , with basis T α as in the previous section. The Gaudin model depends in general on M complex numbers w j , called its sites, to which are attached M independent representations of the algebra g. To obtain the vertex system we address in this section, we restrict our attention here to the case M = 3 and associate with these three sites the representations of g corresponding to the three fields attached to the vertex ρ. More precisely, using the notation defined in Section 1 and in particular the partition N = I ρ,1 ∪ · I ρ,2 ∪ · I ρ,3 constructed from the vertex ρ, we will attach to the three sites w j , j = 1, 2, 3, of the Gaudin model the generators which define representations of g in terms of first-order differential operators in the insertion points x i .
A key ingredient in the construction of the Gaudin model is its so-called Lax matrix, whose components in the basis T α are defined here as where z is an auxiliary complex variable called the spectral parameter. In the above equation, we have denoted the Lax matrix as L ρ α (z) to emphasize that this is the matrix corresponding to the vertex ρ. For any elementary symmetric invariant tensor κ p of degree p on g, there is a corresponding z-dependent Gaudin Hamiltonian of the form where . . . represent quantum corrections, involving a smaller number of components of the Lax matrix L ρ α and their derivatives with respect to z. These corrections are chosen specifically to ensure that the Gaudin Hamiltonians commute for all values of the spectral parameter and all degrees: The existence of such commuting Hamiltonians was first proven in [46], using some previously established results [47] on the so-called Feigin-Frenkel center of affine algebras at the critical level. The explicit expression for the quantum corrections was obtained in [48,49] for Lie algebras of type A and in [50] for types B, C, D, which is the case we are concerned with here (indeed, g = so(d + 1, 1) is of type B for d odd and of type D for d even). We refer to [51] for a summary of these results.
The properties of the Feigin-Frenkel center were further studied in the recent work [52], the results of which imply that the quantum corrections in H (p) ρ (z) are sums of terms of the form where q < p, τ α1···αq is a completely symmetric invariant tensor of degree q on g and r 1 , · · · , r q are positive integers such that r 1 + · · · + r q = p. For what follows, it will be useful to consider the leading part of the Hamiltonian (2.7) alone, without quantum corrections, which we will denote as Let us finally note that the quantum corrections are absent for both the quadratic Hamiltonian H To make the link with the vertex operators defined in Section 1, we will make a specific choice of the parameters w j of the Gaudin model. More precisely, we set In particular, the Lax matrix (2.6) reduces to Let us now study the Gaudin Hamiltonians H (2.10), we simply find where D p,ν ρ,12 is the vertex operator defined in eq. (1.14). To obtain this expression, we have used the fact that T to relabel the Lie algebra indices as in eq. (1.14). Noting that the fractions z −ν (z − 1) ν−p for ν = 0, · · · , p are linearly independent functions of z, it is then clear that one can extract all the vertex operators D p,ν ρ,12 from H (p) ρ (z). Note that the "extremal" operators D p,0 ρ,12 and D p,p ρ,12 coincide with the Casimir operators of the fields at the branches I ρ,1 and I ρ,2 of the vertex. The same equation also holds for the Pfaffian operators H Our goal in this section is to prove the commutativity of the vertex operators D p,ν ρ,12 using known results on Gaudin models. This would follow automatically from the commutativity (2.8) of the Gaudin Hamiltonians H (p) ρ (z) if these Hamiltonians contained the operators D p,ν ρ,12 . But we have already proven above that the latter are naturally extracted from the leading parts H with prefactors composed of powers of z and z − 1, and where τ α1···αq is a completely symmetric invariant tensor on g of degree q < p, as in eq. (2.9). In particular, τ α1···αq decomposes as a product of elementary symmetric invariant tensors κ k , symmetrized over the indices α i , with k ≤ q < p. The correction in the above equation can thus be re-expressed as an algebraic combination of lower-degree vertex operators D k,ν ρ,12 . Since these are the coefficients of the non-corrected Hamiltonian H recursion on the degrees shows that the quantum corrections can indeed be expressed in terms of lower-degree Hamiltonians, as anticipated.
Let us end this subsection with a brief discussion on the role played by the choice of labeling of the branches I ρ,1 , I ρ,2 and I ρ,3 attached to the vertex ρ. As mentioned in Section 1, this choice is arbitrary but enters the definition (1.14) of the vertex operators D p,ν ρ,12 , which in this case contain only the generators T 1 Acting on the correlation function G N , which satisfies the Ward identities (1.11), this Lax matrix then

Restricted vertices and relations between vertex operators
In the previous subsection we have shown that all of the operators listed in eq. (1.14) commute with each other. As we have pointed out before, we did not include operators with ν = 0 and ν = p in the list since these coincide with Casimir differential operators, Here r i denotes the link that is attached to the i th leg of the vertex ρ, i.e. for which I ri,j = I ρ,i with either j = 1 or j = 2. The remaining p − 1 operators satisfy one more linear relation since Let us note that this relation also applies to the Pfaffian vertex operators that exist for p = d/2 + 1 when d is even. We have used this relation to drop one of the vertex differential operators. Once these obvious relations are taken care of, the total number of commuting vertex differential operators is given by eq. (1.15) and matches precisely the maximal number of cross ratios that can be associated to a single (generic) vertex, see upper bound of eq. (1.16) in the introduction. But restricted vertices carry fewer variables, so their corresponding differential operators (constructed in the previous section) must obey further relations. It is the main goal of this subsection to discuss these relations. We will also check that, once these are taken into account, the number of remaining vertex differential operators matches the number (1.16) of cross ratios at restricted vertices.
Our arguments are based on an important auxiliary result concerning the differential operators T (I) α that are associated to some subset I ⊂ N of order |I| ≤ N/2. To present this requires a bit of preparation. Up to this point there was no need to spell out the precise form of the symmetric invariant tensors κ p that we used to construct our differential operators. Now we need to be a bit more specific. As is well known, such tensors can be realized as symmetrized traces, where η AB is the Minkowski metric with signature (d + 1, 1) and η AB is its inverse. This makes it now easy to compute κ p explicitly. The only issue arises in even d. In this case the symmetrized traces in the fundamental representation do not generate all the invariants. In order to obtain the missing invariant, one has to include the trace in a chiral representation. The standard construction employs the spinor representation in which generators T [AB] are represented as where γ A are the d + 2-dimensional γ matrices and the matrix indices are σ, τ = 1, . . . , 2 d/2+1 . One can then project to a chiral spinor representation with the help of γ c ∼ γ 0 · · · γ d+1 .
Let us now introduce the symbol T (I) to denote the following Lie-algebra valued differential operators Upon evaluation in some finite-dimensional representation, such as the fundamental or the spinor representation, these become matrix valued differential operators. With this notation we write our set (1.14) as where we take the trace in the spinor representation and include the factor γ c in the argument.
In finding relations between the vertex differential operators for restricted vertices we actually work with the total symbols of the differential operators rather than the operators themselves. This means that we replace the partial derivatives ∂ (i) µ with commuting coordinates p i µ . The associated matrices of functions of x µ i and p i µ will be denoted byT (I) . As before we shall add a subscript f, s to denote the matrices in the fundamental and the spinor representation. After passing to the total symbol the entries of the matrices commute and we can drop the symmetrization prescription when taking traces.
As a result, the total symbols of the vertex differential operators are simply traces of powers of the matricesT (I) .
We are now ready to state the main result needed to elucidate the relations between vertex differential operators. It concerns the matrix elements of the n th power of the matricesT We are now prepared to discuss relations between vertex differential operators. Let us consider a vertex ρ inside our OPE diagram. As we have explained before, ρ splits the set N into three subsets I ρ,i with i = 1, 2, 3. Each of these sets determines an integer d i = d(I ρ,i , d). Let us suppose that we construct the vertex differential operators using T (Iρ,i) for i = 1, 2 as in eq. (2.21). If one of the integers d 1 or d 2 is smaller than r d we immediately obtain relations among the vertex differential operators.
In fact, when applied to the matricesT (Iρ,1) , our claim (2.23) implies that all operators we obtain when ν > 2d 1 can be expressed in terms of Casimir and vertex differential operators of lower order.
But this does not yet include the full set of relations that appears whenever d 3 is smaller than min(d 1 +  Figure 6. Snowflake OPE diagram. Here the tensor product of any two branches around the vertex ρ would allow a mixed symmetry tensor of depth d = 4 to appear in d ≥ 7, but diagonal conformal symmetry constrains this to match with the symmetric traceless tensor produced on the third branch. To take restrictions from the third branch I 3 into account, it is sufficient to use the total symbol of eq. (2.3) and impose once more our dependence statement (2.23) for product matrices, this time for 3) . This tells us that the matrix elements of powers with n > 2d 3 can be written in terms of lower order terms. In order to convert this observation into relations among the vertex differential operators of order p, we can multiply the expression (2.25) with some appropriate powers ofT (Iρ,i) , i = 1, 2, and consider the matrix elements of the products T (Iρ,1) +T (Iρ,2) n T (Iρ,2) ν T (Iρ,1) p−n−ν (2.26) for n > 2d 3 , any allowed value of p > n and ν = 0, . . . , p − n. By binomial expansion, we can write the expression (2.26) as a linear combination of our basic vertex differential operators. After taking relations on the branches I ρ,1 and I ρ,2 into account, we obtain an additional nontrivial relation from the third branch I ρ, 3 . For example, in the cases where n is odd, contracting (2.25) with T (Iρ,2) B A leads to a relation between vertex and Casimir operators of order n + 1 and further lower order operators.
This effectively reduces the amount of vertex operators at order p by up to p − 2d 3 − 1, though the actual number can be lower in case there are less than p − 2d 3 − 1 vertex operators at order p left after imposing the constraints from the first two branches. If we are interested in counting the number of vertex operators at a given even order p, the procedure we just outlined is summarized in the following counting formula where Θ 0 is the Heaviside step function with Θ 0 (0) = 0, the factor (p−2) gives the maximal amount of vertex operators at order p, the factor (p − 2d i − 1) corresponds to the number of relations introduced for every d i < p/2, and the maximum enforces the number to be 0 if there are more than (p − 2) relations in total.
Our description of relations between vertex differential operators exploited the auxiliary statement (2.23) and did not include the Pfaffian vertex differential operators. It is clear, however, that precisely the same reasoning also applies to the latter using eq. (2.24) instead of eq. (2.23). The counting formula (2.27) gets also slightly modified for this Pfaffian case Note that the arguments we have outlined here exhibit relations between vertex operators, but we have not shown that these relations are complete, i.e. that the remaining vertex differential operators are in fact independent. A priori, it could in fact happen that the exceeding relations we get from here, or additional relations obtained from a different reasoning, could provide additional dependencies.
We checked however that summing the counting formulas (2.27) and (2.28) for all allowed orders p in a given dimension d, gives rise to a number of vertex differential operators equal to the number of cross ratios (1.16) associated with every allowed vertex. This provides strong evidence in favor of the independence of our vertex differential operators. In some particularly relevant cases in lower dimensions we can also prove independence.
Example: The N = 6 snowflake channel for d = 7. Let us see how all of this works in the example of a snowflake channel in d = 7, presented in Figure 6. We enumerate the internal links by r = 1, 2, 3. The associated index sets I r,i are I 1,1 = {1, 2}, I 2,1 = {3, 4}, . . . . Here we have two symmetric traceless tensors associated with r = 1, 2. In a more general OPE diagram these could produce a mixed symmetry tensor with maximal depth d = r d = 4, but in the snowflake diagram the field on the third link must also be a symmetric traceless tensor of depth d(I 3 , d = 7) = 2. Our prescription tells us to consider operators (2.21) up to p = 8. Eliminating powers of T 1 = T (12) and T 2 = T (34) higher than 4, it immediately follows that there are no vertex operators of order 8T while there could be up to two operators of order 6T and two operators of order 4 Here and in the following steps we are using notation for which stroked terms are dependent on lower order operators. Let us also recall that the operators with ν = p/2 have been omitted to account for the relation between the vertex and Casimir differential operators for the third leg. The reduction of T 3 = T (56) to a symmetric traceless tensor implies the existence of p − 2d 3 − 1 relations between p order monomials. The only useful relation in this case is the one produced for p = 6, coming from the which can be contracted with either T 1 or T 2 and traced over to get a relation between the sixth order monomials (including the one associated to the Casimir of the third leg): This reduces the amount of independent vertex operators by one, bringing us to a total of three independent operators, which matches with the number of associated cross ratios.

OPE channels and limits of Gaudin models
At this point we have defined a set of differential operators associated with the intermediate fields and the individual vertices of a given OPE diagram with N external fields. The new vertex operators were constructed in Subsection 2.2 from a Gaudin model with three sites, which was crucial in proving their commutativity. Our construction of the vertex operators has been local in its focus on a particular building block, namely a single vertex that is associated to a local element of the OPE diagram.
The purpose of this section is to adopt a more global perspective by showing that the whole set of

N sites Gaudin model and OPE limits
Let us first define the Gaudin model that we will use in this section. Since this construction is similar to the one of the 3-sites Gaudin model considered in Section 2.2, we refer to that section for details and references. As before, we consider a Gaudin model based on the conformal Lie algebra g = so(d + 1, 1) but now with N sites, whose positions w 1 , · · · , w N ∈ C are for the moment arbitrary. We naturally associate these sites with the N external fields in the correlation function under consideration and more precisely attach to each site i ∈ {1, · · · , N } the representation of g defined by the generators T (i) α , which describe the action of the conformal transformations on the scalar field φ i (x i ) in terms of first-order differential operators. Then we define the (components of the) Lax matrix of the model as α are the diagonal conformal generators that also appear in the Ward identities (1.11). This means that Gaudin Hamiltonians descend to correlation functions G.
The In order to make a precise statement, we need a bit of preparation. The limits we are about to discuss must depend on the choice of the OPE channel. So let us assume we are given such a channel C. In order to define the limits we pick an (arbitrary) external edge in the diagram, which will serve as a reference point and which, up to reordering, we can suppose to have label N . As this edge is external, it is attached to a unique vertex, which we will denote by ρ * . Such a choice of reference vertex defines a so-called rooted tree representation of the diagram. We then draw the OPE diagram on a plane, with the vertex ρ * situated at the top and with each vertex having two downward edges attached.
Such a representation on a plane forces us to make a choice of which edges are pointing towards the left and which edges are pointing towards the right: this choice is arbitrary, and gives rise to what is called a plane (or ordered) representation of the underlying rooted tree. We give an example of such a plane rooted tree representation for an 8-point OPE diagram in Figure 7 below. Recall from Section 1 that each vertex ρ of the OPE diagram defines a partition N = I ρ,1 ∪ · I ρ,2 ∪ · I ρ,3 , with the sets I ρ,j formed by the labels of the external fields attached to the three branches of the vertex. Although the choice of labeling of these branches was arbitrary in Section 1, we will now fix it using the plane rooted tree representation of the diagram picked above: choose the branch I ρ,1 to be the one pointing to the bottom left and the branch I ρ,2 to be the one pointing to the bottom right.
By construction, the last branch I ρ,3 then always points to the top and contains the reference point N .
Each vertex ρ in the diagram is thereby associated with a sequence s ρ = (s ρ 1 , s ρ 2 , . . . , s ρ nρ ) of elements s ρ a ∈ {1, 2}. This sequence s ρ encodes the path from ρ * to ρ. It tells us whether we have to move to the left (for s ρ a = 1) or right (for s ρ a = 2) every time we reach a new vertex until we arrive at ρ after n ρ steps. We shall also refer to the length n ρ of the sequence as the depth of the vertex and to s ρ as the binary sequence of ρ. Note that the top vertex ρ * has depth n ρ * = 0. Let us point out that the notion of depth used in this section refers to the distance from the root ρ * and is very different from the depth d introduced in eq. (1.12) of the introduction.
In order to construct the limit of the Gaudin model that we are interested in, we will need to assign a polynomial g ρ ( ) to each vertex ρ. If s ρ is the binary sequence associated with the vertex ρ, the polynomial g ρ is defined as Obviously the top vertex ρ * is assigned to g ρ * ( ) = 0. The vertices of depth n ρ = 1 are associated with g ρ1 ( ) = 0 or g ρ2 ( ) = 1, depending on whether they are reached from ρ * by going down to the left (ρ 1 ) or to the right (ρ 2 ).
Similarly, we can assign polynomials f i ( ) to each external edge 1 ≤ i < N at the bottom of the plane rooted tree. Once again, we can encode the path from ρ * down to the edge i by a binary sequence s i = (s i 1 , s i 2 , . . . , s i ni ). The length n i of the sequence s i is also referred to as the depth of the edge i. Now we introduce To this end, we will first construct the vertex Lax matrices (2.12) from the Lax matrix (3.1) before studying the associated Hamiltonians (2.7) in the limit. As it turns out, we can recover the parameter free Lax matrix L ρ that is associated with the vertex ρ as Let us note that in the limit, the site w N = −1 associated with the reference field N goes to infinity, while the sites of the other external fields approach z = 0 or z = 1 depending on whether they are located at the right or left branch of the the plane rooted tree, i.e. whether their binary sequence s ρ starts with s ρ 1 = 1 or s ρ 1 = 2. We shall prove eq. (3.6) in the third Subsection 3.3 through a recursive procedure that will also offer insight into the construction of the polynomials g ρ and f i .
Let us now turn to the limit construction for the Gaudin Hamiltonians. We claim that the Hamiltonians The fact that this statement holds for the leading part of the Hamiltonians, without quantum corrections, follows directly from the corresponding limit (3.6) of the Lax matrix pnρ κ α1···αp p L α1 nρ z + g ρ ( ) · · · L αp nρ z + g ρ ( ) But it requires a bit of work to argue that the quantum corrections also have the required behaviour under the limit. Consider a term of the form (2.9) in the correction: by appropriately distributing the powers pnρ , using the fact that r 1 + · · · + r q = p, and performing the change of spectral parameter z → z ρ ( ) = nρ z + g ρ ( ) in the derivatives, we find that such that this correction term reduces in the OPE limit to the corresponding correction in H  Hamiltonians can be written as for arbitrary , and is therefore preserved in the limit → 0,

Examples
Before we prove our main result, let us illustrate the construction of the operators from limits of Gaudin models with two examples. The first one addresses the so-called comb channel OPE diagrams for which we have already outlined the limit in [1]. The second example deals with the snowflake OPE channel of the N = 6-point function.
Comb channel. Let us consider the comb channel OPE diagram with N external fields. To apply the construction of the present section, we first need to pick a plane rooted tree representation of this diagram. We will choose to represent it with all internal edges pointing towards the bottom left. We then label the external edges of the tree as follows: we let N be the top edge of the tree, 1 be edge furthest to the left and label by 2, · · · , N − 1 the external edges pointing to the bottom right at each vertex, from the bottom to the top. Moreover, following our conventions in [1], we enumerate the vertices ρ = [r] by an integer r = 1, . . . , N − 2, from bottom to top. We represent this plane rooted tree in Figure 8 below, with external edges indicated in black and vertices in blue.
associated with the vertex [r] then reads In sum, the vertex Gaudin Hamiltonians of the comb channel OPE limit are [r] (z). (3.14) The above limit coincides exactly with the one introduced in [1] to describe the comb channel, thus showing that the results of [1] are contained in the more general construction discussed here.
Snowflake channel. The results of the present article allow us to discuss more general topologies of OPE diagrams than the comb channel. The first example of such a topology is the snowflake channel of 6-point functions. We represent this OPE diagram as a plane rooted tree following the conventions of Figure 9, where the external edges are labeled in black from 1 to 6 and the vertices are labeled in blue from [1] to [4]. Note in particular that the internal vertex of the diagram corresponds here to the label [3]. [3] [4] Figure 9. Choice of plane rooted tree representation of the snowflake OPE diagram.
We can immediately read off the depth of the four vertices as We can now encode the positions of all four vertices in a binary sequence and apply the formulas (3.3) to construct the polynomials g ρ , yielding g [1] ( ) = g [3] ( ) = g [4] ( ) = 0 and g [2] ( ) = . We will now prove that if the limit construction of Subsection 3.1 holds for the subtrees T and T , then it also holds for the initial tree T , thus proving that it holds for any tree by induction. Let us first introduce some useful notation. We will denote by E and E the external edges of T that belong to the subtrees T and T respectively. Note that if T is not trivial, i.e. if e is not an external edge of the initial tree T , then the full set of external edges of T is E ∪ · {e } (since e is external in T but not in T ). On the other hand, if e is external in T , then we simply have E = {e } and the subtree T is trivial. Let us also denote by V , V and V the set of vertices of T , T and T , such Recursion relations. Recall the polynomial g ρ ( ) defined in eq. (3.3) for any vertex ρ ∈ V of T in terms of the binary sequence s ρ = (s ρ 1 , · · · , s ρ nρ ). Let us suppose that this vertex is contained in the subtree T and thus belongs to V : it is then associated with a binary sequence s ρ = (s ρ 1 , · · · , s ρ n ρ ) in T . By construction, the depth n ρ of ρ in T is given by n ρ − 1. Moreover, it is clear that the binary sequence of ρ in T is related to that in T by s ρ = (1, s ρ 1 , · · · , s ρ nρ−1 ). Indeed, since ρ belongs to T , the path from ρ * to ρ starts by going to the bottom left (s ρ 1 = 1), and is then given by the path from e to ρ, encoded by s ρ . It is then clear that the polynomial g ρ ( ) defined by eq. (3.3) is related to the corresponding polynomial g ρ ( ) defined for T by g ρ ( ) = g ρ ( ). Similarly, if ρ belongs to V , we have s ρ = (2, s ρ 1 , · · · , s ρ nρ−1 ) and thus g ρ ( ) = 1 + g ρ ( ). In conclusion, the polynomials g ρ ( ) satisfy the recursion relation (3.20) A similar analysis can be performed for the sites w i = f i ( ) associated with the external edges i ∈ N through eq. (3.4), distinguishing three cases. If i = N is the top reference edge, then we recall that w N = f N ( ) = −1 . If i ∈ E is an edge belonging to the subtree T , then one can relate the binary sequences s i and s i describing i in T and T in a similar way as for the vertices in the paragraph above. We then find that the polynomial f i ( ) satisfy the recursion relation Note that in the first case, the subtree T is trivial and the index i is then necessarily equal to e , while in the second case i is different from e and the subtree T is therefore non-trivial. Finally, if i ∈ E belongs to the subtree T , then we similarly find Here, i = e in the first case and i = e in the second one.
Induction hypotheses. We will now suppose that the limit procedure defined in Subsection 3.1 holds for the subtrees T and T . To phrase these induction hypotheses more precisely, let us focus first on the subtree T . If it is non-trivial, i.e. if e is not external in T , the external edges of T are E ∪ · {e }. We then introduce the Gaudin Lax matrix associated with T as As T is assumed to be non-trivial here, its vertex set V is non-empty. The induction hypothesis that we make in this subsection is then that eq. (3.6) holds for the subtree T , that is to say In this equation, we used the Lax matrix L ρ α (z) associated with the vertex ρ as defined in the initial tree T : indeed, it is clear that this vertex Lax matrix coincides with the one associated with ρ in the subtree T (in particular, the subsets of external edges I ρ,1 and I ρ,2 associated with the left and right branches of ρ are the same when defined for T as when defined for T ).
Similarly, if T is non-trivial, we consider the associated Lax matrix , (3.25) and suppose that it satisfies the induction hypothesis Proof of the induction. We are now in a position to prove that the induction carries from the subtrees T and T to T . For that, we will show that the limit (3.6) holds for every vertex ρ ∈ V , with three cases to distinguish. If ρ ∈ V belongs to the subtree T , we will use the induction hypothesis Finally, if ρ is the reference vertex ρ * , then the limit will follow without having to use any induction hypothesis. As these proofs are rather technical, we gather them in Appendix A.

Example: 5-point conformal blocks
As an example of our construction of commuting differential operators, let us consider a correlator of five scalar fields and fix the OPE decomposition as in Figure 2.
This correlator can be be written schematically as where Ω (∆i) 5 (x i ) is a prefactor that takes into account the covariance of the correlator with respect to conformal transformations, while ψ (∆i) (u 1 , . . . , u 5 ) is a conformally invariant function which depends on five cross ratios and admits a conformal block decomposition. In order to obtain differential equations for 5-point conformal blocks, one first needs to determine which Casimir and vertex operators characterize these blocks, and then compute their action on the space of cross ratios u i .
For the OPE decomposition of Figure 2, the recipe of Section 2 instructs us to construct four Casimir operators, two for each internal leg and one vertex operator Note that -in agreement with the general recipe -the vertex operator is not uniquely defined and (4.7) can be shifted by terms proportional to (T 1 + T 2 ) 2 (T 3 ) 2 or to the Casimir operators. 3 To make explicit computations, we will use the embedding space formalism of [25], with which one can efficiently compute the action of the differential operators in cross ratio space; we briefly review this formalism in Appendix B. Note here that the dimension of spacetime only appears in our computations as a parameter when contracting Kronecker deltas δ A A = d + 2, and can be therefore kept generic. The first step is the choice of prefactor Ω (∆i) 5 and cross ratios u i . We chose to use the same conventions as [28], where the author computed 5-point blocks in the case of scalar exchange; the expression for Ω (∆i) 5 in physical space coordinates can be easily translated into one in embedding space through the simple relation and one obtains (up to an overall normalization) the prefactor (4.9) Regarding the cross ratios, it is natural to build four of these using the same construction of the standard 4-point cross ratios (u, v) introduced in [7], but supported on two different sets of points while an interesting choice for the fifth cross ratio is In comparison with a potentially more natural parameterization using five independent 4-point cross ratios, as in e.g. [31,43], this parameterization of cross ratio space has the advantage of presenting all of our differential operators with polynomial coefficients in the u i .
Using the scalar representation for generators in the embedding space (B.3), the operators (4.3)-(4.7) can be easily expressed as objects D (Xi) acting on the coordinates X A i . To obtain their action on the space of cross ratios D (ui) , one simply conjugates the D (Xi) by the prefactor as follows F (u 1 , . . . , u 5 ) . and the u i 's in terms of scalar products. The LHS is then obtained by solving (4.10-4.11) for five different scalar products and substituting them in the RHS after the conjugation has been done; the remaining scalar products will drop out and the final answer for the LHS will be expressed only in terms of the cross ratios.
To implement this procedure more concretely, we attach to this publication a Mathematica notebook 4 , where we present the explicit computation of the quadratic Casimirs, as well as the final expressions one gets for the fourth-order Casimirs and the vertex operator (4.7) by trivially extending the same algorithm.
As an attempt to simplify the analytic expressions for the differential equations, it is natural to try to extend the 4-point change of coordinates of Dolan and Osborn [8,9] to the 5-point case. A very good candidate for this purpose is the change of coordinates which leads to the simplest expressions for the quadratic Casimirs that we could find. Indeed, by 1 .

(4.22)
We have also attempted similar types of factorizations for the quartic Casimirs and the vertex operator, in the spirit of the decomposition in equation (4.14) of [9]; so far to no great avail.

Conclusions and Outlook
In this paper we have initiated the construction of multipoint conformal blocks for correlation functions of N scalar fields. More concretely, we constructed a number of independent commuting differential operators on N -point functions that matches the number of conformally invariant cross ratios. In consequence, these operators uniquely determine conformal blocks as bases of joint eigenfunctions into which one can expand any correlation function of N scalar fields. The set of commuting differential operators depends on the choice of an OPE channel. We constructed them for all channel topologies, extending the results we announced in [1]. Our construction relies on the identification of these operators with the Hamiltonians of certain OPE limits of integrable Gaudin models, see Section 3.
Let us note that, in the special case of d = 2, the relation between conformal blocks and Gaudin models is consistent with recent observations in [58] concerning ambitwistor realizations of scattering equations in AdS 3 , given the close connection between conformal blocks and scattering amplitudes in anti-de Sitter space [59]. But it is important to stress that the N -site Gaudin Hamiltonians simplify significantly in the OPE limits relevant to the theory of conformal blocks. While the 4-site Gaudin system for d = 1 is related to the Heun equation, for example, one obtains hypergeometric differential equations in the OPE limits.
In any given channel, the set of differential operators decomposes into two subsets formed by Casimir differential operators on the one hand, and the novel vertex differential operators on the other hand.
Individual vertices can give rise to multivariable integrable systems with up to n v (d) independent commuting Hamiltonians in the most generic case. All such Hamiltonians are of order higher than two and they have not received much attention in the physics literature so far. In order to develop a theory of these vertex systems it seems reasonable to start with the single-variable cases. According to formula (1.16), all vertices that appear in the comb channel of an N -point function in d = 3, 4 contribute one single degree of freedom. The same is true for the nontrivial vertices in the comb channel of N = 5, 6-point functions in any d > 2. Therefore, single-variable vertex systems already cover many cases of interest. The associated vertex differential operator has already been constructed in [1], though not for the most general case. In the sequel to this work we will revisit this operator, construct it for all single-variable vertex systems, and show that it can be mapped to the fourth order Hamiltonian of an elliptic Calogero-Moser-Sutherland model that was first discovered about 10 years ago by Etingof, Felder, Ma and Veselov in [42].
There are two important directions that we will address in order to develop the theory of multi-point blocks. The first one concerns the extension of what we described in the previous paragraph to vertices with more variables, such as the central vertex in the snowflake channel of a 6-point function in d > 2.
We will discuss this in a forthcoming paper. The second important direction concerns the study of (joint) eigenfunctions for the vertex differential operator(s). Unlike the single variable of a 4-point function in d = 1, for which the Casimir differential equation is the ordinary hypergeometric equation, even the single variable case of the vertex system seems to be uncharted, let alone its multi-variable extensions. We plan to address this challenge in the future before we combine the knowledge about vertex systems back into a full theory of multi-point conformal blocks. In this context, it would also be interesting to examine the relations of our approach with the weight shifting technology of multi-point blocks with non-scalar exchange recently exhibited in [41].
While executing the entire program still requires a significant amount of work, applications can be pursued in parallel. This applies in particular to the multi-point bootstrap that was also advocated recently in [43]. One can expect the study of multi-point correlation functions to give access to new dynamical information that is not visible in scalar 4-point functions. The analytic bootstrap study of 4-point functions has provided powerful access to anomalous dimensions of double twist operators, as well as OPE coefficients for the operator product of two fundamental scalars into double twist operators [60,61]. Similarly, crossing symmetry for N = 6-point functions constrains the anomalous dimensions of triple twist operators and various OPEs of the fundamental field with double twist operators, depending on channel topologies that are analyzed.

A Proof of the induction in the limits of Gaudin models
In this appendix, we detail the induction in the proof of the limit procedure of Section 3. We thus refer to this section for notation and definitions. Our main goal is to show that for every vertex ρ ∈ V , the Lax matrix of the N -sites Gaudin model satisfies the limit (3.6). For the purposes of this appendix, it will be useful to rewrite this limit as In the left-hand side of (A.1), and in all this appendix, we fix the sites w i of the Gaudin model to their value w i = f i ( ) prescribed by the limit procedure of Subsection 3.1. Recall that ρ is either the reference vertex ρ * , in V or in V . We will treat these three cases separately.

A.1 Reference vertex
Let us first consider the reference vertex ρ * . The definition of n ρ * and g ρ * made in Subsection 3.1 can be simply rewritten as n ρ * = 0 and h ρ * (z, ) = z.
To conclude, we observe that the labeling of branches at vertices made in Subsection 3.1 using the plane rooted tree representation T of the diagram implies that E = I ρ * ,1 and E = I ρ * ,2 . This shows that the limit (A.1) is satisfied for the reference vertex ρ * .

A.2 Vertices in V .
Let us now consider the case of a vertex ρ ∈ V in the subtree attached to e . Note that the existence of this vertex requires the subtree T to be non-trivial and thus the edge e to be intermediate in the initial diagram. The definition of n ρ and the recursion relation for g ρ ( ) in the first line of eq. (3.20) can be rewritten as where we have defined h ρ (z, ) = n ρ z + g ρ ( ) as the equivalent of h ρ (z, ) for the subtree T . We thus have We will once again separate this expression in three parts, coming from the contributions of the reference edge N and the two subtrees.
Contribution of the reference edge. Let us start with the edge N , with associated site w N = f N ( ) = −1 . Its contribution to (A.13) in the limit → 0 is then given by Contribution of the subtree T . Let us now consider the contribution coming from the subtree T attached to e . If e is external in the initial tree T , then we simply have w e = 1 according to the first line of eq. (3.22) and the contribution to (A.13) is If e is initially intermediate, then the contribution comes from the external edges i ∈ E ⊂ N whose The Gaudin Lax matrix of the subtree T (the one of the full tree, not the ones associated with vertices) is given by eq. (3.23) and thus can be rewritten as where we used the fact that f e ( ) = −1 since e is the reference edge of T . Thus, we can rewrite the above contribution (A.17) as i∈E n ρ +1 T The second term vanishes in the limit → 0. Moreover, by the induction hypothesis (3.24) for the subtree T , the first term tends to L ρ α (z) in this limit. In conclusion, we thus get i∈E n ρ +1 T when → 0, we thus get in the end that as required. This indeed shows that the limit (A.1) is satisfied for vertices ρ ∈ V , using the induction hypothesis (3.24) that a similar property also holds in the subtree T .

A.3 Vertices in V
Let us finally consider a vertex ρ ∈ V . The procedure here will resemble the one in the previous subsection so we will not describe it in detail. The definition of n ρ and the recursion relation for g ρ (  B Embedding space formalism and index-free notation In this section we will briefly review the embedding space formalism of [25], as well as a slight generalization of their index-free notation that we are going to use in the proof of Appendix C. We will in particular focus on scalar and tensorial representations. In short, the embedding space formalism makes use of the isomorphism between the conformal group in d dimensions and the Lorentz group in d + 2 dimensions -both corresponding to SO(d + 1, 1) -to map objects from one space to the other. Points x ∈ R d are mapped to (d + 2)-dimensional vectors X A = (X + , X − , X µ ), here presented in light-cone coordinates with metric ds 2 = −dX − dX + + δ µν dX µ dX ν , and to make the degrees of freedom match, these vectors are constrained to the projective light-cone With this mapping, scalar fields become homogeneous objects of degree −∆ with respect to the coordinates X A they depend on This has the big advantage that the conformal generators, some of which are non-linear, once mapped to the embedding space, acquire all the same linear form, namely the one of Lorentz generators: where X A is obtained by lowering the index of X A by contraction with the flat metric η AB . This not only simplifies the analytic expressions that one works with, but also allows one to carry out computer algebra computations in a much faster and general way through the implementation of indices; this proves very useful for the computation of operators for scalar multipoint blocks like those of Section 4.
When one considers fields in spinning representations, labeled by L spins where indices in the same bracket are symmetrized over, one also requires them to be transverse with respect to any index and traceless with respect to any pair of indices Expressions involving these tensorial representations can be further simplified by the use of an indexfree notation in embedding space. One can implement this by either making the antisymmetries of the tensors manifest with fermionic variables, as was first done in [62], or by making the symmetries manifest with bosonic variables, in a form like Appendix B of [63]. In agreement with our choice of representation for (B.4), we will keep the symmetries manifest and use the latter formalism. The core idea is to remove indices from the tensorial expressions (B.4) by converting them to polynomials on m auxiliary vectors Z i ∈ C d+2 through the contractions To keep everything consistent with the conformal covariance, tracelessness and transversality of the tensors, the coordinate vectors need to satisfy The action of the conformal generators on these fields is also modified to take into account the Z i dependence, so (B.3) becomes (B.11) Using this type of representation for generators and fields, it becomes much simpler to obtain Casimir and vertex differential operators for spinning external fields.

B.1 Classical Embedding space
For the proof in Appendix C we are going to need a classical version of the embedding space formalism for the spinning fields that we have just discussed. For this purpose, we introduce P A as the conjugate momentum to X A , and Q i A as conjugate momenta of the auxiliary variables Z A i . To see which constraints are imposed on these classical variables, one can check the action of the operators on fields or correlation functions to get conditions that need to be satisfied by scalar products of coordinates and momenta. Some of the operators in (B.12) are in fact evaluated to constants on correlation functions. By imposing the same behavior when replacing derivatives with momenta, one 6 this is a "partial" flag manifold if L + 1 < r d , and a "full" flag manifold if L + 1 = r d .
obtains the following relations for phase space variables: (B.13) One can directly verify that the conditions (B.13) combined with the classical version of the generators (B.11)T lead to the correct classical Casimirs associated to mixed-symmetry tensors:

C Relations among vertex differential operators
Our goal in this appendix is to justify the relations (2.23) and (2.24) for the total symbols of our vertex differential operators. There are many ways to derive these relations. Here we shall follow a more pedestrian approach that does not require much background from representation theory.
In order to derive the relations (2.23), (2.24) we first note that these were formulated in terms of the coordinates and momenta of the external scalar fields. The representation of the conformal algebra that is associated to the index set I decomposes into an infinite number of spinning representations.
Each of the irreducible components can be prepared in embedding space formalism (see appendix B). Here we shall study the relations (2.23) for a given irreducible component so that the coefficients f,s are functions of the associated weight and spins rather than functions of symbols of the Casimir differential operators.
After these introductory comments let us approach relation (2.23) by considering the simplest example in which the intermediate irreducible representation is scalar. We can construct explicitly the first three matrices (T n ) A B 7 in the classical embedding space that we introduced in Appendix B.1: (C.1) 7 Here and in the discussion below we shall drop the subscript f .
It is then clear that, when considering scalar representations, the powers of generatorsT n with n ≥ 3 will be dependent on lower powers of the generators. This directly implies that any vertex operator of the type (2.21) that contains a power of a scalar generator higher than two will become dependent on operators of lower order, e.g.T 3 · B ∝T · B . (C.2) We would now like to prove that something analogous to (C.1) is valid for representations of higher depth d and higher powers, respectively. Let us then consider the generators for a mixed symmetry tensor of depth d,T we expect the n-fold contractions of generatorsT n to be independent up to power n = 2d, with the first dependent object produced at power n = 2d + 1. To prove this, let us focus on a specific matrix entry (T n ) A B of the powersT n and construct the following submatrix of the Jacobian for fixed indices A, B and C: If what we argued above holds, we should be able to see that some 2d-minors of (C.4) are equal to zero, while the same matrix with the last row dropped out has all nonzero 2d-minors. We checked this with Mathematica symbolically for the case of symmetric traceless tensors, and numerically for mixed symmetry tensors of depth d ≤ 6, exhausting all tensorial representations that are allowed in the range of dimensions of known CFTs. One can use a similar reasoning to show eq. (2.24).
Let us finally comment on a more conceptual interpretation of the results of this appendix. We consider first the case of a scalar representation, whose generators can be gathered in the matrix (T f ) A B = X A P B − X B P A in the fundamental representation. Using the relations X A X A = 0 and X A P A = −∆, one can show that the matrixT f is diagonalisable, with eigenvalues ∆, −∆ and 0 (with multiplicity d). It is a standard result of linear algebra thatT f is then annihilated by the polynomial with simple roots equal to these eigenvalues, namely T (T − ∆)(T + ∆) = T 3 − ∆ 2 T : we recover this way thatT 3 f = ∆ 2T f . A similar argument can be formulated for a representation with higher depth d = L + 1, characterized by a weight ∆ and L spins l 1 , . . . , l L . The (symbols of the) generators of this representation are given by eq. (B.14) and can be gathered in a matrixT f valued in the fundamental representation.
The traces of odd powers ofT f vanish, while the traces of even powers are given by the classical Casimirs (B.15). These traces are the Newton sums d+2 i=1 λ p i of the eigenvalues λ 1 , · · · , λ d+2 ofT f and thus determine these eigenvalues uniquely (up to permutation). More precisely, we find thatT f has eigenvalues ∆, −∆, l 1 , −l 1 , . . . , l L , −l L and 0 (with multiplicity d − 2L). If we suppose thatT f is diagonalizable, it is then annihilated by the polynomial with simple roots equal to these eigenvalues, This shows that the powerT 2L+3