Thermal Correlation Functions of KdV Charges in 2D CFT

Two dimensional CFTs have an infinite set of commuting conserved charges, known as the quantum KdV charges, built out of the stress tensor. We compute the thermal correlation functions of the these KdV charges on a circle. We show that these correlation functions are given by quasi-modular differential operators acting on the torus partition function. We determine their modular transformation properties, give explicit expressions in a number of cases, and give a general form which determines an arbitrary correlation function up to a finite number of functions of the central charge. We show that these modular differential operators annihilate the characters of the (2m+1,2) family of non-unitary minimal models. We also show that the distribution of KdV charges becomes sharply peaked at large level.

We will therefore just write equations for the left-moving sector explicitly, remembering that the corresponding right-moving results can be obtained by complex conjugation. The trace in equation (1.1) can be taken over the entire Hilbert space of the CFT on a circle, or over a particular Virasoro module; we will give results in both cases. The dependence of the partition function (1.1) on the chemical potentials encodes the consequences of Virasoro symmetry for finite-temperature physics in a useful way. The purpose of this paper is to explore the structure of this ensemble, at finite temperature and finite central charge.
The KdV charges are obtained by integrating a conserved current built from positive powers of the stress tensor around a spatial circle. The first three charges are the zero modes of the operators J 2 ≡ T , J 4 ≡ (T T ) and J 6 ≡ (T (T T )) + (c+2) 12 (T T ), where T ≡ T zz is the left-moving stress tensor, the round brackets denote conformal normal ordering, and the prime is a spatial derivative. 2 The additional term in J 6 is required to ensure that the zero modes commute. There will be similar terms in all the higher spin operators. The zero modes of these operators can be written explicitly in terms of the Virasoro modes of the stress tensor, as L −n L n + L 2 0 − (2k + 1 6 )L 0 + k k + 11 60 , : L n 1 L n 2 L n 3 : + We have defined k = c/24 for convenience. With these explicit expressions one can verify (with some work) that the I 2m−1 are indeed mutually commuting values, so we will use the notation X = T r Xq L 0 −k for the thermal average. We also use Z to denote the thermal partition function Z [β] with no further chemical potentials. We study the calculation of these thermal expectation values and their modular transformation properties.
One motivation for understanding the GGE is that it is important in the application of the eigenstate thermalization hypothesis (ETH) [15][16][17][18][19] to two-dimensional conformal field theories (CFTs), and more broadly in the understanding of chaos in two-dimensional CFTs. 4 Two dimensional CFTs present an especially interesting case as they have an infinite-dimensional symmetry algebra which does not (except in certain special cases, such as the minimal models) trivialize the dynamics. Thus 2D CFTs can exhibit chaotic dynamics [30], and hence could be expected to obey the ETH. The implications of our results for ETH will be discussed in a companion paper [1]. A related discussion will appear in [31].
We will obtain exact results for the thermal expectation values of the KdV charges at finite temperature and finite central charge. We will first discuss the general features of these expectation values and see how these features allow us to determine the form of some of these thermal expectation values without any explicit calculation. We will then discuss the explicit calculation using two different methods. The explicit calculation can in principle be carried out for any KdV charge, but in practice the calculation quickly becomes cumbersome. We will describe certain techniques (based on the structure of Virasoro null states) which give exact results when a naive-brute force computation is not feasible.
The key features of these thermal expectation values are: • The n-point function I 2m 1 −1 . . . I 2mn−1 in a given Virasoro module is a differential operator acting on the character of that module, regarded as a function of the torus modular parameter (i.e. of temperature). The form of this differential operator depends only on the central charge. As a consequence, the thermal expectation values in the full CFT are given by the same differential operators acting on the ordinary thermal partition function Z[β]. This feature expresses the kinematic nature of the KdV charges; the values of KdV charges on descendent states are entirely determined by the structure of the Virasoro algebra.
• The one-point functions I 2m−1 are modular forms of weight 2m. Thus the differential operator is a modular differential operator. It can be written in terms of the modular covariant Serre derivative, with coefficients which are themselves modular forms (and hence linear combinations of the Eisenstein series E 4 and E 6 ). 5 • The n-point functions I 2m 1 −1 . . . I 2mn−1 are quasi-modular forms of weight 2 i m i and depth n − 1. Thus the differential operator is a quasi-modular differential operator which is a linear combination of modular differential operators with coefficients containing up to n − 1 powers of the Eisenstein series E 2 .
We establish this using general results of [32] on the relation between line integrals and surface integrals of operators on the torus. When applied to n−point functions this leads to a certain "anomaly" under the modular S-transformation, which leads to correlation functions which are quasi-modular rather than modular.
• I r vanishes as an operator in the (s+2, 2) minimal model when r is a multiple of s [4,33]. This means that the differential operators for I r . . . , where . . . denotes any combination of KdV charges, will annihilate the characters of this minimal model. In particular, the differential operator for I 2m−1 is an m th order differential operator which annihilates the m characters of the (2m+1, 2) minimal model when the central charge is set to the appropriate value. This is a remarkable relation, but it only constrains expectation values for a particular value of the central charge, so will not be particularly useful in computations at general central charge. We will therefore postpone the discussion of its proof to the end of the paper.
Before diving into details, let us illustrate these statements in the simplest case. For the first non-trivial KdV charge I 3 , we find where we introduce the notation ∂ = q∂ q = −∂ β . In the second form of the equation we have used the Serre derivative, defined as where E 2 is an Eisenstein series. These derivatives are useful because they map modular forms to modular forms: D r applied to a modular form of weight r gives a modular form of weight r + 2. Thus, since the partition function is modular invariant, the one point function I 3 will transform as a modular form of weight 4. Similarly, we find The statement that I 2 3 is a quasi-modular form of weight 8 and depth 1 reflects the fact that the modular differential operators in the square brackets have weight 8 and 6, respectively, and that only one power of E 2 appears in this expression. The appearance of E 2 in this expression reflects the existence of an "anomaly" under modular S-transformations. Indeed, the coefficient of E 2 in this expression is precisely the thermal one-point function of an operator appearing in OPE of J 4 with itself. This will be described in more detail in section 2. 6 The above modular differential operators can be applied to the partition function of a CFT or to an individual Virasoro character, depending on whether one wants to obtain these expectation values in a full CFT or in a particular representation.
The relationship with minimal models is found by evaluating these differential operators in the (5, 2) minimal model, i.e. the Yang-Lee model with k = − 11 60 . One can indeed check that these modular differential operators annihilate the characters of the Yang-Lee model. Indeed, at k = − 11 60 the differential operator (1.4) is precisely the one associated with the null state in the vacuum representation of the Yang-Lee model. 7 In section 2 we will give the general argument that the n-point functions are modular derivatives of the partition function, involving up to n − 1 powers of E 2 . This reduces the problem of determining a particular thermal expectation value to fixing a finite number of undetermined coefficients.
In section 3 we show how some of these coefficients can be fixed from the action of the KdV charges at low levels, and by the structure of null states in Virasoro representations. This allows us to determine a number of the thermal expectation values exactly. The complete list of results so obtained is summarized in appendix C. Although the resulting differential operators appear unfamiliar, in fact they can be written in terms of more familiar hypergeometric differential operators, as explained in appendix B.
Sections 4 and 5 describe two other, more straightforward approaches to the computation that can be used to verify the results of section 3. In these sections we carry out the explicit calculations of the thermal expectation values directly from the definition. This allows us to explicitly confirm the expected structure, and also provides an approach which can in principle be pushed to arbitrary order, although in practice the calculations become laborious and far less efficient than the methods described in section 3. We carry out explicit calculations in two ways: • We can write the KdV charge in terms of Virasoro modes, and cycle the modes around the commutator to relate the thermal expectation value to derivatives of the ordinary thermal partition function of the CFT. We will illustrate this in detail for the computation of I 3 in section 4. This approach is the most straightforward; however, dealing 6 A similar structure has actually appeared in the literature before in the context of W-algebras. The authors of [34] computed the modular transformation properties of Tr W 2 0 q L−0−k , where W0 is the zero mode of the spin 3 generator W (z) of the W3 algebra. They discovered that this trace transformed like a modular form, except for an anomaly term that is proportional to E2 times the differential operator appearing in equation (1.4). In light of the discussion in section 2 this is easy to understand: the operator J4, whose zero mode is I3, appears as the coefficient of the 1/z 2 term in the W (z)W (0) OPE. 7 We refer to [35] for a more detailed discussion of the construction of these modular differential operators.
It would be interesting to investigate whether our results are relevant for the approach to the classification of CFTs based on the classification of modular differential operators, advocated in [36] and used recently in [37,38].
with all the terms rapidly becomes unwieldy as we move to higher powers or charges of higher degree. Some computational details are relegated to appendix E.
• We can use Zhu's recursion relations [39], as extended in [40] (see also [35]). These recursion relations give us expressions for the thermal n-point functions of the stress tensor, which we can then appropriately integrate to get thermal expectation values of the KdV charges. The recursion relations are based on the same kind of cycling of operators around the trace, but they provide a powerful general encoding of these results, which make some of the modular properties of the result more manifest. Using this technique allows us to push the calculation to higher order. We discuss the calculation in section 5, though we relegate some of the details to appendices F, G.
In section 6, we discuss the relation to the minimal models, showing that the KdV charges I r vanish as an operator in the (s + 2, 2) minimal model whenever r is a multiple of s. This makes explicit certain properties of the KdV charges which were asserted in [4,33]. The relation fixes the coefficients in the n-point functions for the particular values of the central charge corresponding to the minimal models. It is particularly useful for the one-point functions, which get related to the differential operators associated to null states [35,41].
In the final section, we will discuss an important application of these results: to determine the statistics of KdV charges. In particular, by explicit evaluation of the differential operators on a Virasoro character we can determine the moments of the KdV charges at arbitrary level. For example, by applying equations (1.4) and (1.6) to the Verma module character we can determine the mean and the variance of the I 3 at each level in a generic Virasoro representation. The result of this is that at high level the distribution of KdV charges will become very sharply peaked around its average value. We will make a few remarks on this in section 7, and expand further on the statistics of KdV charges in [1].
In summary, we have determined the general structure of the thermal expectation values of n-point correlation function of the KdV charges, and found explicit expressions for a number of cases. These expressions are exact results at finite temperature and finite central charge. This is a useful step towards understanding the GGE for two-dimensional CFTs; the structure is however too complicated to allow for a simple exponentiation to obtain an explicit formula for the GGE partition function at finite c and finite temperature.
A key element in our derivation is the determination of the modular transformation properties of the n-point functions. In our companion paper [1], we will apply these results to explore the high-temperature behaviour of the GGE at finite central charge and consider the comparison to expectation values of KdV charges in a typical high-energy microstate.
Many open problems remain. These include finding more efficient methods for determining the remaining free parameters in the thermal expectation values, and applying similar techniques to theories with extended symmetry algebras.

Modular differential form of the thermal expectation values
In this section we explain how the thermal expectation values of KdV charges can be written as modular derivatives of the partition function. This fixes the form of the thermal expectation values up to a finite number of coefficients, which are polynomials in the central charge.

Differential operator form
For ordinary conserved charges, such as a U (1) charge, the partition function with a chemical potential for the charge carries additional information about the theory above and beyond what is available in the thermal partition function Z(β). That this is not the case for the KdV charges should not be surprising; the KdV charges evaluated on primary states are simply functions of the L 0 eigenvalue of the primary, and the KdV charges evaluated on descendents can in principle be determined by applying the Virasoro commutation relations One consequence is that the expectation values of KdV charges in a Virasoro module can be written as derivatives of the character with respect to the modular parameter (and hence thermal expectation values in the full CFT are the same derivative of the partition function).
To understand this a bit more explicitly, consider the thermal expectation value of a product of KdV charges evaluated in a Virasoro module built on a primary state of dimension h: Here the trace runs over the Virasoro descendents of a primary state with L 0 = h. Using the explicit expressions for the KdV charges, the RHS can be rewritten as a polynomial in the Virasoro modes. Virasoro modes L r for r = 0 can be cycled around the trace using the commutation relations. We pick up a factor of q −r when L r commutes past q L 0 , and commuting it past the other modes gives terms which are lower order polynomials in the Virasoro modes. Thus, we can iteratively rewrite the expression in terms of functions of q times traces of lower polynomials in the Virasoro modes, leaving us ultimately with an expression involving functions of q times traces of polynomials in L 0 . The factors of L 0 can then be rewritten in terms of derivatives of the character χ h (q) = Tr h [q L 0 −k ] with respect to q. Since I 2m−1 involves at most m factors of L 0 , this expression contains derivatives up to order M = i m i , and our expectation value is an M th order differential operator where we use the notation ∂ = q∂ q . The coefficients are functions of q and k (the dependence on k coming from the central term in the Virasoro commutator (2.1)) but do not depend on the conformal dimension h of the character: the differential operator is universal, depending only on the choice of KdV charges. We can obtain the thermal expectation value in the full CFT by summing over primaries, giving us an expression where the same differential operator acts on the partition function: In section 4, we will discuss the explicit evaluation of the thermal expectation values using this procedure, but for now, the key point is that it determines a general form for the expectation values, as a differential operator acting on the character or on the partition function. 8

Modular transformation properties
Since the expression obtained in the previous subsection is a function of the modular parameter, it is natural to ask about its behaviour under modular transformations.
Let us begin with the expectation value of a single KdV charge, where the answer is simple. The charge is the contour integral of a conserved chiral current J 2m around the spatial circle: If we describe the torus by the standard z coordinate, with z ∼ z + 1 ∼ z + τ , then J = 1 0 dz 2π J(z). The integration around the spatial circle appears to break modular symmetry, since we have picked a particular cycle of the torus around which to integrate. However, since the current J 2m is conserved we can freely translate the spatial integral around the torus. This means that we can average over the location of the spatial circle to turn the one-dimensional integral over the spatial circle into an integral over the entire torus. The result is that where, following [32], we have introduced the following notation for integrals over the torus: The factor of 1/τ 2 arises because when we average over the torus, we must divide by the area of the torus. An important point is that with this normalization the integral 2πτ 2 is itself modular invariant. The final step in the argument is to note that the operator J 2m is (up to total derivative terms) a quasiprimary with dimension 2m. Thus under a modular transformation its one point function J 2m (z) transforms like a modular form of weight 2m.
The result of equation (2.6) is therefore that I 2m−1 transforms like a modular form of weight 2m.
We saw in the previous section that the expectation value is a combination of derivatives of the partition function. It is now natural to rewrite these derivatives in terms of the Serre modular derivative. To simplify the notation, we will denote the Serre derivative of a modular form A of weight r by DA = D r A, where D acts as a derivation on products of modular forms as D(AB) = (DA)B + A(DB).
The statement that I 2m−1 is a modular form of weight 2m implies that we can write it as a combination of the derivatives D k Z where k = 0, . . . m, with coefficients which are modular forms of weight 2m − 2k. We can then write these modular forms in terms of the Eisenstein series E 4 and E 6 to obtain a general expression for the one point function: where the coefficients c i are functions only of the central charge k of the CFT. We will see that they are polynomials in k of increasing order. 9 The modular properties thus suffice to determine the q-dependence in the differential operator, determining the operator up to a finite number of coefficients.
For the higher-point functions, there is a subtlety which implies that the differential operators are not exactly modular covariant. For example, consider the two-point function We wish to convert these contour integrals into integrals over the torus by using the fact that the currents J 2n and J 2m are conserved. But averaging the two operators over the torus will give a singular contribution where the operators coincide, which is absent in the expression with a pair of contour integrals. This contribution was analyzed carefully by Dijkgraaf [32], who showed that where [J 2m J 2n ] 2 denotes the coefficient of (z − w) −2 in the OPE J 2m (z)J 2n (w). 10 Again, rewriting the contour integrals in terms of torus integrals allows us to understand the modular transformation properties. The correlation functions on the RHS are modular forms of weight 2m + 2n and 2m + 2n − 2 respectively, but the explicit factor of τ 2 in the second term makes it transform inhomogeneously under a modular transformation. The result is that the two-point function I 2m−1 I 2n−1 is a quasi-modular form. 9 Note that for sufficiently high weight, there is more than one modular form of a given weight; for example E 3 4 and E 2 6 are both modular forms of weight 12. So at higher orders there will generically be more than one term in the sum at a given order in derivatives. 10 We use square brackets here to avoid confusion with the round bracket notation for conformal normal ordering.
When we express this correlation function in terms of a differential operator, this inhomogeneous term is reproduced by allowing a single factor of the Eisenstein series E 2 to appear, as E 2 has an inhomogeneous transformation under modular transformations which reproduces the transformation of (2.10). Hence the expectation value is a combination of a modular form of weight 2m + 2n, and E 2 times a modular form of weight 2m + 2n − 2. The general form for the differential operator is (2.11) where again the coefficients c i , d i are functions only of k. The first term in parentheses is a modular form of weight 2m + 2n. The second term in parentheses is a modular form of weight 2m + 2n − 2, and is the one point function of the operator [J 2m J 2n ] 2 on the torus.
For higher-point functions of the KdV charges, more coincident singularities will appear when we convert contour integrals to torus integrals. The approach of [32] can be applied to express these in terms of inhomogeneous terms involving the second-order OPE singularity of the coinciding operators. The result is that the n-point correlation function of KdV charges I 2m 1 −1 . . . I 2mn−1 is a quasi-modular form of weight 2 i m i and depth n − 1, meaning that the differential operator will contain terms with up to n − 1 powers of E 2 , with each term having total weight 2 i m i .
These modular transformations will be explored in more detail in specific cases in [1], where we use them to relate the high-temperature behaviour of the correlation functions to their low-temperature limit. Here the key point is that they strongly constrain the form of the differential operator appearing in the KdV n-point functions.

Determining thermal expectation values: q expansion and null states
In the previous section we derived a general form for the thermal expectation value of a product of KdV charges in terms of a differential operator which is determined up to a finite number of constants. In this section, we discuss to what extent we can use data from low levels in a Virasoro module, and the structure of Virasoro representations, to determine these coefficients. We will discuss only a few representative computations. A complete list of results obtained using this method is given in Appendix C.

One-point functions
We consider first the one-point functions I 2m−1 h . The character is χ h (q) = q h−k (1 + q + . . .) for h = 0; the action of the differential operator on this character will then have a q-expansion where I L=0 2m−1 is the value of the KdV charge on the primary, and I L=1 2m−1 is its value on the unique state at level one. 11 It is easy to determine the value of I L=0 2m−1 from the basic expression for the KdV charges (1.2); it is simply a function of the L 0 eigenvalue of the primary, h. It will be a polynomial of order m in h, as the KdV charge I 2m−1 involves up to m powers of L 0 . We can then compare this to the result that one would obtain if one acts on the character χ h (q) with a general differential operator of the form (2.8), and use this to fix the undetermined coefficients in the differential operator. When this differential operator acts on the character, the m powers of h come from the m derivatives with respect to q acting on the q h−k term in the character. Thus, knowledge of I L=0 2m−1 will fix one coefficient at each order in derivatives. This will fix the first few terms in the differential operator uniquely. For sufficiently small values of m this is enough to determine the differential operator uniquely.
The values of KdV charges I 3 up to I 15 on primary states were tabulated in appendix B of [33]. We can use this to obtain, for example,

2)
One can check that the values in minimal models vanish as advertised. In the Yang-Lee (5, 2) minimal model, k = − 11 60 , and there are two primaries: the vacuum with h = 0 and a scalar of weight h = − 1 5 . The vacuum representation has a null vector at level 4. In [35], this null state was shown to imply a differential equation for the characters of this theory, This is precisely our differential operator I 3 for k = − 11 60 , so we see I 3 h = 0 in the Yang-Lee model. Similarly the (7, 2) minimal model has k = − 17 42 . The differential operator in (3.3) for this value of k reduces to the differential operator which annihilates the characters of the (7, 2) minimal model [41].
This method allows us to determine the differential operators for I 3 up to I 9 . However, for I 11 and higher, we are confronted with another problem. The terms in the differential operator are modular forms, and for weight w ≥ 12 there is no longer a unique modular form of weight w. Thus knowledge of I L=0 2m−1 alone is no longer sufficient to determine the modular differential operator uniquely. Fortunately, we can use the structure of null states in Virasoro representations to obtain further constraints. For example, we know that when h = 0 the representation has a null state at level 1. Thus when applied to the vacuum character χ 0 = q −k (1 + q 2 + . . .) our differential operator must not give a term proportional to q −k+1 . This provides an extra constraint, which allows us to determine the differential operators up to I 13 . One could make further progress, for example, by calculating I L=1 2m−1 at general h using the Virasoro commutation relations, but we have not pursued this.
A complete tabulation of results is given in appendix C.
diagonal. In particular, at levels zero and one, where there is a single state, this state is necessarily an eigenstate of the KdV charges, and it makes sense to talk about the value of the KdV charge at levels zero and one.

Two-point functions
We can also use constraints from the low orders in the q expansion to fix parameters in our expression for I 2m−1 I 2n−1 . For a generic primary, we have the q expansion Thus, we can use the known values of the KdV charges on the primaries, and the values of the KdV charges on the level 1 states, which can be obtained from the differential operators for the one-point functions we worked out in the previous subsection, to constrain the form of the differential operator for the products. Note that we cannot yet determine the individual eigenvalues at level 2, as the preceding calculation only determines the sum i I i,L=2 2m−1 . This gives us enough data to determine the simplest of the quadratic expressions: and More complicated differential operators, such as I 2 5 , are not completely fixed by the constraints at levels 0 and 1 so more constraints are necessary. In this case, we can now use the fact that when there is a null state at level 2. Thus for this value of h there is just one state at level 2, and (I 2 from the previously obtained differential operator, we have one more constraint on the coefficients in I 2 5 , which allows us to determine In appendix C, we show that these expressions vanish on the minimal models as expected.

Higher point functions
We can proceed in this way for correlation functions of higher power of the KdV charges. For example, knowledge of I 2m−1 h and I 2 2m−1 h gives us enough data to determine the two eigenvalues I i,L=2 2m−1 individually. We can then use the constraints from levels 0, 1, and 2 to constrain expressions for three-point functions. In addition, the generic character has three states at level 3, except when where there is a null state at level 3 and only two states are present. Thus we can determine the I i,L=3 2m−1 for this value of h. This gives us enough data to determine In appendix C, we also calculate I 2 3 I 5 , and show that these expressions vanish on the minimal models as expected.
We can proceed in this manner to obtain expectation values of higher powers of the KdV charges, although this procedure becomes more difficult. One major obstruction is that the number of states at level n grows rapidly with n. For example, at level 4 we have 5 states, so even if we consider the value of h for which there is a null state at level 4, we still cannot determine the eigenvalues of I 3 at level 4 from knowledge of I 3 , I 2 3 and I 3 3 alone. Nevertheless it is still possible to determine I 4 3 by using constraints from minimal models which have a more complicated structure of null states. This is described, and the explicit result for I 4 3 is given, in appendix C. It would be interesting to implement this strategy more systematically, and see whether it is possible to obtain the differential operators for all KdV charges recursively in this manner. If this is not possible, then more direct (but less efficient) computations may be necessary. We describe these more direct methods in the next two sections.

Explicit evaluation by commutators
The methods described in the previous section are by far the most efficient, but can be checked using a more straightforward evaluation of the thermal expectation of KdV charges which does not rely on the detailed modular structure of the differential operator. In this and the following section we will describe two such methods for evaluating the KdV charges.
In this section we will simply consider the evaluation of the KdV charges (1.2) as polynomials in the Virasoro modes. As discussed in section 2.1, using the Virasoro commutation relations (2.1), we can simplify these polynomial expressions and rewrite them in terms of lower polynomials. The key idea is that in the terms which involve non-zero modes, we can move a generator L r past the q L 0 −k in the trace, picking up a factor of q −r from the commutator [L r , L 0 ] = rL r . Then we can move it back past the other generators in the trace to recover the original expression with a factor of q −r , plus additional terms coming from the right hand side of the commutator. We can then solve for the original expression in terms of the terms coming from the right hand side, which involve fewer generators.
Explicitly, for the simplest case of I 3 , The factors of L 0 can be obtained as derivatives of the partition function, Thus For the first term, we apply the method described above, moving the L −n around the trace: Thus, we can write where for any odd positive integer s, we define These sums are related to the Eisenstein series through where B 2k are the Bernoulli numbers. In particular, So we can finally write reproducing the result obtained in section 3.1.
We could continue to explore further examples by applying this commutation method. The calculation of I 2 3 using this method is described in appendix E. However, the calculation quickly becomes lengthy, and the appearance of the Eisenstein series appears somewhat miraculous from this point of view. In the next section we describe an alternative approach which organizes the calculation more cleanly.

Evaluation using Zhu's recursion relations
A second approach to the evaluation of thermal expectation values of KdV charges is to compute an n-point function of the stress tensor on the torus, and then perform the necessary contour integrals to extract the desired KdV charges. Expressions for the thermal n-point functions of the stress tensor can be obtained by using the recursion relations of Zhu [39], extended in [40] (the extension in [40] allows us to consider correlation functions also involving zero modes, which arise in the recursion starting from a pure operator correlation function). These recursion relations are reviewed in appendix F. The basic idea is the same as in the previous calculation: we can obtain an expression for an n-point function in terms of n − 1 point functions by taking one operator in the correlator, expanding it in modes, and commuting each of the modes around the correlator. The method of [39] uses the vertex operator algebra formalism to organize this calculation in a convenient way, which makes some of the modular properties more manifest.

One-point functions
The calculation of one-point functions in this method is fairly straightforward. For example, I 3 is the zero mode of the operator (T T )(u), where the round brackets denote conformal normal ordering, and u is a coordinate on the cylinder with compact real part. Given an expression for the thermal two-point function of T , T (u 1 )T (u 2 ) , we could implement the conformal normal-ordering by a contour integral, where C is a contour encircling u 1 , and then take the zero mode in u 1 by a further integral, where we adopt a convention that u has period 1 in the real direction, u ∼ u + 1, to match the relation between plane and cylinder coordinates usually used in the vertex operator algebra formalism.
Zhu's recursion relations provide an efficient method for computing the thermal n-point functions of the stress tensor and its derivatives. For the case of I 3 , the relevant torus amplitude in the vertex operator algebra notation is F ((ω, z 1 ), (ω, z 2 ); τ ), where we write q as q = e 2πiτ ,ω = ω − kΩ is the state corresponding to the stress tensor on the cylinder, and z i = e 2πiu i is the corresponding planar coordinate. Here ω is the state corresponding to the stress tensor operator on the plane, and Ω is the global ground state, dual to the identity operator. Note that although we are interested in calculations on the cylinder, it is conventional in the vertex operator formalism to express correlators as functions of the coordinates on the plane.
The calculation of this correlator is discussed in appendix F; the result is where P m (x, q) are the alternative Weierstrass functions introduced in [39], which are defined in appendix A. We will subsequently usually omit the q argument for compactness of notation. This correlator is a function of z 2 z 1 = e 2πi(u 2 −u 1 ) , as we would expect on general grounds; it just depends on the separation of the two operators on the cylinder. We see that it is a sum of derivatives of the partition function, whose coefficients are Weierstrass functions of ratios of the positions of the insertion points. This pattern continues for higher n-point functionsthe coefficients involve functions of the different ratios z j z i for i < j. The leading derivative term in (5.3) arises from taking the zero mode in each of the stress tensor operators, and rewriting the thermal expectation value of the zero mode in terms of derivatives, as in the previous section. Any stress tensor n-point function will similarly always start with a term with n derivatives of the partition function.
For I 3 , the integral (5.1) is just extracting the zero mode in the Laurent expansion in u 2 − u 1 . This Laurent expansion is obtained by rewriting P m (e 2πi(u 2 −u 1 ) , q) in terms of the Weierstrass function ℘ m (u 2 − u 1 , τ ) as in (A.10). Using the expressions in appendix A, we find that Thus the result is in agreement with the previous results. In this calculation, it is a little clearer how the Eisenstein series appear, as the zero modes of the Weierstrass functions under the conformal normal ordering. I 5 is the zero mode of the operator J 6 = (T (T T )) + (k + 1 6 )(2π) 2 (T T ), where again round brackets denote conformal normal ordering, and the derivatives are with respect to the cylinder coordinate u. 12 We can proceed by obtaining expressions for T (z 1 )T (z 2 )T (z 3 ) and T (z 1 )T (z 2 ) from the recursion relation. We then normal order, letting z 3 → z 2 and then z 2 → z 1 in the first term, and z 2 → z 1 in the second term. The result is the one-point function of the operator J 6 (z 1 ) , which is already independent of z 1 as expected.
This gives a result Writing this in terms of modular covariant derivatives reproduces the result in (5.7). Given an expression for the operator corresponding to I 2m−1 , it is easy to extend this calculation to find an expression for I 2m−1 . Since there is some confusion in the literature as to the precise form of these operators, we have instead used the calculation backwards; for small m, we can take a generic form for the operator with arbitrary coefficients and evaluate the one-point function, and then fix the coefficients in the expression for the operator by requiring agreement with the expressions in section 3.1. This calculation is described in appendix D.

Higher-point functions
For higher-point functions, the calculation of the correlation functions of KdV charges from the correlation functions of the stress tensors and their derivatives becomes more complicated. There are two issues; the conformal normal ordering to obtain a correlation function of the currents J 2m from the expression for the stress tensors becomes more involved, and we need to take zero modes of the expression for the current correlation function with respect to the angular coordinates on the spatial circles to obtain the correlation functions of KdV charges. The conformal normal ordering is relatively straightforward, while taking the zero modes is more troublesome. These issues are described in detail in appendix G, we will just illustrate the issues here.
Let us consider I 2 3 . Applying the recursion methods, we obtain a thermal four-point function of the stress tensor, where the coefficients are functions of the ratios z j z i with i < j. These are now products of Weierstrass functions, and the related functions g i k defined in appendix A; for example, (5.9) In each term, a given z i appears in the denominator in the argument of at most one Weierstrass function.
To calculate I 2 3 from this, we first need to conformal normal order the stress tensors in pairs, by taking z 4 → z 3 , z 2 → z 1 . For the Weierstrass functions of z 4 z 3 or z 2 z 1 , we need to make a Laurent expansion using the expression in terms of For the Weierstrass functions of other arguments, we need to make a Taylor expansion, with coefficients which are functions of z 3 z 1 . Subleading terms in the Taylor expansion may contribute if they multiply the singular terms in the Laurent expansion of functions of z 4 z 3 or z 2 z 1 . The conformal normal ordering is described in more detail in appendix G. This normal ordering will give us an expression for (T T )(z 3 )(T T )(z 1 ) , again as a sum of ∂ 4 Z . . . Z, with coefficients involving functions of z 3 z 1 . To evaluate I 2 3 , we now want to integrate over the spatial circle, that is, we want to integrate over the phase of z 3 and z 1 , which picks out the zero mode in the Laurent expansion of these functions in terms of z 3 z 1 . For I 2 3 , the functions are products of Weierstrass functions, and we can show that this zero mode can be rewritten in terms of derivatives of Eisenstein series, see appendix G. The result reproduces (3.6). Similarly, for I 3 I 5 , we can obtain a result which reproduces (3.7).
In evaluating I 2 5 , the most difficult term is the six-point function of the stress tensor. Conformal normal ordering gives us a two-point function (T (T T ))(z 1 )(T (T T ))(z 4 ) . This is a function of the ratio z 4 z 1 , which will be a sum of ∂ 6 Z . . . Z. There are three arguments in the original six-point function which become z 1 in this expression, so terms can involve products of up to three Weierstrass functions of z 4 z 1 . The terms with products of two functions we can deal with analytically, but for the terms with three functions, the zero mode in the Laurent expansion involves a double sum, which we can't analytically simplify in terms of Eisenstein series.
Similarly, in I 3 3 , conformal normal ordering the six-point function of the stress tensor, we obtain an expression for (T T )(z 1 )(T T )(z 3 )(T T )(z 5 ) . Now we need to take the zero modes in the arguments z 1 , z 3 , z 5 . In any given term, there are at most two Weierstrass functions where the denominator in the argument is z 1 , and at most two where the denominator in the argument is z 3 , so we need to deal with products of up to four Weierstrass functions. For the terms with up to three Weierstrass functions, the zero mode can be rewritten analytically in terms of derivatives of Eisenstein series. But when we have a product of four, we again encounter terms where the zero mode is a double sum. We have not found any explicit conversion of these double sums into expressions in terms of Eisenstein series.
We are able to confirm that the results of the recursion relation agree with the previous results, by expanding the double sums we encounter in a q expansion and checking order by order that they agree with a specific combination of Eisenstein series. However, these issues with the zero modes make it difficult to extend this recursion relation approach to obtain higher order results in an algorithmic way.

KdV charges for minimal models
In this section, we demonstrate that for the minimal models (2n + 1, 2) for n ≥ 1, the KdV charges I m with m divisible by 2n+1 vanish as operators. These non-unitary minimal models have central charges and have n + 1 Virasoro primary operators with conformal dimensions We have already seen evidence for this vanishing in the results for the thermal expectation values in the previous sections. The vanishing of I 2n+1 in the minimal model (2n + 3, 2) can be understood using a relation to the null states of the theory (following the work of [44,45]), which we will describe first. We will then give a general argument using integrability techniques from [4] that all the KdV charges with spin divisible by 2n + 1 indeed vanish in the (2n + 3, 2) model.

Relation to null states
This section reviews the discussions in [44,45], which elucidate the relations between the vanishing of KdV charges and the null states in the minimal models. We first review the arguments in the case of Yang-Lee model in detail, and then summarize the basic arguments for the general (2n + 3, 2) model. The detailed arguments can be found in [44,45] and references therein.
For the (5, 2) (i.e. Yang-Lee) theory, the vacuum has a level 4 null state, corresponding to the operatorJ T . 3) The KdV charge I 3 is I 3 = J 4 with J 4 = (T T ). Since the difference betweenJ 4 and J 4 is a total derivative, the KdV charge can also be written as a zero mode ofJ 4 , ButJ 4 should be set identically to zero as an operator statement since it corresponds to a null state. Hence I 3 = 0 as an operator statement in the Yang-Lee model. For the general (2n + 3, 2) theory, a similar argument can be constructed. The vacuum has null states at level 1, 2n + 2 and 2n + 5. We will call the operators corresponding to the null states at 2n + 2 and 2n + 5J 2n+2 andJ 2n+5 respectively. We would like to write I 2n+1 in terms of the zero mode ofJ 2n+2 , finding the difference between J 2n+2 (whose zero mode is I 2n+1 ) andJ 2n+2 as above. However, the explicit expression for the current J 2n+2 is complicated, so we proceed instead by calculating the commutator ofJ 2n+2 with I 3 . Roughly speaking, one finds that [ J 2n+2 , I 3 ] ∼ L −3J2n+2 . Realizing thatJ 2n+2 itself has a null state at level 3, which is preciselyJ 2n+5 ∼ L −3J2n+2 , we see that the RHS of this commutator vanishes. Thus, the zero mode ofJ 2n+2 will indeed commute with I 3 in the minimal model, so it can be identified with I 2n+1 . Hence I 2n+1 = 0 in the (2n + 3, 2) minimal model. In principle, it seems that this argument should be generalisable to show that all KdV charges with spin divisible by 2n + 1 vanishes, but we did not pursue this further.

Relation to integrability
We will now prove the general statement using the connection to integrability. There is a long-standing connection between the study of CFT in terms of the KdV charges and minimal models. In the work of Bazhanov-Lukyanov-Zamolodchikov (BLZ) [4], an approach to CFT based on quantum versions of the monodromy matrices of KdV theory was introduced. 13 Although it is impossible to review the full details here, we will try to motivate the various objects which play a role in our story, while referring the interested readers to the vast literature of integrability for more details. 14 A central object that we will need is the transfer matrix/operator T j (λ). This arises from the consideration of the Lax matrices, L and M . The Lax matrices are linear operators acting on some auxiliary internal space V with some fixed but arbitrary dimension (not necessarily related to the dimensions of phase space). In a classical system with finitely many degrees of freedom (DOF), the entries of the Lax matrices are some functions of phase space while the defining property of the Lax pair is that the equation of motion of the system is equivalent to the Lax equation dL/dt = [M, L]. From this, one can show that all moments of L, i.e. Tr V L k (with the trace being a matrix trace), are constants of motions and that they commute among themselves. 15 For a system with finitely many DOFs, we will only have finitely many conserved charges, and so some of these moments are redundant. Therefore, finding all eigenvalues of L provides all conserved charges. There are two generalizations of this formalism. First of all, one has to extend this construction to the case of classical field theory, i.e. infinitely many DOFs. To do so, in the case of finite DOFs, one constructs a generating polynomial x k TrL k which captures all conserved charges. In the case of field theory, one would need infinitely many conserved charges. So instead of a generating polynomial, it is convenient to promote the generator of conserved charges to a generating function (of some parameter called the spectral parameter λ). The Lax matrices L(λ) and M (λ), still satisfying the Lax equation, now have entries which are functions of λ, and the moments of L(λ) provide conserved charges. We then expand near λ → ∞, and the coefficients in the series expansion provide the (infinitely many) conserved charges. Finally, in a quantum field theory, we shall promote the Lax matrix to a Lax "matrix-operator". What we meant by that is that the Lax "matrix-operator" L(λ) is first and foremost a matrix in the sense of linear operator acting on internal space V , just as in the classical case. It is also a linear operator acting on the Hilbert space. We shall, following the traditions in the field, simply call this the Lax operator. The transfer matrix (or the quantum monodromy matrices) is then given by Tr V L(λ), which still remains a quantum operator. In a CFT 2 , we take V as the spin-j (where j is a positive half integer) representation of SL(2, R). The transfer matrix is then defined as the trace of the Lax operator, T j (λ) ≡ Tr j L(λ) and by definition T 0 (λ) = I [4,33].
The non-trivial result of [4] is that the T j (λ) are not all independent, but instead they satisfy "fusion relations": so T j (λ) for j > 1/2 are determined in terms of T(λ) ≡ T 1 2 (λ). This basic operator also acts as a generating function for the KdV charges, in the sense that there is an asymptotic expansion as λ → +∞ (6.8) .

(6.9)
For the (2n + 3, 2) minimal models, there are further relations: , . . . , n. (6.10) In particular, this implies that T n+ 1 2 (λ) = I, and it leads to a truncation of the fusion relations to finitely many relations. For the (2n + 3, 2) minimal models, we have 2n+3 , ξ = 2 2n + 1 . (6.11) As a notation, we shall denote t j (λ) as the eigenvalues of T j (λ) in the Fock space of the free field representation. The fusion relations are where we have defined x ≡ λ 2n+3 2n+1 . We wish to show that the KdV charges I 2k−1 vanish as operators whenever 2k − 1 is divisible by 2n + 1. We shall prove this for n = 1 and n = 2 before proceeding to the case of general n. Now it is useful to introduce f (x) ≡ log t(x) − mx, (6.15) such that the asymptotic expansion is

Example 1: Yang-Lee minimal model
In terms of f (x), the constraint is Expanding near x → ∞, the LHS is where we need to sum over k such that 2k − 1 is divisible by 3. The RHS will vanish exponentially since t ≈ e mx−O(1/x) . More explicitly, since the argument of the log is close to unity, we can use the taylor expansion of log to write Thus we conclude that I 2k−1 vanishes if 3 divides 2k − 1. Note that in this example, we take x → ∞ along real x. However, from [33], we know t could have zeroes when arg x = ±π/3. The careful way to do the calculation is to factor out the zeroes and apply the fusion relations and asymptotic expressions on the remaining object. The readers are referred to [33] for a more careful discussion.

(6.22)
We eliminate f 1 using the first equation, so the second equation becomes .

(6.23)
Expanding near infinity, the LHS nicely becomes The RHS can be expanded and is vanishing since t 1 (x) and t 1 2 (x) are exponentially vanishing. Since t 1 (x) ≈ e mϕx (6.25) even with the phases, for real x the real part of the exponent is given by Re[ϕe ±iπ/5 ] = 1 + cos 2π 5 > 0. Thus we conclude that I 2k−1 vanishes if 5 divides 2k − 1.

General n
Generalizing to arbitrary n is now straightforward. Writing the fusion relations in terms of the f j and eliminating all but f 1 where Q = e iπ/(2n+1) . Expanding at infinity, the LHS becomes while the RHS will be exponentially suppressed. Thus one concludes that I 2k−1 vanishes when 2n + 1 divides 2k − 1.
As noted in the n = 1 calculation, we take x → ∞ along real x for the asymptotic expansion while the fusion relations require a wedge around the real axis. From [33], we know that this fine up to the zeroes of t. The readers are referred to [33] for a more careful discussion about how to factor out the zeroes of t and apply the arguments to the remainder.

Application: statistics of KdV charges
We have shown that the correlation functions of the KdV charges in a particular Virasoro module can be written as modular differential operators acting on the character. In this section we discuss an application of this result. So far we have worked in canonical ensemble, where we turn on a potential q = e −β conjugate to energy. However, for a Verma module (the generic Virasoro representation with c > 1) the modular derivatives take a particularly simple form. This allows us to extract results in the micro-canonical ensemble, and in particular allows us to write exact formulas for the moments of the KdV charges as a function of level. The statistics of I 1 = L 0 − k are trivial, since I 1 has the definite value I 1 = h + n − k at level n. The statistical distribution of the higher KdV charges, however, is non-trivial. We illustrate this below by computing the mean and the variance of I 3 as a function of level.
The character of a generic Virasoro representation is We can then compute where we have introduced the shifted levelh = h − k + 1 24 to simplify subsequent formulae. For example, we have These results are for the generic character; we can obtain similar results for characters with null states by taking differences. For example the vacuum character χ vac = χ h=0 − χ h=1 , so which gives for example (7.6) If we were to study the KdV charges in a full CFT we would sum the above expressions over all representations appearing in the theory. In this case it is generally more useful to keep the expression for the KdV correlation function as a differential operator, as this does not include any dependence on the conformal dimension h of the character, so the thermal correlation function is given by the same differential operator acting on the partition function.
However, if our interest is in an individual character, the above replacement is a significant simplification. It can be used to study the statistics of the KdV charges on the P (n) states at level n in the Virasoro representation In section 3, we used knowledge of the value of the KdV charges at low level to fix the form of the differential operator in the one-point functions. We can also use this relation in the inverse sense: given the differential operator, we can act with it on the character and expand the result in a q-expansion to obtain an expression for the average value of the KdV charge at each level. This approach was already used implicitly in section 3, where we used the differential operator form of the one-point functions to determine the values of the KdV charges at level one, which was then used as input in the calculation of two-point functions. But given (7.3), we can now proceed more systematically.
For example, expanding (7.4) in a q-expansion gives lσ 1 (l)P (n − l), (7.7) where P (n) is the number of partitions of n and the divisor functions σ i (n) = d|n d i appear because they are the expansion coefficients of the Eisenstein series expansions: E 2 = 1 − 24 ∞ n=1 σ 1 (n)q n and E 4 = 1 + 240 ∞ n=1 σ 3 (n)q n . In writing this expression we have used some properties of the divisor functions.
This expression gives the sum of I 3 over all the states at level n, that is the trace of the matrix elements between states in this level. To compute the average value I 3 (n) of the KdV charge at level n we must divide by the number of states, to get: lσ 1 (l) P (n − l) P (n) .
(7.8) Applying similar analysis to the other one-point functions gives explicit expressions for the average value of the KdV charges on the states at a given level in the Virasoro representation. Similarly studying the correlation functions will give higher moments and correlations in the distribution of KdV charges on the descendant states at a given level; this is interesting information characterizing the action of the KdV charges within a given Virasoro representation.
For example, (7.10) This allows us to compute the variance ∆I 3 (n) of the KdV charge as a function of level: (7.11) Again, identities of the divisor function have lead to significant cancellations.
We will just make a two brief comments about the interpretation of these results. The first is that the mean value of I 3 (n) increases quadratically in both h and n when either of these are taken to be large. This is not a surprise, since I 3 is quadratic in the stress tensor. Then second (and less obvious) statement is that the normalized variance ∆I 3 (n)/I 3 (n) vanishes at large n. This means that as we increase the level, the distribution of KdV charges becomes very sharply peaked around its average value I 3 (n). This is illustrated in Figure 1 for typical values of h and k. We will study the statistics of KdV charges at high level in more detail in [1].

A Eisenstein series & Weierstrass functions
In our discussion we use the Eisenstein series E 2k (τ ), whose expansion is It is useful to record the formula We also recall the Ramanujan identities The Eisenstein series E 2k (τ ) for k ≥ 2 are modular forms of weight 2k. Products of E 4 and E 6 provide a basis for the space of modular forms, so all the higher Eisenstein series can be written as polynomials in E 4 , E 6 . E 2 is a quasimodular form, transforming as under the modular S transformation.
In the recursion relations, we use the Weierstrass functions defined in [39], defined for k ≥ 1, which converge for |q| < |x| < 1. For compactness, we often omit the q argument in P k . Derivatives with respect to x can be simplified using the relation If one takes the argument x = q z = e 2πiz , this becomes simply ∂ z P k (e 2πiz , q) = kP k+1 (e 2πiz , q).
In the recursion of [40], we also encounter functions defined for k ≥ 1, where we use the notation ∂ = q∂ q . For k > i, these can be rewritten as q derivatives of the P k , The P k (x, q) are related to the familiar Weierstrass functions ℘(z, τ ) by [39] where q z = e 2πiz , q = e 2πiτ , and G 2 (τ ) is an Eisenstein series, see below. Using the relation (A.7) and the similar relation ∂ z ℘ k (z, τ ) = −k℘ k+1 (z, τ ), this determines the relation for all higher k, Note that the arguments of P k are exponentials of the arguments of ℘ k . To evaluate the behaviour of P k (x, q) near x = 1, we will need the Laurent expansion of ℘ k (z, τ ) for small z, Note that the combinatorial factor in the sum vanishes for 2n + 2 − k < 0, so the only singular term in the expansion is the explicit 1 z k . The G 2k (τ ) are Eisenstein series, but with a different normalization, G 2k (τ ) = 2ζ(2k)E 2k (τ ).

B Hypergeometric form of differential operators
The differential operators appearing in the thermal expectation values may appear exotic, but -at least for the one point functions I 2m−1 -they can be put in a more familiar hypergeometric form by a change of variables.
The essential point is the following. In the above parameterization, we have represented the KdV charges as derivatives with respect to q = e 2πiτ where τ is the standard modular parameter of the torus defined by the identifications z ∼ z + 1 ∼ z + τ . However, we can alternatively parameterize a torus as an algebraic curve, as the set of solutions to the equation Here we view the torus as a double cover of the Riemann sphere C * parameterized by the w coordinate, with branch points at w = 0, λ, 1 and ∞. The two sheets of the cover are the two roots of equation (B.1). Indeed, any two-fold cover of C * which is branched over four points can put in this form, with λ being equal to the cross-ration of the four points. The parameter λ now sweeps out the torus moduli space, and is related to the usual torus modular parameter by This is the usual modular lambda function. The important point is that the modular SL(2, Z) symmetries of the torus now act on λ in a very simple way, as the anharmonic transformations generated by This means that, when written in terms of the parameter λ, all of the modular properties of our differential operators will be replaced by transformation properties with respect to (B.3).
It is straightforward to write all of the ingredients of our differential operators in terms of λ, suing When we compute the expectation value of a single power of a KdV charge all of the factors of ∂ λ θ 3 θ 3 will cancel among each other; this is a consequence of the fact that I 2m−1 is a modular form. For example, 8) It is then a straightforward exercise to put these in the standard form of a hypergeometric differential equation with regular singular points at λ = 0, 1, ∞. For example, in terms of the variableZ = (λ(1 − λ)) Z the differential operator for I 3 is the usual second order hypergeometric differential operator with parameters For quasi-modular forms like I 2 3 the results are rather more complicated, and cannot be written as hypergeometric differential operators. For example, (B.10)

C Thermal expectation values in q expansion
Here we give the full results on the determination of the differential operators using the data from the q-expansion.

C.1 One-point functions
We can write where I L=0 2m−1 is the value of the KdV charge on the primary, and I L=1 2m−1 is its value on the unique state at level one. Knowledge of I L=0 2m−1 will fix one coefficient at each order in derivatives in the differential operator determining the one-point function. The values of KdV charges on primary states were tabulated up to I 15 in appendix B of [33]. Using that data allows us to fix For higher m, this calculation leaves undetermined coefficients in I 2m−1 , as there are now combinations which vanish at leading order in q. We can obtain some further results by using the fact that for the ground state, h = 0, there is no state at level 1, χ 0 = q −k (1+q 2 +. . .), so for h = 0 we must have I L=1 2m−1 = 0. This provides one extra constraint on the parameters, and allows us to determine To go further, we would need more data, for example, from evaluating I L=1 2m−1 at general h. One can check that the values in minimal models vanish as advertised: • In the Yang-Lee (5, 2) minimal model, k = − 11 60 , and there are two primaries: the vacuum with h = 0 and a scalar of weight h = − 1 5 . The vacuum representation has a null vector at level 4. In [35], this null state was shown to imply a differential equation for the characters of this theory, This is precisely our differential operator I 3 for k = − 11 60 , so we see I 3 h = 0 in the Yang-Lee model. One can also verify that I 9 h = 0; this is most explicitly seen by using the replacement D 2 χ h = I 3 h − k 60 E 4 χ h to rewrite • The (7, 2) minimal model has k = − 17 42 . The differential operator in (3.3) for this value of k reduces to the differential operator which annihilates the characters of the (7, 2) minimal model [41].
• Similarly (C.4) reduces to the operator which annihilates the characters of the (9,2) minimal model for k = − 23 36 , (C.9) reduces to the operator which annihilates the characters of the (11, 2) minimal model for k = − 29 33 , and (C.6) reduces to the operator which annihilates the characters of the (13, 2) minimal model for k = − 175 156 .

C.2 Two-point functions
As discussed in section 3, the differential operators for two-point functions are determined by writing and using the known values of the KdV charges on the primaries, and the values of the KdV charges on the level 1 states, which can be obtained from the differential operators for the one-point functions we worked out in the previous subsection, to constrain the form of the differential operator for the products. This determines It is useful to rewrite these expressions in terms of derivatives of I 3 , using the replacement D 2 Z = I 3 − k 60 E 4 Z. This simplifies the form of the coefficients somewhat, and makes explicit the vanishing for the Yang-Lee model. We have We see that this is a combination of modular derivatives of I 3 and terms which vanish when k = − 11 60 . Thus, this is zero in the Yang-Lee model. Similarly we can write We would also expect I 3 I 5 to vanish for the (7, 2) minimal model, which can be made explicit by rewriting it in terms of derivatives of I 5 , We can also write it in terms of derivatives of I 5 : Thus, we see that it vanishes for the (7, 2) minimal model. We have also obtained the differential operator for I 3 I 7 by these methods: We have written this in a form which makes it explicit that it vanishes on the (9, 2) minimal model; it also vanishes on the Yang-Lee model.

C.3 Three-point functions
In section 3, we obtained This expression also vanishes on the Yang-Lee model, which can be made manifest by rewriting the expression in terms of modular covariant derivatives acting on I 3 , (C. 21) This has been written in a form that makes it explicit that it vanishes on the (7, 2) minimal model; it can also be shown to vanish on the Yang-Lee model.

C.4 Further constraints from particular representations
We can make some further progress by considering the action of the differential operator on the characters in particular minimal models, where the number of states at low levels is reduced. For example, we can apply the differential operators for I 3 , I 2 3 and I 3 3 to any Virasoro character, which fixes completely the eigenvalues of I 3 of this representation at any level with 3 or fewer states. We can then take the general expression for I 4 3 as a quasi-modular differential operator and consider its action on the same characters, and use this to constrain the unknown coefficients in the differential operator. We start by using the action on the general Verma module character up to level 3, as well as the vacuum character up to level 5 and the character with h at the value (3.8) (for which the representation has a null state at level 2) up to level 4. This still leaves several unfixed numerical coefficients, which can be determined by computing the action on the characters of the Ising, Tri-critical Ising, Potts and Tri-critical Potts model. The result is:

D Recovering operators from one-point functions
The KdV charges are the zero-mode of local operators constructed from the stress tensors and their derivatives. For I 2m−1 , the operator J 2m is weight 2m, and has a component with m stress tensors. The components with fewer stress tensors are determined by requiring that I 2m−1 commutes with the lower order KdV charges. 16 This is an involved calculation, and there is a relative lack of explicit expressions in the literature, so we have taken expressions for the operators with arbitrary coefficients and performed a recursion calculation to obtain their one-point functions, and then compared to the expressions we obtained in section 3.1 to fix the coefficients. This confirms the form J 6 = (T (T T )) + (2k + 1 6 )(2π) 2 (T T ), and gives us Note that these operators are only defined up to the addition of total derivative terms, as it is really the zero modes I 2m−1 that we are interested in. For J 8 and J 10 , these are useful cross-checks on the calculation, as there are fewer free parameters in the expression for the operator than in our expression for I 7 h and I 9 h . In J 12 , we are starting to get different independent operators at the same derivative order, so the expression has the same number of free parameters as I 11 h . For J 14 , there are too many independent operators for us to fix all the coefficients uniquely using the expression for I 13 h .

E Brute-Force calculation of thermal expectation values
In this appendix we describe some further calculations using the straightforward cycling of the Virasoro operators around the trace.

E.1 I 3 with a primary
We consider the thermal expectation value of a primary field V = m e 2πimu v m of dimension h, where u is the coordinate on the cylinder, with an insertion of the KdV charge I 3 , We can calculate the first term by moving L's around the trace, To calculate the last term in (E.2) we expand V in modes and use the commutation relation This gives us By translation invariance, V does not depend on the location of V on the torus, therefore the derivative terms vanish in the last expression. The result is Note that the derivative operator appearing in the first term is precisely the same one that appears in I 3 in (1.4). It is interesting that for conserved currents, with h = 1, we would get just this term. The torus one-point function of a chiral primary operator is a modular form of weight h, so it is instructive to re-write our result in terms of the Serre derivative D h = ∂ − h 12 E 2 , As with the calculations quadratic in KdV charges, we get a modular form of weight h + 4 plus E 2 times a modular form of weight h + 2.
We will describe the brute-force calculation of I 2 3 without the use of the recursion relations. Using the definition (1.2), we have We write I ≡ n≥1 L −n L n . From (4.6), we have and hence We will calculate the term n,m≥1 L −n L n L −m L m by moving the L's around the trace, For the terms involving L n−m L m−n and L −(m+n) L (m+n) we re-label the indices as m − n ≡ l and m + n ≡ l correspondingly for the first two summands (for some fixed n ≥ 1), and then split the summation ranges, The finite sums involving L −l L l , cancel and we are left with L −m L m + 8k(n 4 − n 2 )I 1 + 4k 2 (n 6 − n 2 )1 = 12n 2 I + 4n 2 I 2 1 + 8kn 4 I 1 + 4k 2 n 6 1 + − 5 6 n 2 − n 3 + 11 6 n 4 I 1 + k n 2 10 + n 4 6 − n 5 + 11 15 n 6 1 − n 2 L −n L n .
Putting everything together, and using the relation (4.9) between the sigmas and the Eisenstein series, we recover (3.6).

F Recursion relations for stress tensors and derivatives
In our second approach to the calculation of the thermal expectation values of KdV charges, we make use of the thermal n-point functions of stress tensors and derivatives of the stress tensor, calculated using the recursion relations of [39,40]. The basic idea of the recursion relation is the same as in the direct computation: commuting the modes around the trace.
We take an n-point function Tr[V 1 (z 1 ) . . . V n (z n )q L 0 −k ], (F.1) pick one of the operators in the n-point function, expand it in modes V r . For each non-zero r, we can commute the mode around the trace. The commutators with the other operators give a series of terms involving n − 1 point functions, and commuting the mode past q L 0 −k will give the original expression back multiplied by q r . Thus, we can obtain an expression for the non-zero modes V r in terms of a sum of n − 1 point functions where one of the operators in the n − 1 point function has the operator mode V r acting on it, and we sum over the n − 1 other operators. The zero mode cannot be treated in this way, so it is left as a separate contribution. Thus, the full n point function can be written as an n − 1 point function with an insertion of the zero mode V 0 , and a sum over r of the sum of n − 1 point functions where one of the operators in the n − 1 point function has the operator mode V r acting on it. To continue the recursion, we need also to commute the modes of the operators past any zero modes present in the correlation function.
In [39,40], the calculation is carried out in the vertex operator algebra formalism; to describe this we need to introduce some notation. For a state a we have a vertex operator V (a, z). We will often omit the explicit vertex operator label; so for example we write the modes of the operator V (a, z) as a n , with V (a, z) = n a n z n+ha . (F.2) We also introduce the square bracket modes where the expansion coefficients c are defined by We write the torus n-point function of the vertex operators associated to states a 1 , . . . a n as F ((a 1 , z 1 ), . . . (a n , z z ); τ ). Note that the coordinates are coordinates on the plane; for the torus, the fundamental region can be taken to be 1 < |z| < |q|. In the recursion, we will also encounter n-point functions with additional insertions of zero modes of some operators. We write F (b 0 ; (a 1 , z 1 ), . . . (a n , z z ); τ ) for the n-point function with powers of the zero mode of the vertex operator corresponding to the state b.
The key result from [40] is the recursion relation z 1 ), . . . (a n , z z ); τ ) = F (b 0 a 1 0 ; (a 2 , z 2 ), . . . (a n , z n ); τ ) For the case = 0, this reduces to the recursion relation of [39]. This has the structure sketched at the beginning of this section: we have expanded the operator V (a 1 , z 1 ) in modes a 1 m . The first term is the zero mode a 0 . The second line is a sum of n−1 point functions, with a mode a 1 [m] acting on one of the operators, summed over modes and over which operator it acts on. Note that the sum over modes a 1 m has been re-organised in terms of the square bracket modes. The functions appearing in this case are simply g 0 k = P k . When we generalise to = 0, dealing with commuting the a 1 m past the zero modes b 0 turns out to be nicely expressed in terms of the sum over i above; see [40] for details.
We want to apply this recursion relation to calculate the n-point functions involving just stress tensors and their derivatives. The stress tensor on the cylinder corresponds to the statẽ ω = ω − kΩ, where Ω is the global vacuum state, and ω = L −2 Ω is the state corresponding to the stress tensor on the plane. The zero modeω 0 = L 0 − k. The square bracket modes form a representation of the Virasoro algebra, up to normalization; the Virasoro modes are given by L [n] = (2πi) 2ω [n + 1]. Explicitly, we havẽ and so on. Note that the first expression only involves two terms, as there are just two nonzero coefficients c(2, j, 0), and the coefficient of the term involving L 2 in the RHS of the last expression vanishes.
In the recursion relation, we will have states obtained by acting with modes of the stress tensor on the stress tensor. where d n = 2k (2πi) 2n+4 n!((n + 2) 2 − 1)(n + 2). (F.13) Applying recursion will also require the action of modes of the statesω[0] nω on other states of the same form. This can be determined recursively; for m ≤ n + i + 1, Since all the states obtained are of the formω[0] nω , or the vacuum state Ω corresponding to the identity operator, the most general correlation function we need to consider in the recursion is F (ω 0 ; (ω[0] n 1ω , z 1 ) . . . (ω[0] nnω , z n ); τ ). Applying the recursion formula (F.5) to this case, we have where the upper limit on the sum on m is M = n j + n 1 + i + 1, and in the last term the jth argument of the correlator is absent because the state is the vacuum Ω, so the operator is the identity operator. This expresses the n-point function of the stress tensor in terms of n − 1-point and n − 2 point functions. This recursion can be straightforwardly implemented in Mathematica. The simplest case is the two-point function of the stress tensor. Applying the recursion relation, we get (F.18) The one-point functions of the operator on the torus get a contribution only from the zero mode of the operator, so this simplifies to (using the fact that (ω[0]ω) 0 = 0) F ((ω, z 1 ), (ω, z 2 ); τ ) = F (ω 2 0 ; τ ) + 2 (2πi) 2 P 2 z 2 z 1 F (ω 0 ; τ ) + 12k (2πi) 4 P 4 z 2 z 1 F (τ ). (F. 19) Theω 0 = L 0 − k can be rewritten in terms of q-derivatives, to give F ((ω, z 1 ), (ω, z 2 ); τ ) = ∂ 2 Z + P 2 z 2 z 1 2 (2πi) 2 ∂Z + P 4 (F.20) For the higher-point functions, if we start with an n-point function just of stress tensors, F ((ω, z 1 ), . . . , (ω, z n ); τ ), the leading derivative term comes from just taking the zero mode in each of the operators, so it will have n q-derivatives, F ((ω, z 1 ), . . . , (ω, z n ); τ ) = ∂ n Z + . . . . (F.21) If we have a more general correlator, theω[0] nω do not have zero modes, so the leading derivative order is the sum of the number of stress tensor operators plus half the number of other operators. At each stage in the recursion, we eliminate the dependence on one of the arguments z i in the correlators, generating coefficients which are functions of the ratio z j z i for j > i. Thus, the final form is a sum of products of functions g i k z j z i for j > i, where in a given term z i can only appear in the denominator of one of the functions in the product.

G Extracting the KdV correlation functions
The recursion relations give us expressions for correlation functions of the stress tensors and their derivatives. This method can be applied iteratively to efficiently obtain such expressions up to quite high order. However, to obtain the desired correlation functions of KdV charges, we have two further steps: first, to obtain a correlation function of the currents J 2m , we need to appropriately conformal normal order the stress tensors. Second, we need to integrate each of the currents over the spatial circle, taking zero modes of the expression for the current correlation function with respect to the angular coordinates. The conformal normal ordering is relatively straightforward, while taking the zero modes is more troublesome.

G.1 Conformal normal ordering
Having obtained an expression for the relevant n-point correlation functions of the stress tensor and its derivatives using the recursion relation, to extract the correlation functions of the operators corresponding to the KdV charges, we have to conformal normal order products of operators. Operators are conformal normal ordered in pairs; if the correlation function involves V (z i ) and W (z j ), we obtain a correlator involving (V W )(z i ) by taking z j → z i , and keeping the zero mode term in the Laurent expansion in u, where we write z j = e 2πiu z i . the real parts of the coordinates u i on the cylinder. This integration picks out the part independent of z i in a Laurent expansion of the resulting expression in powers of z i . Evaluating this zero mode is the most non-trivial part of our computation from the recursion relation.
For I 3 or I 5 , this final stage of the computation is trivial; after conformal normal ordering, we are considering the one-point function (T T )(z) or J 6 (z) , and this one-point function is already independent of z.
For I 2 3 , we obtain an expression for (T T )(z 3 )(T T )(z 1 ) . There were two arguments z 2 , z 1 in the original four-point function corresponding to z 1 in this correlator, so a given term can involve at most two Weierstrass functions of z 3 z 1 . Thus, we are interested in zero modes of products. The product is x n 2 , (G.6) so the zero mode comes from terms with n 2 = −n 1 , For m 1 + m 2 odd, this vanishes, as the terms with n < 0 cancel the terms with n > 0. For m 1 + m 2 even, (1 − q n 1 )(1 − q n 4 ) (n 1 + n 4 ) p 3 q n 1 +n 4 (1 + q n 1 +n 4 ) (1 − q n 1 +n 4 ) 2 n p 1 1 n p 2

4
(1 − q n 1 )(1 − q n 4 ) (n 1 − n 4 ) p 3 q n 1 (1 + q n 1 −n 4 ) (1 − q n 1 −n 4 ) 2 (1 − q n 1 )(1 − q n 4 ) (n 4 − n 1 ) p 3 q n 4 (1 + q n 4 −n 1 ) (1 − q n 4 −n 1 ) 2 . (G.18) We can then do a q expansion; comparing to the q expansions of quasimodular forms of the expected order, we were able to show that the particular combination which appears in I 3 3 can be rewritten in terms of Eisenstein series. For I 2 5 , after conformal normal ordering we have an expression for J 6 (z 1 )J 6 (z 4 ) . This can involve products of up to three Weierstrass functions of z 4 z 1 . The products of pairs of Weierstrass functions which appear are cases we have treated already. The product of three Weierstrass functions has a zero mode which involves a double sum: x n 1 +n 2 +n 3 , (G.19) so the zero mode is the terms with n 1 + n 2 + n 3 = 0. Again, we were not able to rewrite this in terms of Eisenstein series explicitly, but we were able to reproduce the q expansion of the result from a q expansion of a combination of Eisenstein series.

G.3 Tracing the modular weight through the calculation
The argument of [32], reviewed in section 2.2, implies that the results for the expectation values of the KdV charges are combinations of powers of E 2 with coefficients which are modular forms of particular weights. One of the motivations for developing the recursion relation approach was to have a method of calculation which makes this structure manifest. Unfortunately, while the relation to Eisenstein series is clearer in this approach, there are still issues.
In the recursion relation calculation for an n-point function of the stress tensor, in each term 2n = 2n ∂ + i m i , where n ∂ is the number of q-derivatives (acting either on the partition function or a Weierstrass function) and i m i is the sum of the indices of the Weierstrass functions appearing. This property is usually preserved under conformal normal ordering, as P 2k z j z i is replaced by E 2k , or contributes a singular term which, if it multiplies functions like P m z j z k , will pick out corresponding higher terms in the Taylor expansion. However, there is an exception; the −iπ term in (G.1) does not obey this rule.
In the zero modes of operators, the zero mode mostly has a modular weight corresponding to the sum 2n ∂ + i m i for the product whose zero mode we're taking, but not always. In the simple products of pairs of Weierstrass functions this rule of thumb is obeyed: the resulting zero mode in (G.9) has weight m 1 + m 2 , and similarly the zero mode in (G.10) has weight m 1 + m 2 + 2. However, in (G.11) and (G.12), we have combinations of terms of weight m 1 + m 2 + 4 and m 1 + m 2 + 2. Similarly, in the product of triples of Weierstrass functions, (G.13) has weight m 1 + m 2 + m 3 as expected, but (G.14) has weight m 1 + m 2 + m 3 − 1, and (G.15) is a combination of weight m 1 + m 2 + m 3 + 2 and m 1 + m 2 + m 3 . The zero mode in the product of four Weierstrass functions also has contributions of weight m 1 + m 2 + m 3 + m 4 and m 1 + m 2 + m 3 + m 4 − 2.
Thus, the expression for a KdV charge obtained from the n-point function of the stress tensor can have components of different modular weight. The first place these difficulties arise is in the computation of I 3 3 , where there are components of modular weight 12 and 10, but the components of modular weight 10 cancel. This cancellation is required to reproduce the expected structure, but it is unfortunate that the structure is not more manifest.