Local charges in involution and hierarchies in integrable sigma-models

Integrable $\sigma$-models, such as the principal chiral model, ${\mathbb{Z}}_T$-coset models for $T \in {\mathbb{Z}}_{\geq 2}$ and their various integrable deformations, are examples of non-ultralocal integrable field theories described by (cyclotomic) $r/s$-systems with twist function. In this general setting, and when the Lie algebra ${\mathfrak{g}}$ underlying the $r/s$-system is of classical type, we construct an infinite algebra of local conserved charges in involution, extending the approach of Evans, Hassan, MacKay and Mountain developed for the principal chiral model and symmetric space $\sigma$-model. In the present context, the local charges are attached to certain `regular' zeros of the twist function and have increasing degrees related to the exponents of the untwisted affine Kac-Moody algebra $\widehat{{\mathfrak{g}}}$ associated with ${\mathfrak{g}}$. The Hamiltonian flows of these charges are shown to generate an infinite hierarchy of compatible integrable equations.


Introduction
One of the hallmarks of integrability in a (1 + 1)-dimensional field theory is the existence of an infinite number of conserved charges. At the classical level, this property can be attributed to the existence of a Lax connection, depending on some auxiliary complex spectral parameter λ, whose zero curvature equation for all λ is equivalent to the equations of motion of the field theory. In principle, the conserved charges can all be obtained by expanding the monodromy of the Lax connection in λ around specific points. Depending on the chosen point of expansion, the resulting charges can be either local or non-local in the fields entering the Lax connection.
When passing to the Hamiltonian formalism, one requires additionally that the conserved charges be in involution with respect to the Poisson bracket of the theory, namely that they Poisson commute not only with the Hamiltonian but also between themselves. It is known since the mid-eighties that a sufficient condition guaranteeing the involution of the charges built from the monodromy is that the Poisson bracket of the Lax matrix L(λ, x), the spatial component of the Lax connection, be of Maillet's general r/s-form [1,2]. In most cases of interest, the pair of matrices r 12 (λ, µ) and s 12 (λ, µ), giving the form its name, are rational functions on C 2 valued in the two-fold tensor product g⊗g of some finite-dimensional Lie algebra g. The mathematical formalism underlying this particular case was pinned down by Semenov-Tian-Shansky in [3] where the r-and s-matrices were understood to be the skew-symmetric and symmetric parts of a single solution R 12 (λ, µ) of the classical Yang-Baxter equation. In general, the latter is related to the standard solution R 0 12 (λ, µ) on the (twisted) loop algebra over g by where ϕ(λ) is a rational function [4][5][6] (see also [7]), known as the twist function. The loop algebra over g is twisted by some automorphism σ of order T ∈ Z ≥1 , with the non-twisted case corresponding to T = 1. Recall that the twist function plays an essential role in characterising a given integrable field theory. In fact, it was shown recently in [8] to form an integral part of the Lax matrix itself, viewed as a rational function valued in the untwisted affine Kac-Moody algebra associated with g.
The principal chiral model serves as the prototypical example of an integrable σ-model which fits the general formalism of r/s-systems with twist function [5]. As such, it has proved extremely fruitful over the past couple of years to try and reinterpret some of its properties in terms of its twist function. Indeed, once a given property has been understood algebraically at the level of the twist function, it almost immediately generalises to other integrable field theories described within this formalism. The first illustration of this general philosophy came about from the desire to generalise the Faddeev-Reshetikhin construction [9], initially developed in the context of the SU(2) principal chiral model, to other integrable σ-models with twist function. Specifically, it was shown in [10] that the key initial step of this construction can be naturally reformulated in the language of twist functions. This led to a proposal for extending the Faddeev-Reshetikhin approach to a wide range of other models of interest, including the symmetric and semi-symmetric space σ-models, as well as the Green-Schwarz superstring on AdS 5 × S 5 in [11]. The prospect of extending the subsequent steps in the Faddeev-Reshetikhin construction to these other models remains an exciting open problem. It is interesting to note that, despite its usefulness in the case of the principal chiral model, the relevance of the twist function was first appreciated in [6], following [12,13], on the much more elaborate AdS 5 × S 5 superstring both within the Green-Schwarz [14] and the pure spinor [15] formulations.
Another important application of the formalism of the twist function ties in with the great effort made in recent years [16][17][18][19][20][21][22][23][24][25][26][27] towards deforming some well known integrable field theories, such as the principal chiral model as well as symmetric and semi-symmetric space σ-models, while preserving their integrability. Specifically, it was realised in [18] that the so called Yang-Baxter σ-model, first introduced by Klimčík in [16] as a certain one-parameter deformation of the principal chiral model on any real Lie group G 0 , could be naturally obtained by deform-ing the poles of the twist function of the principal chiral model. This led to an immediate broadening of the landscape of Yang-Baxter type deformations to include also one-parameter deformations of the symmetric and semi-symmetric space σ-models, incorporating, in particular, the Green-Schwarz superstring on AdS 5 × S 5 in [21]. In fact, many other deformations were also understood a posteriori to arise in this fashion [28]. It is worth noting in passing that the bi-Yang-Baxter σ-model [17,29] is special in this regard. Although it was originally devised as a two-parameter deformation of the principal chiral model on any real Lie group G 0 , it can equally be regarded as a two-parameter deformation of the symmetric space σ-model on G 0 ×G 0 /G 0,diag with G 0,diag the diagonal subgroup of G 0 ×G 0 . It was shown in [30] that it is this latter formulation which fits within the general framework of r/s-systems with twist function. The original description of the model as a double deformation of the principal chiral model can be recovered by fixing the G 0,diag gauge symmetry, at the expense of losing the formulation in terms of a twist function [30] (see, however, subsection 7.2 below).
Aside from providing a systematic way of constructing integrable deformations, the idea of deforming the pole structure of the twist function of a given integrable field theory has also played a pivotal role in establishing and characterising the symmetry algebras of the resulting deformed models. Specifically, the principal chiral model and (semi-)symmetric space σ-models all have a double pole in their twist function which splits up into a pair of simple poles when their Yang-Baxter type deformation is switched on. It was shown in [18,21], based on earlier work [31][32][33][34][35][36][37][38] (see also more recent related results [39][40][41]), that the charges extracted from the leading order in the expansion of the monodromy at this pair of simple poles satisfy all the relations of a Poisson algebra U q (g), the semiclassical counterpart of the quantum group U q (g) with q = q . This general feature of Yang-Baxter type deformations of double poles in the twist function was subsequently related to Poisson-Lie G-symmetries in [42]. Amongst the charges spanning the Poisson algebra U q (g), those associated to non-Cartan generators are all non-local. In the example of the Yang-Baxter σ-model, the level zero charges together with two additional non-local charges coming from the next order in the expansion of the monodromy around the simple poles of the twist function, have been shown [43] to satisfy all the defining relations of the semiclassical counterpart of the quantum affine algebra U q ( g).
The purpose of the present article is to provide another application of the general formalism of r/s-systems with twist function. Specifically, we will describe how, in this general framework, infinite towers of local charges can be associated with certain zeros of the twist function, all of which are in pairwise involution. Following the same spirit as recalled above, the starting point of our approach was to reinterpret the construction of local charges in the principal chiral model due to Evans, Hassan, MacKay and Mountain [44] in the present language of twist functions. In fact, this construction had soon been generalised to include also the (supersymmetric) principal chiral model with a Wess-Zumino term in [45], symmetric space σ-models in [46] as well as supersymmetric coset σ-models in [47]. Each of these generalisations can be regarded as further evidence that such a construction should hold for any integrable field theory with twist function, while at the same time providing indications on how to do so. In the remainder of this introduction we will briefly summarise the main results of the paper.
Let us first note that in all of the integrable σ-models with twist function described above, every zero of ϕ(λ) is such that ϕ(λ)L(λ, x) is regular there. In a general integrable field theory with twist function ϕ(λ) we shall say that any zero of ϕ(λ) with this property is regular. We denote by Z the set of regular zeros of ϕ(λ) in C. As discussed in subsection 7.4, the regularity property of the zeros of the twist function in an integrable σ-model is related to a general condition used for describing these models as dihedral affine Gaudin models [8]. We shall further distinguish between two types of zeros: cyclotomic ones and non-cyclotomic ones. This notion depends on the order T of the automorphism σ. In a model with T = 1, every point is by definition non-cyclotomic, whereas in a model with T > 1, every point is non-cyclotomic except for the origin and infinity. As explained in subsection 2.2, throughout our analysis the point at infinity will be treated in much the same way as the origin by using an inversion of the spectral parameter.
To every λ 0 ∈ Z, or every λ 0 ∈ Z ∪ {∞} if infinity is also a regular zero, we will associate a subset of integers E λ 0 ⊂ Z ≥2 and a corresponding tower of local charges Q λ 0 n labelled by n ∈ E λ 0 . The first main property of these charges which we will establish is that any two such charges Q λ 0 n and Q µ 0 m for any λ 0 , µ 0 ∈ Z and n ∈ E λ 0 , m ∈ E µ 0 are in involution. Moreover, if infinity is a regular zero and either λ 0 or µ 0 is taken to be the point at infinity, the corresponding local charges will only Poisson commute up to a certain field C(x) which will coincide with the coset constraint in Z T -coset σ-models. Following the standard terminology from the theory of constrained Hamiltonian systems, we will refer to equalities as being weak when they hold only after setting this particular field to zero, see subsection 5.1. Furthermore, we show that in every example of integrable σ-model considered, the Hamiltonian can be expressed as a particular linear combination of the collection of quadratic local charges Q λ 0 2 for λ 0 ∈ Z ∪ {∞} and the momentum of the model. It then follows that all of the local charges are conserved.
Let us briefly outline the construction of the local charges by considering first the case when λ 0 ∈ Z is non-cyclotomic. If the Lie algebra g is of type B, C or D then the density of the local charge Q λ 0 n is obtained simply by evaluating Tr ϕ(λ) n L(λ, x) n (1.1) at the regular zero λ 0 . When g is of type A, on the other hand, the density of the local charge Q λ 0 n is given instead by a certain polynomial in the above expressions, determined as in [44] with the help of a generating function. In either case, E λ 0 is given here by the set of exponents of the affine Kac-Moody algebra g associated with g, shifted by one (we do not treat the case of the Pfaffian in type D). In the example of the principal chiral model on a real Lie group G 0 treated in [44], the twist function has simple zeros at ±1 and the evaluation of ϕ(λ)L(λ, x) at λ = ±1 produces the currents j ± = g −1 ∂ ± g of the theory, where g is the G 0 -valued field of the principal chiral model and ∂ ± are the partial derivatives along light-cone coordinates on the worldsheet. We recover in this way the higher spin local charges in involution of the principal chiral model constructed in [44].
When the regular zero λ 0 ∈ Z is cyclotomic, i.e. λ 0 = 0, it may happen, as a result of the equivariance properties of both the Lax matrix and twist function, that the evaluation of (1.1) at the point λ 0 vanishes identically. More precisely, the first non-vanishing term in the power series expansion of (1.1) around λ = 0 is of order λ rn for some 0 ≤ r n ≤ T − 1. If the Lie algebra g is of type B, C or D, or also of type A with an inner automorphism σ, then we define the density of the local charge Q 0 n as the coefficient of this leading term. The case when g is of type A and the automorphism σ is not inner is treated in a similar fashion to the case of a non-cyclotomic point in type A, with the densities of the local charges Q 0 n being obtained by means of a generating function. In each case it turns out that we need to restrict attention to indices n such that 0 ≤ r n < T − 1. As a result, and in contrast to the case of a non-cyclotomic regular zero, some exponents of the affine Kac-Moody algebra g are 'dropped' in the construction of the subset E 0 , specifically those such that r n = T − 1. In the case of a symmetric space σ-model, for which T = 2 so that only charges for which r n = 0 are kept, we recover in this way the local charges found in [46].
The collection of local charges Q λ 0 n , λ 0 ∈ Z, n ∈ E λ 0 in involution generate an infinite set of Poisson commuting Hamiltonian flows Q λ 0 n , · on the phase space of the model. To every such flow we then associate a corresponding g-valued connection ∇ λ 0 n = Q λ 0 n , · + M λ 0 n (λ, x) for some g-valued matrix M λ 0 n (λ, x) depending on the spectral parameter λ. The second main property of the local charges Q λ 0 n , λ 0 ∈ Z, n ∈ E λ 0 which we establish is that the connection ∇ λ 0 n for any λ 0 ∈ Z and n ∈ E λ 0 commutes with the connection ∇ x = ∂ x + L(λ, x). In this sense, the local charges generate a hierarchy of integrable equations. We use this result to deduce that the local charges Q λ 0 n , λ 0 ∈ Z, n ∈ E λ 0 are in involution with the non-local charges extracted from the monodromy of L(λ, x). Moreover, we go on to show that when g is of type B, C or D, any two such connections ∇ λ 0 n and ∇ µ 0 m for λ 0 , µ 0 ∈ Z and n ∈ E λ 0 , m ∈ E µ 0 also commute with one another. Finally, we have also checked these results in the case of type A for low values of n and m and on this basis we conjecture it to hold in general. If infinity is a regular zero then the majority of these results still hold in the weak sense when we consider also the local charges associated with infinity.
This article is organised as follows. The general framework of r/s-systems with twist function which we employ throughout the article is reviewed in section 2. In particular, we introduce the notion of a regular zero in the complex plane which plays a central role in our analysis. In subsection 2.1, we define the R-matrices entering the r/s-systems of interest and discuss their equivariance properties, as well as those of the Lax matrix and the twist function. We will make extensive use of these properties when discussing local charges extracted from cyclotomic regular zeros of the twist function. The list of examples we shall consider is given in paragraph 2.1.4. In subsection 2.2, we define the notion of a regular zero at infinity and relate it to that of a regular zero at the origin by inversion of the spectral parameter. Finally, we establish some general results in subsection 2.3. Section 3 is devoted to the procedure for extracting local charges in involution in the case of a non-cyclotomic regular zero. In particular, we present in subsection 3.4 an explicit construction of the currents K λ 0 n for type A algebras using gen-erating functions in the spirit of [44]. Section 4 deals with charges at cyclotomic zeros. We explain how the equivariance properties of the various objects affect the construction of local conserved charges in involution. Here the Lie algebras of type B, C and D can still be treated uniformly but in type A we need to consider separately the cases when the automorphism σ is inner or not. The generating function for Lie algebras of type A with non-inner automorphism is presented in subsection 4.6. A list of properties of these local charges is collated in section 5, including the fact that the local charges extracted from different regular zeros Poisson commute (weakly when the point at infinity is involved). Moreover, we show that all the local charges commute with the field C(x) which will play the role of the constraint in Z T -coset σ-models, therefore showing that they are gauge invariant. We also discuss the reality conditions of all the local charges. The Hamiltonian flows of the local charges Q λ 0 n are studied in detail in section 6. The main result that any two of the g-valued connections ∇ λ 0 n and ∇ µ 0 m satisfy a zero curvature equation is established in subsection 6.3. Finally, in section 7 we apply all these results to the examples listed in paragraph 2.1.4. In paragraph 7.4 we interpret some of our results and assumptions in terms of dihedral affine Gaudin models [8]. We end with an outlook and two appendices.
2 Framework and general results 2.1 Non-ultralocal models with twist function 2.1.1 Non-ultralocal models In this section we begin by outlining the general framework which we will work in. We consider a two dimensional field theory with spatial coordinate x taking values on the circle S 1 or the real line R. The dynamics of the model is described by a Hamiltonian H and a Poisson bracket {·, ·} on the phase space. Let P denote the conserved momentum of the model, whose Poisson bracket generates the derivative ∂ x with respect to x on the phase space.
We suppose that the model is integrable with Lax matrix L valued in a Lie algebra g and depending on a complex spectral parameter λ. Moreover, we also suppose that the Poisson bracket of the Lax matrix with itself assumes the form of Maillet's r/s-system, i.e.
using the standard tensorial notations i. In this equation, R 12 is a g⊗g-valued matrix depending on the spectral parameters λ and µ, δ xy is the Dirac distribution and δ ′ xy = ∂ x δ xy . When the R-matrix is skew-symmetric, the term containing δ ′ xy vanishes and the model is said to be ultralocal. In general, R is non skew-symmetric and the model is then non-ultralocal.
In this article, we will focus on the case where the Lie algebra g is simple. More precisely, we will restrict to the classical types A, B, C and D of the Cartan classification, seen in their defining representations 1 : Let us recall the family of solutions of this equation that we shall consider. Let T be a positive integer, σ : g → g an automorphism of g of finite order T and ω a primitive T th root of unity. We fix a basis T a of g and its dual basis T a normalised such that Tr(T a T b ) = δ a b . We define the Casimir tensor on g to be where here and throughout we use summation convention on repeated Lie algebra indices. The standard R-matrix is then a solution of the CYBE (2.2). Note that σ 1 C 12 = σ −1 2 C 12 , as the Killing form on g is σ-invariant. In this article, we will consider a matrix R obtained by "twisting" this cyclotomic matrix, namely where ϕ is a rational function, called the twist function of the model. It is easy to check that R is also a solution of the CYBE (2.2). We shall call a Poisson bracket (2.1) with such an R-matrix an r/s-system with twist function. As σ T = Id, the eigenvalues of σ are of the form ω p where, by convention, we take p ∈ {0, . . . , T − 1}. We denote by g (p) the corresponding eigenspace and by π (p) the projection on g (p) in the direct sum g = T −1 p=0 g (p) . One then has the identity Here J n is the standard symplectic structure on C 2n given by J n = 0 Id −Id 0 and M T denotes the transpose of M .

Equivariance properties
As σ is of order T , it defines an action of the cyclic group Z T = Z/T Z on g. On the other hand, Z T can be seen as acting on the complex numbers C via multiplication by ω. We then remark that the matrix R 0 is equivariant under these two actions, in the sense that We will suppose that the Lax matrix L possesses a similar equivariance property, namely The compatibility of these two properties with the Poisson bracket (2.1) imposes that from which we deduce that λϕ(λ) is invariant under the action of Z T . Thus, there exists a rational function ζ such that λϕ(λ) = ζ(λ T ). (2.10) In this article, we will be interested in the zeros of the twist function ϕ in C. We will say that such a zero λ 0 is regular if ϕ(λ)L(λ, x) is holomorphic at λ = λ 0 . By virtue of the equivariance properties (2.8) and (2.9), if λ 0 is a regular zero, all points of the orbit Z T λ 0 are also regular zeros. Let us pick (arbitrarily) one of them. We then form a set Z of regular zeros of ϕ such that for every pair of distinct points λ 0 and µ 0 in Z, the orbits Z T λ 0 and Z T µ 0 are disjoint. As explained in subsection 2.2, we will also be interested in the case where the differential form ϕ(λ)dλ has a zero at infinity, i.e. where has a zero at α = 0. We will also see that the appropriate notion of a regular zero at infinity corresponds to requiring that 1 α ϕ 1 α L 1 α , x be holomorphic at α = 0.

Examples
To end this subsection we list some examples of models which fit the framework of this article. We consider a real Lie group G 0 whose Lie algebra g 0 is a real form of g. The space of fields valued in the cotangent bundle T * G 0 is naturally equipped with a canonical Poisson bracket. One can parametrise this phase space by two fields g and X, respectively G 0 -valued and g 0valued. The Poisson bracket then reads Models with T = 1. The first class of examples of non-ultralocal models with twist function consists of the principal chiral model (PCM) and its integrable deformations (dPCM). Among these models are the so-called Yang-Baxter deformation (or η-deformation) [16][17][18], the PCM with Wess-Zumino term and the combination of these two deformations [23]. The Lax matrix of these models all have the form The g 0 -valued fields j 0 and j 1 are expressed in terms of (g, X), in a model dependent way. In the simplest case, the PCM, we have For the deformed models, these relations are modified and involve two parameters η and k (which are zero for the PCM). In all these models, the Lax matrix (2.13) satisfies the Maillet bracket (2.1) with an Rmatrix of the form (2.5). The matrix R 0 is given by equation (2.4) for T = 1 (and thus σ = Id) and the twist function is [23] where A is a certain function of η and k, which vanishes when η = 0. The Hamiltonian and the momentum of the model are given by with B a global factor depending on A and k via the relation B = − 1 4 ϕ ′ dPCM (1)ϕ ′ dPCM (−1). The twist function ϕ dPCM has two zeros at +1 and −1. Moreover, these are regular zeros, i.e. ϕ(λ)L(λ, x) is regular at λ = ±1. Note that the evaluation of the latter at ±1 gives There exists another two-parameter deformation of the PCM, the so-called bi-Yang-Baxter model [17]. It depends on two parameters η andη, such that the case η =η = 0 corresponds to the PCM. It is integrable [29,30] and has a Lax matrix which satisfies a Poisson bracket of Maillet type (2.1). The associated R-matrix is not of the form (2.5) and so the bi-Yang-Baxter model does not quite fit within the above framework. However, we will see in subsection 7.2 how the method discussed here can be adapted to apply also to the bi-Yang-Baxter model.
Models with T > 1. The second class of examples that we will consider are the σ-models on Z T -coset spaces [48]. Consider first the case of symmetric spaces, i.e. of Z 2 -coset spaces, which possess a family of one-parameter deformations (that we shall call dZ 2 models). These models are constructed as follows. Let σ be an involutive automorphism of G 0 and let H = G σ 0 denote the fixed-point subgroup of G 0 under the action of σ. The symmetric space is then the quotient G 0 /H. Differentiating σ at the identity, we induce an automorphism of g 0 . We extend it linearly to the complex algebra g and thus obtain an automorphism of g of order 2, that we shall still denote σ. As a vector space, g then decomposes as the direct sum of the eigenspaces g (0) and g (1) , corresponding to the eigenvalues +1 and −1 (cf. paragraph 2.1.2).
The symmetric space σ-model and its deformation are expressed in terms of fields A (k) and Π (k) valued in the eigenspaces g (k) , for each k = 0, 1. These fields are defined in terms of the canonical fields g and X in a model dependent way (in particular, this definition depends on the deformation parameter η). The Lax matrix of the model (deformed or not) is then [18] The Poisson bracket of this Lax matrix with itself is of Maillet type (2.1) with a R-matrix of the form (2.5). The matrix R 0 is given by equation (2.4) with T = 2 and σ the above involutive automorphism of g. The twist function is [18] (2.18) The regular zeros of these models are thus the origin 0 and the infinity ∞.
To conclude the list of examples let us discuss briefly the general Z T -coset model, for any T ∈ Z ≥2 . Let σ be an automorphism of the complexified group G, of order T . The coset space is then G 0 /H, where H = G σ ∩ G 0 . We will still denote by σ the corresponding automorphism of g, which induces an eigenspace decomposition g = T −1 p=0 g (p) (see paragraph 2.1.2 for details). The phase space of the model is expressed in terms of fields A (k) and Π (k) , for k = 0, . . . , T − 1, which are defined in terms of the canonical fields g and X and belong to g (k) . The Lax matrix is then [49,8] Note that in this equation, and in general, we consider the exponents (k) only modulo T , so that A (T ) = A (0) for example. All Z T -coset models possess a gauge symmetry under the action of the subgroup H of G 0 . In the Hamiltonian formulation presented above, A (0) plays the role of a gauge field and Π (0) is the constraint associated with the gauge symmetry. In those models also, the Poisson bracket of the Lax matrix (2.19) with itself is non-ultralocal with a twist function. The associated matrix R 0 is the one in equation (2.4), with the automorphism σ introduced above. The twist function is [49,8] (2.20) As for the Z 2 -coset, the regular zeros are 0 and ∞.

Infinity and inversion of the spectral parameter
In this article we will construct a tower of local charges associated with each regular zero of the twist function. As mentioned in paragraph 2.1.3, the set of regular zeros can include the point at infinity, although the sense in which infinity can be a regular zero is slightly different from the definition of finite regular zeros. In this subsection, we show how the notion of a regular zero at infinity is related to that of a regular zero at the origin through inversion of the spectral parameter, i.e. by the change of parameter λ → α = λ −1 . Under such a change of spectral parameter we have where ψ(α) is defined in equation (2.11). Suppose that infinity is a zero of the twist function, i.e. that ψ(0) = 0, and define We will say that infinity is a regular zero if P (α, x) is regular at α = 0. In the remainder of this subsection we will assume this to be the case. We then set From the equivariance properties (2.8) and (2.9) of L and ϕ, we deduce that C is valued in the grading g (0) . Let us note here that in the Z T -coset models, described in paragraph 2.1.4, this field C coincides with the gauge constraint Π (0) . Starting from the Poisson bracket (2.1) and using the form (2.5) of the R-matrix, we find Using the expression (2.4) of R 0 , we have As P (α, x) and αψ(α) are regular at 0, taking the limit α → 0 in the above Poisson bracket, we then obtain Applying the same kind of reasoning we also find Let us define a new Lax matrix From the fact that [C (k) 12 , Z 2 ] for any Z ∈ g (0) , we find that Using this identity and the Poisson brackets (2.1), (2.24) and (2.25), we prove that the Poisson bracket of L with itself is also of the r/s-form, namely We now define The following theorem is the main result of this subsection.
Proof. Using equation (2.6), we find that 12 . (2.29) The theorem follows from the Poisson bracket (2.26) and the identity which is a consequence of equation (2.29).
To interpret Theorem 2.1, let us note that the matrix R 0 21 is nothing but the matrix R 0

12
for the automorphism σ −1 . Moreover, from the equivariance properties (2.8) and (2.9), we find that the corresponding properties of L ∞ and ψ are The Poisson bracket of L ∞ is thus an r/s-system with twist function ψ, automorphism σ −1 and spectral parameter α = λ −1 . Moreover, the point α = 0 is a regular zero of this r/s-system. Indeed, we supposed that α was a zero of ψ(α) and one can check explictly that ψ(α)L ∞ (α, x) is regular at α = 0.
It is worth noting that the procedure just described is involutive, in the following sense. If ϕ(λ)L(λ, x) is regular at λ = 0, one can check that α = ∞ (which corresponds to λ = 0) is a regular zero of the r/s-system of L ∞ and, moreover, that the corresponding field C ∞ obtained by evaluating λ −1 ψ(λ −1 )L ∞ (λ −1 , x) at λ = 0 is equal to C. Re-inverting the spectral parameter α to λ = α −1 , we can thus construct a "new" Lax matrix L ∞ (λ −1 , x) − λψ(λ −1 ) −1 C(x). According to Theorem 2.1, this Lax matrix should satisfy an r/s-system with twist function ϕ and automorphism σ. A direct computation reveals that this Lax matrix is actually equal to the initial Lax matrix L.
Let us end this subsection by illustrating the inversion of spectral parameter on the example of Z T -coset models. As noted above, for these models the field C coincides with the constraint Π (0) . After performing the change of spectral parameter λ → α = λ −1 , we find a twist function Note that the property ψ(α) = −ϕ(α) is also true for the twist function (2.18) of the η-deformed Z 2 -model. The new Lax matrix is Comparing this to the initial Lax matrix (2.19), we see that it simply corresponds (up to a minus sign on terms involving Π (k) ) to changing every grading (k) to (T −k), which is equivalent to considering the automorphism σ −1 instead of σ.

Poisson brackets of traces of powers of L
Recall that we consider the Lie algebra g in its defining matrix representation (see Table 1). We may therefore take powers of elements of g and traces of these matrices. In the following sections, we will extract local charges in involution from the traces of powers of the Lax matrix L. In this subsection, we will establish general results on the Poisson brackets of powers of L and their traces.
Suppose that X and Y are g-valued quantities such that Then the Poisson brackets of powers of X and Y are where, for t = a, b, c, we defined Proof. The Poisson bracket being a derivation, we can use the Leibniz rule yielding We conclude observing that X k 1 Y l 2 and X n−1−k 1 Y m−1−l 2 commute with X 1 and Y 2 and using the identity Corollary 2.3. Suppose that X and Y are g-valued quantities such that Then we have Proof. Starting with Lemma 2.2, the corollary follows from the cyclicity of the trace and the vanishing of traces of commutators.
Let us now apply these results to the Lax matrix L. We work in the framework described in subsection 2.1. We define S n (λ, x) = ϕ(λ) n L(λ, x) n (2.33) and T n (λ, x) = Tr S n (λ, x) . (2.34) Starting with the Poisson bracket (2.1) and the expression (2.5) of the R-matrix, we apply Corollary 2.3. We find that

Charges at non-cyclotomic zeros
The purpose of this section is to describe the procedure for extracting local charges in involution from non-cyclotomic regular zeros of the twist function ϕ. Let us first explain what we mean here by a non-cyclotomic point. If T = 1, i.e. if σ = Id and there is no cyclotomic invariance, we define any point as being non-cyclotomic. If T ∈ Z >1 , a non-cyclotomic point is a point which is not fixed by the action of the cyclic group Z T , i.e. which is not the origin or infinity.
Throughout this section we fix a non-cyclotomic regular zero λ 0 . We will focus here on the case where λ 0 is different from infinity. The case λ 0 = ∞ is treated by the same method, just replacing L by L ∞ and ϕ by ψ (cf. subsection 2.2). The fact that λ 0 is a regular zero implies that S n (λ, x) and T n (λ, x), defined in equations (2.33) and (2.34), are both holomorphic at λ = λ 0 . Thus, we can define the current Let us briefly comment on the explicit expression of these currents in the case of the PCM. As explained in paragraph 2.1.4, the PCM has two regular zeros at +1 and −1. The corresponding currents are J ±1 n,PCM (x) = Tr j n ± (x) , where j ± (x) = 1 2 j 1 (x) ± j 0 (x) . These currents are the one investigated in [44], from which local charges in involution for the PCM are constructed. In this section, we will follow the method developed in [44], generalising it to any current (3.1) associated with a non-cyclotomic regular zero λ 0 of the model.

Poisson algebra of the currents
We begin by computing the Poisson bracket of the currents J λ 0 n (x) and J λ 0 m (y). Specifically, we would like to evaluate equation (2.35) at λ = µ = λ 0 . Since λ 0 is a regular zero, S n−1 (λ 0 , x) and S m−1 (λ 0 , y) are well defined. Thus, it remains to determine U 12 (λ 0 , λ 0 ). Starting with the definition (2.36) of U and using ϕ(λ 0 ) = 0, one has U 12 (λ, λ 0 ) = ϕ(λ)R 0 12 (λ, λ 0 ). Recall from equation (2.4) that R 0 12 (λ, λ 0 ) is not regular at λ = λ 0 , so that we cannot simply evaluate the above equation at λ = λ 0 . However, as λ 0 is a non-cyclotomic point, the matrix R 0 has the following local behaviour where A λ 0 12 (λ) is regular at λ = λ 0 . Using again ϕ(λ 0 ) = 0, we then obtain where ϕ ′ denotes the derivative of ϕ with respect to the spectral parameter λ. Thus, one has Recall the completeness relation for Z in g. We cannot directly apply this identity to equation (3.3) as S m−1 (λ 0 , y) does not belong to g in general (recall that S m−1 is defined as the (m − 1) st power of a matrix in g).
Following [44], we will show in the next subsections how to circumvent this difficulty. We will treat separately the case where g is of type B, C or D and the case where g is of type A.

Type B, C and D algebras
Let us first consider the case where g is of type B, C or D, i.e. where g is an orthogonal or a symplectic algebra (cf. Table 1). One can check that, for these algebras, if X belongs to g, X n also belongs to g if n is odd. Moreover, all matrices in g are traceless. We then deduce that the currents J λ 0 n are zero for n odd. Thus, we will only extract local charges from the traces of even powers of L, i.e. from the currents J λ 0 2n . The Poisson bracket of such currents is given by equation (3.3). The right hand side contains Tr 2 C 12 S 2m−1 (λ 0 , y) 2 , and since 2m − 1 is odd we have S 2m−1 (λ 0 , y) ∈ g. Hence, we can apply the completeness relation (3.4), which yields Using the definition (2.33) of S, one has (3.6) Define the local charges where the integration is over the whole domain of the spatial coordinate x (i.e. the real line R or the circle S 1 ). Once integrated over y, the right hand side of (3.6) is a total derivative with respect to x. Assuming the periodicity of the fields if x ∈ S 1 or that they decrease at infinity if x ∈ R, we then conclude that Q λ 0 2n , Q λ 0 2m = 0. In conclusion, we have constructed a tower of local charges Q λ 0 2n in involution, as integrals of the currents J λ 0 2n (x). These currents are polynomials in the fields appearing in the Lax matrix L(λ, x). More precisely, the current J λ 0 2n is a homogeneous polynomial of degree 2n.
Up to a global factor, the Poisson bracket (3.6) is the same as the bracket (4.16) of [44]. Thus, we can apply the methods developed in [44]. In particular, this allows to construct a more general tower of local charges Q λ 0 2n (ξ) in involution, depending on a free parameter ξ ∈ R. These charges are defined as integrals of some currents K λ 0 2n (ξ). These currents are given by homogeneous polynomials in the J λ 0 2k 's, depending on the free parameter ξ ∈ R. In particular, the first currents K λ 0 2n (ξ) are given by: The expression of the current K λ 0 2n (ξ) is determined (up to a global factor) recursively from equation (3.6) by demanding that the charge Q λ 0 2n (ξ) be in involution with all the charges Q λ 0 2m (ξ) (m = 2, . . . , n − 1) constructed thus far. It can also be found without recursion with the help of a generating function, which allows a general proof of the involution of the charges Q λ 0 2n (ξ): we refer the reader to the subsection 3.4 for more details.
Taking ξ = 0 in equation (3.8), we get K λ 0 2n (ξ = 0) = J λ 0 2n . Hence, we recover the local charges Q λ 0 2n introduced in equation (3.7) as a special case of this one-parameter family of local charges. For different parameters ξ and ξ ′ , the towers of charges Q λ 0 2n (ξ) and Q λ 0 2n (ξ ′ ) do not commute with one another. We thus have to work with a fixed value of ξ: in the rest of this article, we will mainly focus on the simplest case ξ = 0. This choice is justified first by simplicity, but also because the proof of the existence of an integrable hierarchy associated to the charges Q λ 0 2n (ξ), presented in section 6, works only for the case ξ = 0.

Type A algebras
Let us now consider the case where g is of type A, i.e. where g = sl(d, C) for some d ∈ Z ≥2 (see Table 1). If X ∈ g, we have Tr(X) = 0 by definition, but in general X n / ∈ g and Tr(X n ) = 0 for n ≥ 2. Thus, we consider the currents J λ 0 n for n ≥ 2. The Poisson bracket between two such currents is given by equation (3.3). Since in general S m−1 (λ 0 , y) does not belong to g, we cannot use the completeness relation (3.4) to simplify this equation. However, a variant of the identity (3.4) exists for any matrix Z ∈ M d (C). Indeed, using the facts that Z − 1 d Tr(Z)Id belongs to g and that Tr 2 (C 12 ) = 0, we find that Applying this relation to equation (3.3) and using the identities Integrating both sides over x and y, we see that the right hand side does not vanish identically as it did in subsection 3.2. Nevertheless, following the method of [44] we will be able to construct new currents K λ 0 n such that the charges Poisson commute with one another.
The Poisson bracket (3.10) is to be compared to equation (4.5) of [44], from which it differs only by an overall factor. We can therefore directly apply the procedure developed in [44] to the present case so as to construct the desired currents K λ 0 n 's. The expression for the first K λ 0 n 's read These currents are similar to the currents K λ 0 n (ξ) described in (3.8) for g of type B, C or D. More precisely, the current (3.12) coincide with the currents K λ 0 n 1 d , recalling that for type B, C and D, the J λ 0 2k+1 's vanish. As for K λ 0 n (ξ) in type B, C and D, the expression of the current K λ 0 n for type A is determined (up to a global factor) recursively from equation (3.10) by demanding that the charge Q λ 0 n be in involution with all the charges Q λ 0 m (m = 2, . . . , n − 1) constructed thus far. However, in the present case, one does not have the freedom of a free parameter ξ in the definition of K λ 0 n : there is a unique tower of charges in involution Q λ 0 n . As in the case of type B, C and D algebras, the current K λ 0 n (x) is a homogeneous polynomial of degree n in the fields appearing in the Lax matrix L(λ, x). And as explained in [44], the degrees n for which the current K λ 0 n (x) is non-zero are the exponents of the untwisted affine Kac-Moody algebra g plus one.
At this stage, we do not have a proof that the recursive algorithm described above can be applied indefinitely. We shall now recall from [44] how to construct explicitly the current K λ 0 n without a recursive algorithm, using generating functions.

Generating functions
In the previous subsections 3.2 and 3.3, we introduced currents K λ 0 n (ξ) (for types B, C and D) and K λ 0 n (for type A), constructed recursively from the currents J λ 0 n (and which depended on a free parameter ξ for types B, C and D). In this subsection, we will show how to construct these currents using generating functions.
We will mainly focus on the case where g is of type A and will briefly comment on types B, C and D at the end of the subsection. Let us then suppose that g = sl(d, C), so that we can use the notations and results of subsection 3.3. We introduce and so that A(λ, µ, x) = exp F (λ, µ, x) . By expanding the matricial logarithm in (3.13) as a power series in µ one finds We are interested in the evaluations of F (λ, µ, x) and A(λ, µ, x) at λ = λ 0 , which are well defined as λ 0 is a regular zero. Following [44], we look for for some rational number p n , where f (µ)| µ n denotes the coefficient of µ n in the power series expansion of f (µ). The Poisson brackets of the currents T n (λ 0 , x) = J λ 0 n (x) are given by equation (3.10). This allows one to compute {F (λ 0 , µ, x), F (λ 0 , ν, y)} and {A(λ 0 , µ, x), A(λ 0 , ν, y)}. As equation (3.10) coincides with the equation (4.5) of [44] up to a global factor, these Poisson brackets are the same as in [44] (equations (4.13) and (4.14)), still up to the global factor. Thus, the procedure of [44] applies and we conclude that the Poisson bracket of the local charges (3.11) defined in terms of the currents (3.16) is It follows that the charges Q λ 0 n are in involution if we choose, for any k ∈ Z ≥2 , p k = k−1 d . The corresponding currents are given by One can check that the first currents defined by this generating function are given by equation (3.12), up to overall global factors. The current K λ 0 n (x) is the evaluation at λ = λ 0 of the more general current which we will need later. The equation allows one to compute W n (λ, x) as a polynomial in the T k (λ, x). More precisely, W n is related to the T k 's in the same way that K λ 0 n is related to the J λ 0 k 's.
We end this subsection by saying a few words on Lie algebras g of type B, C or D. In this case, we saw in subsection 3.2 that the local charges in involution can be taken as integrals of currents K λ 0 2n (ξ), depending on a free parameter ξ (see equation (3.8)). These currents can be obtained from the J λ 0 2k 's using a generating function, similar to the one presented above for type A. We will not enter into details here and will just present the final result, based on reference [44]. The current K λ 0 2n (ξ) can be computed as: Starting from the Poisson bracket (3.6), one can show that the corresponding charges Q λ 0 2n (ξ) are in involution, using similar techniques as above for type A. We refer the interested reader to reference [44] for details on the proof. An explicit computation shows that the first currents K λ 0 2n (ξ) obtained from the above equation are given by equation (3.8), up to overall global factors.

Summary
To conclude this section, let us summarise the results that we obtained. In particular, we will use this as an opportunity to extend the notations K λ 0 n and W n , defined for a type A algebra in the previous subsections, to other types. This will serve to uniformise the notation in the rest of the article.
When g is of type A, the currents K λ 0 n (x) are given in subsection 3.4 through equation (3.18). We also defined a current W n (λ, x) depending on the spectral parameter λ in equation (3.20). For a Lie algebra g of type B, C or D (as treated in subsection 3.2), we introduced currents K λ 0 n (ξ), depending on a free parameter ξ. However, as explained at the end of susbection 3.2, we will only use the currents J λ 0 n (x) = K λ 0 n (ξ = 0, x) in the rest of this paper. In order to employ uniform notations throughout the paper, we shall define in this case With these conventions, independently of the type of g, the current K λ 0 n (x) is the evaluation of W n (λ, x) at λ = λ 0 and the charge Q λ 0 n is given by Recall also that we restrict the degrees n of the currents K λ 0 n to some subset E λ 0 of Z ≥2 . In fact, independently of the type of g, E λ 0 can (almost) be seen as the set of exponents of the affine algebra g plus one. This was already observed for type A in subsection 3.3, based on the results of [44]. For types B, C and D, we saw in subsection 3.2 that E λ 0 is the set of all even numbers, which turns out to coincide with the exponents of g plus one for types B and C [44]. For type D, there are some exponents missing in this construction (the rank modulo the Coxeter number), which are related to the Pfaffian (see [44]). Although we do not consider the Pfaffian here, we expect that it should be possible to construct a corresponding local charge in the present framework too.
Having introduced these type-independent notations, we can summarise the results of this section by the following theorem.
Theorem 3.1. Let λ 0 be a non-cyclotomic regular zero of the model. Then, for any m and n in E λ 0 , the charges Q λ 0 n and Q λ 0 m are in involution, i.e. we have {Q λ 0 n , Q λ 0 m } = 0. The notations and results summarised above will be generalised to the case of cyclotomic zeros in the following section. Let us note here that there will be some subtlety in the definition of the current W n for an algebra of type A in the case when the automorphism σ is inner, compared to the definition given above. We shall discuss this in subsection 4.7.

Charges at cyclotomic zeros
In this section, we explain how to construct towers of local charges in involution attached to cyclotomic regular zeros of the twist function ϕ. Recall that a cyclotomic point is a point fixed by the action of the cyclic group Z T , i.e. the origin or infinity. Suppose we are considering a model with a regular zero at infinity. As explained in subsection 2.2, working in the new spectral parameter α = λ −1 and with the new Lax matrix L ∞ amounts to treating, instead, a model with a regular zero at α = 0 and automorphism σ −1 . Hence it is sufficient to describe the extraction of local charges at the origin.
Throughout this section we therefore consider a model with T > 1 and a regular zero at the origin. We thus have ϕ(0) = 0 and ϕ(λ)L(λ, x) regular at 0. Using the equivariance property (2.9), we see that the smallest power of λ in ϕ is of the form αT − 1, for some α ∈ Z ≥1 . In terms of the function ζ, defined in equation (2.10), this implies that ζ(λ T ) = O(λ αT ). We will mostly need the fact that ζ(λ T ) = O(λ T ), i.e. that ϕ(λ) = O(λ T −1 ), and more precisely the asymptotic property Recall that in the previous section we extracted local charges by evaluating the traces of powers of ϕ(λ)L(λ, x) at the regular zeros. In the case of a cyclotomic point, this method is not sufficient to extract all charges, as such traces can vanish. To understand how to construct the whole algebra of local charges, we will first need to establish equivariance properties of S n (λ, x).

Equivariance properties
Recall the equivariance properties (2.9) and (2.8) of ϕ and L. In this subsection, we look for a similar relation for S n (λ, x). In general, S n (λ, x) does not belong to the Lie algebra g since it is defined as the power of an element of g seen in the fundamental representation. Thus, one cannot consider directly the action of σ on S n (λ, x).
We refer here to the discussion of appendix A. We will restrict to the case where σ is not one of the special automorphisms of D 4 = so (8, C). In this case, we can extend σ to a linear endomorphism on the space F of all matrices acting on the fundamental representation, that we shall still denote σ (see details in appendix A). Note that this new endomorphism σ of F is still of order T . We will also need the following properties of σ. For any Z ∈ F we have for some ǫ in {1, −1}. Note that ǫ is always 1 except when g = sl(d, C) and σ has a non-trivial outer part, in which case T is even. We shall write ǫ = ω ηT 2 , with η in {0, 1}. From equations (2.9) and (2.8) and the identity (4.2a), we deduce that S n satisfies the equivariance property with κ = 1 + ηT 2 . Let us consider the power series expansion We then find σ(A n,r ) = ω r+κ(n−1)+1 A n,r . Thus, Tr(A n,r ) vanishes except if r ≡ r n [T ], where r n is the remainder of the euclidian division of −nκ by T . We define J 0 n (x) = Tr A n,rn (x) . In particular, the first term in the power series expansion of where ζ was defined in equation (2.10). Finally, let us remark that equation (4.6) implies that T n (λ, x) has the following equivariance property T n (ωλ, x) = ω rn T n (λ, x). (4.7)

Poisson algebra of the currents
One can extract the Poisson brackets of the currents J 0 n (x) and J 0 m (y) as the coefficient of λ rn+rm in the power series expansion of {T n (λ, x), T m (λ, y)}. The latter can be computed from equation (2.35). Specifically, using the identity (2.6) we find with ζ defined in equation (2.10). Taking the limit µ → λ we obtain The first term of this Poisson bracket has the same structure as the Poisson bracket (3.3). The main difference coming from cyclotomy is thus the second term, which involves the partial Casimir C (0) 12 . We recall that we have the partial completeness relation for any Z ∈ g.
The second term in (4.9) will therefore involve the projection S To determine these projections, we can make use of the power series expansion (4.4) and equation (4.5). In particular, one finds that A n−1,r+T is in F (0) if and only if A n−1,r also belongs to F (0) . Let us then define q n to be the unique integer between 0 and T − 1 such that A n−1,r belongs to F (0) if and only if r ≡ q n [T ]. Using equation (4.5) we find q n ≡ r n + 1 [T ]. So q n = r n + 1 if r n ≤ T − 2 and q n = 0 if r n = T − 1.
To simplify the Poisson bracket (4.9), we will need to distinguish between three cases: • g is of type B, C or D, • g is of type A and σ is inner, • g is of type A and σ is not inner.

Algebra of type B, C or D
We first consider g to be of type B, C or D. Recall that in this case S 2n−1 (λ, x) belongs to the Lie algebra so that T 2n−1 (λ, x) is zero and hence we consider only the currents J 0 2n (x). Moreover, we can use the completeness relations (3.4) and (4.10) in (4.9). We then find After integration over y, the first term becomes a total derivative with respect to x by virtue of equation (3.5) and thus vanishes when integrated over x.
Recall moreover that the Poisson bracket of J 0 2n (x) with J 0 2m (y) is obtained from (4.11) by keeping only the term λ r 2n +r 2m in the power series expansion. We note that the smallest power of λ in the second term of (4.11) is αT − 2 + q 2n + q 2m (cf. equation (4.1) and above). As we saw in the previous subsection, In the case where r 2n and r 2m are different from T − 1, the smallest power of λ is then r 2n + r 2m + αT so the second term of (4.11) does not contribute to the Poisson bracket of J 0 2n (x) with J 0 2m (y), as α ≥ 1.
If r 2n or r 2m is equal to T − 1 then there will be a contribution from this term involving other objects than only the J 0 k 's, preventing us from constructing charges in involution. Thus, we will only consider the currents J 0 2k (x) such that r 2k = T − 1. We then have We have thus extracted a tower of local charges in involution from the Lax matrix around the origin. Just as in the non-cyclotomic case, these charges are integrals of some polynomials of even degrees in the fields appearing in the Lax matrix. The main difference with the noncyclotomic case is the fact that, in general, we do not have a current of any even degree. More precisely, we 'dropped' the currents of degree 2n, for all n such that r 2n = T − 1. Recall from appendix A that in the case of an algebra of type B, C or D, we have ǫ = 1 and κ = 1. Thus r 2n is the remainder of the euclidian division of −2n by T , which means that r 2n = T − 1 if and only if 2n ≡ 1 [T ]. In particular, we see that there is no drop of any degrees if T is even.

Algebra of type A and σ inner
Let us now suppose that g is sl(d, C) and σ is inner. In this case, we have the generalised completeness relation (3.9). Moreover, we also have a similar identity for the partial Casimir C (0) 12 , derived as follows. Recall that for any Z ∈ F , Z − 1 d Tr(Z)Id belongs to g. Moreover, we note that the identity Id is in the grading zero F (0) for σ inner (cf. appendix A). Using equation (4.10), we then have Using equations (3.9) and (4.12) in the Poisson bracket (4.9), we obtain The Poisson bracket of J 0 n (x) with J 0 m (y) is obtained by extracting the coefficient of λ rn+rm in the above equation. To treat the second term on the right hand side of this equation, we follow the discussion of the previous subsection 4. 3. The smallest power of λ appearing in this term is αT − 2 + q n + q m and if we restrict to n and m such that r n and r m are different from T − 1, this power is strictly greater than r n + r m . The term then does not contribute to the Poisson bracket {J 0 n (x), J 0 m (y)}. Let us turn to the third term on the right hand side of equation (4.13). It can be seen from equation The smallest power of λ that can appear in this term is thus 2T − 2 + r n−1 + r m−1 , which is always greater than 2T − 2 and therefore strictly greater than r n + r m if r n and r m are different from T − 1.
In conclusion, only the first term of the right hand side of (4.13) contributes to the Poisson bracket {J 0 n (x), J 0 m (y)}, which then has the same structure as in the previous subsection. Integrating this bracket over x and y, we recognise the integral of a total derivative proportional to ∂ x T n+m−2 (λ, x), which then vanishes, assuming appropriate boundary conditions. Thus, for any n and m such that r n and r m are different from T − 1, we have . As in the subsection 4.3, we have ǫ = 1 and κ = 1 for σ inner. It follows that the integers n such that r n = T − 1 (for which we do not consider the charge Q 0 n ) are the ones equal to 1 modulo T .

Algebra of type A and σ not inner
Finally, let us treat the case where g = sl(d, C) and σ not inner. In particular, this implies that T is even and we shall write T = 2S in this subsection. We still have the generalised completeness relation (3.9). As σ is not inner, we have σ(Id) = −Id and hence π (0) (Id) = 0. We deduce that in this case, the partial completeness relation (4.10) actually holds for any Z ∈ F . Equation (4.9) then gives We follow the method of the previous subsections and look for the power r n +r m of λ in the right hand side of this bracket. As explained in subsection 4.3, the second term does not contribute when we restrict to r n and r m different from T − 1.
The first term is treated as in the case of a non-cyclotomic point: The powers of λ appearing in the power series expansion of the first term are then of the form r n+m−2 − 2 + aT , with a ∈ Z ≥1 . One has r n+m−2 ≡ r n + r m + 2 [T ] and 0 ≤ r n+m−2 ≤ T − 1. Moreover, r n + r m + 2 is always between 0 and 2T − 2 if we suppose r n and r n different from T − 1. If 0 ≤ r n + r m + 2 < T , we have r n+m−2 = r n + r m + 2 and the powers r n+m+2 − 2 + aT are then all strictly greater than r n + r m . If T ≤ r n + r m + 2 ≤ 2T − 2 then r n+m−2 = r n + r m + 2 − T and the power r n+m−2 − 2 + aT is equal to r n + r m if and only if a = 1.
Finally, let us consider the third term on the right hand side of (4.14). The powers of λ in its power series expansion are of the form r n−1 + r m−1 − 2 + aT , with a ∈ Z ≥1 . Note that r k−1 ≡ r k +1+S ≡ r k +1−S [T ]. Thus, r k−1 = r k +1+S if 0 ≤ r k < S −1 and r k−1 = r k +1−S if S − 1 ≤ r k ≤ T − 1. We then conclude that the power r n−1 + r m−1 − 2 + aT is equal to r n + r m if and only if r n + 1 − S ≥ 0, r m + 1 − S ≥ 0 and a = 1.
Combining all the above results, we find a closed expression for the Poisson bracket of the currents J 0 n (x) and J 0 m (y) when r n and r m are different form T − 1, specifically As in the case of a non-cyclotomic zero, one can construct charges in involution by taking integrals of new currents K 0 n (x), constructed as polynomials of the currents J 0 m (x) for m ≤ n. The method in the present case is similar: one can construct the currents K 0 n recursively by asking that the corresponding charge Q 0 n Poisson commutes with Q 0 m for all m < n. One of the main difference with the non-cyclotomic case is the fact that we do not consider a current K 0 n when r n = T − 1 (we say that such a degree n drops from the construction). The second difference is the presence of the terms θ k in the Poisson bracket (4.16) compared to (3.10). Since these terms depend on the numbers r k , the construction of the K 0 n 's will depend on T . As an illustration, we give here the expression of the first K 0 n 's for T equal to 2 and 4. In the case T = 2, we have κ = 2 hence r n = 0 for all n ≥ 2. Thus, there is no drop of any current due to the condition r n = T − 1 and the current J 0 n (x) is simply the evaluation of T n (λ, x) at λ = 0. Moreover, since 2 − T and 1 − S are both zero for T = 2, we note that all θ k terms in the Poisson bracket (4.16) are equal to 1. The construction for the K 0 n 's is therefore the same as in the non-cyclotomic case and their expression is given by (3.12).
Let us now consider the case T = 4. We find κ = 3 and r n ≡ n [4] and therefore drop the currents J 0 4k+3 (x). Constructing the K 0 n 's recursively we find where we dropped the currents K 0 3 and K 0 7 and more generally all currents K 0 4k+3 .
Comparing these currents to the ones constructed in the non-cyclotomic case (3.12), one can observe a pattern in the cyclotomic procedure. Here also, the current K 0 n is constructed by correcting J 0 n with monomials J 0 m 1 . . . J 0 mp such that m 1 + . . . + m p = n. We observe that not all the monomials in the non-cyclotomic corrections appear among the cyclotomic ones but the ones that do have the same coefficients (for example 15/4d for the J 2 J 4 correction of J 6 ). Moreover, we note that a monomial J 0 m 1 . . . J 0 mp appearing in the non-cyclotomic procedure is also present in the cyclotomic expression if and only if r m 1 + . . . + r mp = r n .
The above observations are still found to hold for larger values of T and n (although we do not include the corresponding expression of K 0 n for conciseness). This allows one to find the currents K 0 n (x) without going through the recursive procedure if one already knows the result for the non-cyclotomic case. A more systematic approach to constructing higher conserved charges in involution in the cyclotomic case would be to find a generating function for the K 0 n 's, generalising the results of subsection 3.4. This will be the subject of the next subsection.

Generating function for type A with σ not inner
In subsection 3.4 we presented, following [44], the generating function for constructing the currents K n (x) in the non-cyclotomic setting. In particular, we found that the relation between the K n (x)'s and the J m (x)'s is given by equation (3.18). In the previous subsection we showed how the currents K 0 n (x) could be constructed in the cyclotomic case from the knowledge of the corresponding result in the non-cyclotomic case. In particular, starting from the expression of K n as a polynomial of the J m 's, we observed that K 0 n can be constructed in the same way by keeping monomials J 0 m 1 . . . J 0 mp with the same coefficient if and only if r n = r m 1 + . . . + r mp .
This procedure for going from the non-cyclotomic to the cyclotomic setting has a natural interpretation in terms of equation (3.18). Indeed, the current K 0 n constructed above is equal to where the projection onto the term λ rn ensures that we keep only the monomials satisfying the condition r n = r m 1 + . . . + r mp .
Recall that λ r k J 0 k (x) is the first term in the power series expansion of T k (λ, x). Moreover, the next terms start with a (r k + T ) th power of λ. Since r n < T , such terms can be added to the exponent in equation (4.18) without changing the left hand side as they cannot contribute to a λ rn -term. We may therefore also write In terms of the definitions (3.13) and (3.14) of F (λ, µ, x) and A(λ, µ, x) and equation (3.15), we can further re-express equation (4.19) as 20) or again Starting from equation (4.14), using equation (4.15) and the identity f (y)δ ′ where We want to compute the Poisson brackets between the charges Q 0 n defined as integrals of the currents (4.20). To begin with, note that the Poisson bracket between F (λ, µ, x) and F (λ, ν, y) can be obtained from equations (3.15) and (4.21). We then find that Up to a global factor and treating the spectral parameter λ as an external parameter, Ω nm (λ, x, y) has the same structure as the right hand side of equation (3.10), which as we saw already had the same structure as equation (4.5) of [44]. This equation (4.5) is used in [44] to compute the Poisson brackets of the generating functions (equations (4.13) to (4.15)). The method developed in [44] for computing these Poisson brackets then applies to the terms Ω kl in equation (4.23) and gives a similar result, up to the global factor and the dependence on λ that we keep. Specifically, we find where h nm (µ, ν) was defined in equation (3.17).
The first term on the right hand side of (4.24), proportional to h nm (µ, ν), vanishes when p k = k−1 d for all k ∈ Z ≥2 , as expected. It therefore remains to show that the second term also does not contribute when we restrict to the (r n + r m )-th power of λ. From equation (4.22b), we see that the powers of λ appearing in the power series expansion of ∆ kl (λ, x, y) are of the form q k + q l − 2 + aT , with a ≥ α ≥ 1 and q n defined in subsection 4.2.
Recall from subsection 4.2 that q k ≡ r k + 1 [T ], and therefore q k + r n−k ≡ r n + 1 [T ]. Suppose now that r n and r m are different from, and so in particular strictly less than, T − 1. As q k + r n−k is always positive and congruent to r n + 1 modulo T , which is strictly less than T , it then follows that q k + r n−k ≥ r n + 1. Similarly, we have q l + r m−l ≥ r m + 1 and we thus deduce that q k + r n−k + q l + r m−l − 2 + cT is greater than or equal to r n + r m + T . We deduce that the term (4.26) cannot contribute to the (r n + r m ) th power of λ, as required.
In conclusion, we have found that, for any n and m such that r n and r m are different from T − 1, one has {Q 0 n , Q 0 m } = 0, where the charge Q 0 k is defined as the integral of the current K 0 k (x) given by (4.20). Recall that the current J 0 n (x) is constructed as the coefficient of λ rn in the power series expansion of T n (λ, x). Similarly, one can rewrite equation (4.20) as

Summary
Let us summarise the results of this section, as we did for non-cyclotomic zeros in subsection 3.5. In general, we define a charge Q 0 n , associated with a current K 0 n (x), by where the definition of W n (λ, x) and r n depends on the type of g and whether σ is inner or not. For g of type B, C or D (see subsection 4.3) and for g of type A and σ inner (see subsection 4.4), we simply choose W n (λ, x) = T n (λ, x), so that K 0 n (x) = J 0 n (x). In this case, we consider r n as the remainder of the euclidian division of −n by T . When g is of type A and σ is not inner (the case discussed in subsections 4.5 and 4.6), we choose W n (λ, x) as given by equation (3.20). In this case, r n is defined as the remainder of the euclidian division of −n 1 + T 2 by T .
The equivariance properties (4.7) and (4.25) imply that, independently of the type of g and of σ being inner or not, we have for W n defined as above. It therefore follows that the powers of λ appearing in W n (λ, x) are of the form r n + kT , with k ∈ Z ≥0 . In particular, the current K 0 n (x) of equation (4.28) is the coefficient of the smallest power of λ in W n (λ, x).
As in the non-cyclotomic case, we restrict the degree n of the currents K 0 n (x) to some subset E 0 of Z ≥2 . More precisely, n belongs to E 0 if n − 1 is an exponent of the affine algebra g and r n is different from T − 1 (with the exception of the exponents related to the Pfaffian in type D, as in subsection 3.5). The results of this section can be summarised as the following theorem. There is a subtlety in the definition of W n (λ, x) for g of type A and σ inner. Indeed, in this case the current K 0 n (x) is extracted just from T n (λ, x) as recalled above. Yet in section 3, the current K λ 0 n at a non-cyclotomic point was extracted instead from W n (λ, x) which differs from T n (λ, x). Thus, for this case, we choose the appropriate definition of W n (λ, x) depending on whether the regular zeros of the considered model are cyclotomic or not.
We end this section by an open question. For a non-cyclotomic regular zero λ 0 and g of type B, C or D, we considered local charges in involution Q λ 0 n built as the integral of the currents J λ 0 n (x) (see subsection 3.2). However, based on the results of [44], we also exhibited a more general family of local charges in involution Q λ 0 n (ξ), depending on a free parameter ξ ∈ R and whose corresponding currents K λ 0 n (ξ, x) are constructed as polynomials in the J λ 0 k 's. In the subsection 4.3 of the present section, where we deal with a cyclotomic regular zero at the origin for g of type B, C or D, the charges Q 0 n are also constructed simply as integrals of the currents J 0 n (x). It is thus natural to ask if there exists in this case a more general family of charges Q 0 n (ξ), depending on a free parameter ξ, as for the non-cyclotomic case. Moreover, for g of type A and σ inner (as treated in subsection 4.4), the charges Q 0 n are also integrals of the currents J 0 n (x) (we do not need to construct more complicated currents to obtain charges in involution). It is thus also natural to look for a more general family of charges Q 0 n (ξ). This would be an interesting result, as it would exhibit an important qualitative difference between the non-cyclotomic case and the cyclotomic one (with σ inner), for g of type A.
We expect these one-parameter families (for g of type B, C and D or for g of type A with σ inner) to exist. More precisely, we expect them to be given by the first non-zero coefficient in the power series expansions of a suitable generating function, depending on ξ, around the cyclotomic regular zero λ = 0. As for the non-cyclotomic case, we will focus in this article on the local charges which do not depend on a free parameter ξ (as described at the beginning of this subsection), for the same reasons as the ones discussed at the end of subsection 3.2 for a non-cyclotomic zero.

Algebra of local charges in involution
In the previous sections, we constructed a tower of local charges in involution at every regular zero λ 0 ∈ Z. More precisely, we constructed currents K λ 0 n (x), with degrees n in some subset E λ 0 of Z ≥2 , such that the charges Q λ 0 n defined as the integral of K λ 0 n (x) are in involution with one another. In this subsection, we study the whole algebra of local charges in involution, formed by all the Q λ 0 n for λ 0 ∈ Z and n ∈ E λ 0 . More precisely, we prove that currents K λ 0 n (x) and K µ 0 m (y) extracted at different regular zeros are in involution. We establish this result for regular zeros in Z, excluding the point at infinity. If infinity is a regular zero then one can also extract charges in involution Q ∞ n , using the results of subsection 2.2. The Poisson brackets of the corresponding currents with the currents at finite regular zeros involve some subtleties and will be treated separately at the end of the subsection.
In general, the currents K λ 0 n (x) and K µ 0 m (y) are constructed as polynomials of the currents J λ 0 n (x) and J µ 0 m (y). It is therefore sufficient to prove that the Poisson bracket of J λ 0 n (x) and J µ 0 m (y) is zero. The currents J λ 0 n (x) and J µ 0 m (y) are extracted from T n (λ, x) and T m (µ, y), whose Poisson bracket is given by equation (2.35). We can suppose that µ 0 is different from 0 and thus is a non-cyclotomic point, so that J µ 0 m (y) = T m (µ 0 , y). Using the fact that U 12 (λ, µ 0 ) = ϕ(λ)R 0 12 (λ, µ 0 ) since ϕ(µ 0 ) = 0, we can evaluate equation (2.35) at µ = µ 0 to find We will now treat separately the cases λ 0 cyclotomic or λ 0 non-cyclotomic.
Suppose that λ 0 is non-cyclotomic so that J λ 0 n (x) is simply the evaluation of T n (λ, x) at λ = λ 0 . Recall from paragraph 2.1.3 that, by construction of the set Z, as λ 0 and µ 0 are different elements of Z, the cyclotomic orbits Z T λ 0 and Z T µ 0 are disjoint. In particular, this means that R 0 12 (λ, µ 0 ) is regular at λ = λ 0 . Indeed, by equation (2.4) the poles of R 0 12 (λ, µ 0 ) are the points of the orbit Z T µ 0 . Moreover, S n−1 (λ, x) is regular at λ = λ 0 and we have ϕ(λ 0 ) = 0. Thus, evaluating equation (5.1) at λ = λ 0 we find that the currents J λ 0 n (x) and J µ 0 m (y) are in involution, as expected.
Let us now treat the cyclotomic case where λ 0 = 0. In this case, J 0 n (x) is the coefficient of λ rn in the power series expansion of T n (λ, x) (cf. subsection 4.1). The Poisson bracket of J 0 n (x) with J µ 0 m (y) is then the coefficient of λ rn in equation (5.1). Recall from section 4 that for n ∈ E 0 , we have r n < T − 1. Yet, in equation (5.1), R 0 12 (λ, µ 0 ) and S n−1 (λ, x) are regular at λ = 0 and ϕ(λ) = O(λ T −1 ), hence the involution of J 0 n (x) and J µ 0 m (y). In conclusion, we have proved Theorem 5.1. Let λ 0 , µ 0 ∈ Z and let n ∈ E λ 0 and m ∈ E µ 0 . Then, if λ 0 = µ 0 , we have Combining this theorem with the results of previous sections, we conclude that the local charges Q λ 0 n , for all λ 0 ∈ Z and n ∈ E λ 0 , are in involution with one another.
We now turn to the case where one of the regular zeros is the point at infinity. In this case, the current J ∞ m (y) is extracted from the Lax matrix L ∞ (α, y). From the Poisson brackets (2.1) and (2.24), we get with the matrix R defined in (2.27).
In the following, we will say that an equation is true weakly, and we will then use the symbol ≈ instead of =, if the equation holds when one puts the field C(x) to zero. This denomination and its interest will be made clear when studying Z T -coset models, in which case the field C(x) is interpreted as a gauge constraint. Note, in particular, that the last term of equation (5.2) vanishes weakly. Using Corollary 2.3, we may then compute the Poisson brackets of T n (λ, x) with T ∞ m (α, y) = Tr S ∞ n (α, x) = Tr ψ(α) n L ∞ (α, x) n weakly, to find . Suppose first that λ 0 ∈ Z is non-cyclotomic. We have ϕ(λ 0 ) = 0, and hence The Poisson bracket between J λ 0 m (x) and J ∞ m (y) is then obtained, weakly, by extracting the coefficient of α rm in the equation above. For m ∈ E ∞ , we have r m < T −1. Yet ψ(α) = O(α T −1 ) and R 0 21 (α −1 , λ 0 ) and S ∞ m−1 (α, y) are regular at α = 0. Thus J λ 0 n (x) and J ∞ m (y) are weakly in involution.
It remains to consider the case where λ 0 = 0. In this case, J 0 n (x) is the coefficient of λ rn in T n (λ, x) and we restrict to n such that r n < T − 1. We have to find the coefficient of λ rn α rm in equation (5.3). Due to the presence of ϕ(λ) or ψ(α) in the two terms appearing in V 12 (λ, α), we see that either the power of λ or the power of α is greater than T − 1 and thus cannot contribute to the term λ rn α rm . In conclusion, we have the following theorem.
Theorem 5.2. Suppose that infinity is a regular zero of the model. Let λ 0 ∈ Z, n ∈ E λ 0 and m ∈ E ∞ . Then we have Combining this theorem with the results of previous sections, we see that the local charges Q λ 0 n , for all λ 0 ∈ Z ∪ {∞} and n ∈ E λ 0 , are (at least weakly) in involution with one another.
We thus constructed a large algebra of local charges in involution, composed of the charges Q λ 0 n , with λ 0 regular zeros and n ∈ E λ 0 . Since these charges are local, they are also in involution with the momentum P of the theory, whose Poisson bracket generates the derivative with respect to the spatial coordinate x. We have not yet discussed the conservation properties of these charges. For the models that we consider as examples in this article (listed in paragraph 2.1.4), we will see in section 7 that the Hamiltonian H of the theory always belongs to the algebra of local charges described above. It therefore follows that all these charges are conserved. More precisely, we will see that H is always a linear combination of the quadratic charges Q λ 0 2 and the momentum P.

Gauge invariance
In this subsection, we anticipate the application of the construction developed in the previous sections to integrable σ-models on Z T -coset spaces. In those models, infinity is a regular zero and the corresponding field C(x) (cf. subsection 2.2) is a gauge constraint. We prove here that the currents K λ 0 n (x) constructed at regular zeros λ 0 in the previous sections are gauge invariant, in the sense that they Poisson commute with the constraint C(y). As the K λ 0 n 's are polynomials of the J λ 0 n 's, it is enough to prove the following theorem.
Finally, let us treat the case λ 0 = ∞, for which J ∞ n (x) is given by the coefficient of α rn in T ∞ n (α, x) = Tr ψ(α) n L ∞ (α, x) n . Using the definition of L ∞ and the Poisson brackets (2.24) and (2.25), we find This bracket has the same structure as equation (2.24). Therefore, the case λ 0 = ∞ is treated exactly in the same way than the case λ 0 = 0, which ends the proof.

Reality conditions
To close this section let us discuss the reality conditions on the charges Q λ 0 n extracted at regular zeros in the previous sections. In the examples of paragraph 2.1.4, we consider integrable σmodels with target space G 0 or a quotient of G 0 , where G 0 is a real Lie group. If g 0 is the Lie algebra of G 0 , then the Lax matrix of the model is a g-valued field, where g is the complexification of g 0 . In other words, g 0 is a real form of the complex Lie algebra g. The fact that the σ-models we consider are on the real form G 0 (or one of its quotient) is encoded on the Lax matrix as the following reality condition: where the bar denotes complex conjugation and τ is the involutive semi-linear automorphism of g characterizing the real form g 0 . Moreover, we also have a reality condition on the twist function, which simply reads ϕ(λ) = ϕ(λ).

(5.5)
In particular, if λ 0 is a zero of ϕ, its conjugateλ 0 is also a zero of ϕ. Combining the reality conditions (5.4) and (5.5), we also see that if λ 0 is a regular zero (see paragraph 2.1.3),λ 0 is also a regular zero. Thus, the regular zeros can be of two types: real ones λ 0 ∈ R and conjugate pairs λ 0 ,λ 0 .
We will use the reality condition (5.4) in a similar way to the way we used the equivariance property (2.8) in subsection 4.1. In particular, as we consider powers of the Lax matrix, which are not in the Lie algebra g in general, we will need to extend "naturally" the automorphism τ to the whole algebra F of matrices acting on the defining representation of g. This was done for the automorphism σ in subsection 4.1 and appendix A. One can apply similar ideas to τ , using the classification of real forms of the classical Lie algebras A, B, C and D. We do not present the details here and just summarise the results.
There exists an extension of τ on the whole algebra of matrices F , which coincides with τ when restricted to the Lie algebra g, and that we shall still denote τ . This extension is still an involutive semi-linear map of F to itself. However, it is not in general an algebra homomorphism. The main properties of the extension τ that we will need are the following. There exists γ ∈ {1, −1} such that
Consider a regular zero λ 0 . Suppose first that λ 0 is complex: its conjugateλ 0 is then also a regular zero. According to the previous sections, we can extract two towers of (possibly complex) currents J λ 0 n (x) and Jλ 0 n (x) by evaluating T n (λ, x) at λ = λ 0 or λ =λ 0 (note that λ 0 cannot be a cyclotomic point as it is complex). However, according to equation (5.7), these currents are not independent. Indeed, they are related by the reality condition Thus, considering linear combination of Q λ 0 n and Qλ 0 n , we extract from each pair λ 0 ,λ 0 of complex regular zeros two towers of real charges in involution: Q λ 0 n + γ n Qλ 0 n and i Q λ 0 n − γ n Qλ 0 n .
Suppose now that λ 0 is a real and non-cyclotomic regular zero. Equation (5.7) then imposes the reality condition Thus, the current J λ 0 n is either real or pure imaginary. In each case, we can extract only one tower of real local charges. Consider now the case where λ 0 is the origin and thus a cyclotomic real point. The current J 0 n (x) is then the coefficient of λ rn in the power series expansion of T n (λ, x). Yet, this coefficient is also the one ofλ rn in the power series expansion of T n (λ, x). The reality condition (5.7) then implies that equation (5.8) also holds for λ 0 = 0.
Finally, let us discuss the case where λ 0 is infinity, which we consider as a real point. From the reality conditions (5.4) and (5.5), we find that the field C(x) defined in subsection 2.2 is real, in the sense that τ C(x) = C(x). We then obtain reality conditions on the Lax matrix L ∞ (α, x) and the twist function ψ(α) similar to equations (5.4) and (5.5). As a result we can apply the above discussion, since the point at infinity in the variable λ corresponds to the origin in the variable α, and conclude that equation (5.8) also holds for λ 0 = ∞.
To summarise this subsection, we have shown that one can extract: • one tower of real local charges for each real regular zero λ 0 , • two towers of real local charges for each pair λ 0 ,λ 0 of complex regular zeros.
In other words, one can extract as many towers of real charges as there are regular zeros.

Integrable hierarchies and zero curvature equations
In the previous sections, we constructed a infinite set of local charges Q λ 0 n in involution, with λ 0 regular zeros. It induces an infinite set of commuting Hamiltonian flows, defined by Q λ 0 n , · . In this section, we show that these flows generate a hierarchy of integrable equations. More precisely, we associate with each charge Q λ 0 n a connection which commutes with the connection ∇ x = ∂ x + L(λ, x). We show that the connections ∇ λ 0 n also commute with one another for finite regular zeros λ 0 . The commutativity of these connections takes the form of zero curvature equations. In particular, we will use the zero curvature equations involving L(λ, x) and the M λ 0 n (λ, x)'s to prove that the local charges Q λ 0 n are in involution with the non-local charges extracted from the monodromy of L(λ, x).

Zero curvature equations with L
The starting point of this article is an integrable system with Lax matrix L(λ, x) and Hamiltonian H. The dynamical equations of this system are generated by the Poisson bracket with H. They are encoded in the form of a zero curvature equation on the Lax matrix L(λ, x) and a g-valued matrix M(λ, x). In this subsection, we study the dynamics of the Lax matrix under the Hamiltonian flows generated by the local charges Q λ 0 n constructed in the previous sections. More precisely, we show that these dynamics also take the form of a zero curvature equation on L(λ, x): Theorem 6.1. Let λ 0 ∈ Z and n ∈ E λ 0 . There exists a matrix M λ 0 n (λ, x) such that we have the zero curvature equation Proof. Let us apply the second result of Corollary 2.3 to the r/s-system (2.1). Using the form (2.5) of the R-matrix, we find {L(λ, x), T n (µ, y)} = n Tr 2 R 0 12 (λ, µ)S n−1 (µ, y) 2 , L(λ, x) δ xy (6.2) − n Tr 2 R 0 12 (λ, µ)S n−1 (µ, y) 2 δ ′ xy − n ϕ(µ) ϕ(λ) Consider first the case where λ 0 is a non-cylotomic regular zero. Evaluating the equation above at µ = λ 0 and using ϕ(λ 0 ) = 0, we have Suppose now that λ 0 is the origin, which is a cyclotomic point, in which case J 0 n (y) is constructed as the coefficient of µ rn in the power series expansion of T n (µ, y). Moreover, as n ∈ E 0 , we have r n < T − 1 (see section 4). The Poisson bracket {L(λ, x), J 0 n (y)} is thus the µ rn -term in equation (6.2). We have ϕ(µ) = O(µ T −1 ) and r n < T − 1, thus the last term of equation (6.2) cannot contribute to µ rn . Thus, we also have equation (6.3) for λ 0 = 0, with We will say that N λ 0 n is the Lax matrix associated with the charge defined as the integral of the current J λ 0 n . Equation Using the fact that the Poisson bracket is a derivation, we find from equation (6.3) that After integration over y, we get the required zero curvature equation.
Thus, the Hamiltonian flows of the charges Q λ 0 n generate dynamical equations that can be recast in the form of zero curvature equations. In conclusion, we have constructed a hierarchy of integrable systems with Lax matrix L(λ, x) and Hamiltonians Q λ 0 n . The zero curvature equations of Theorem 6.1 can be seen as the commutativity of the connections with ∇ x = ∂ x +L(λ, x). This connection ∇ x can be thought of as the connection associated with the local momentum P of the theory. As already mentioned, we will see in section 7 that for the models we consider, the Hamiltonian is given by a linear combination H = λ 0 ∈Z a λ 0 Q λ 0 2 + bP of the quadratic charges Q λ 0 2 and the momentum P. Therefore, the matrix M(λ, x) of equation (6.1) can be constructed as λ 0 ∈Z a λ 0 M λ 0 2 (λ, x) + bL(λ, x).
Theorem 6.1 only treats the case of finite regular zeros λ 0 . Let us also briefly discuss what happens when λ 0 = ∞. In this case, J ∞ n (x) is extracted from the Lax matrix L ∞ (α, x). Since this matrix satisfies an r/s-system with twist function ψ(α), one can apply the method developed here. Doing so we find that the dynamics of L ∞ (α, x) under the Hamiltonian flow of Q ∞ n takes the form of a zero curvature equation. Moreover, starting with the Poisson bracket (5.2) and working weakly, we also find a weak curvature equation in the same way as M λ 0 n (λ, x) was built from N λ 0 n (λ, x) for a finite regular zero λ 0 . In other words, Theorem 6.1 also applies for λ 0 = ∞ when Poisson brackets are considered weakly.
Let us end this subsection by stating a few properties of the Lax matrix M λ 0 n (λ, x). Using the equivariance property (2.7), we find that The Lax matrix M λ 0 n thus satisfies the same equivariance property (2.8) as L. Recall that the Lax matrix N 0 n (λ, x) is extracted as the µ rn -term in Consider the equivariance properties (4.3) and Combining it with the fact that Tr σ(Y )σ(Z) = Tr(Y Z) for any matrices Y, Z ∈ F (see appendix A), we find that N n (ωµ ; λ, x) = ω rn N n (µ ; λ, x). (6.8) Therefore, the power series expansion of N n (µ ; λ, x) in µ contains powers of the form r n + kT , with k ∈ Z ≥0 . In particular, N 0 n (λ, x) is the coefficient of the smallest power in this expansion, in the same way as J 0 n (x) is in the expansion of T n (µ, x). Let us define M n (µ ; λ, x) from N n (µ ; λ, x) and T n (µ, x) in the same way we constructed M λ 0 n (λ, x) from N λ 0 n (λ, x) and J λ 0 n (x). In particular, M λ 0 n (λ, x) is the evaluation of M n (µ ; λ, x) at µ = λ 0 . From equations (4.6) and (6.8), we find the following equivariance property M n (ωµ ; λ, x) = ω rn M n (µ ; λ, x). (6.9) is the coefficient of the first term µ rn in the power series expansion of M n (µ ; λ, x).

Involution with non-local charges
In this subsection, we use the result of the previous one to prove that the local charges Q λ 0 n are in involution with the non-local charges extracted from the monodromy of the Lax matrix L(λ, x). This monodromy is defined as the path-ordered exponential where the integral is taken on the real line R or the circle S 1 , depending on the coordinate space of the model. Consider also the partial transfer matrices These matrices satisfy the initial condition T (λ ; x, x) = Id and the differential equations x, y), (6.10a) ∂ y T (λ ; x, y) = T (λ ; x, y)L(λ, y). (6.10b) Moreover, the variation of T under a infinitesimal variation δL of L is given by This formula allows one to compute derivatives of T and in particular its Poisson bracket with the local charge Q λ 0 n , for λ 0 ∈ Z and n ∈ E λ 0 . Specifically, we have The Poisson bracket of Q λ 0 n and L(λ, z) is given by Theorem 6.1. Using this together with the equations (6.10) and (6.11), we find x, y). If the spatial coordinate is taken on the real line (from −∞ to ∞) and the fields are assumed to be decreasing at infinity fast enough, we get Q λ 0 n , T (λ) = 0, i.e. the whole monodromy T (λ) is in involution with Q λ 0 n . If the spatial coordinate is taken on the circle (from 0 to 2π) and the fields are assumed to be periodic, we get Q λ 0 n , T (λ) = T (λ), M λ 0 n (λ, 0) . In this case, Q λ 0 n Poisson commutes with any central function of T (λ), e.g. the traces Tr T (λ) k and the determinant det T (λ) . Thus, we have Theorem 6.2. The monodromy T (λ) (resp. the central functions of T (λ)) is in involution with the local charges Q λ 0 n for λ 0 ∈ Z and n ∈ E λ 0 , if the spatial coordinate is taken on the real line (resp. the circle). In particular, it is conserved.
Proof. It just remains to prove the conservation of the non-local charges. This follows from the fact that the Hamiltonian H can be expressed as a linear combination of the quadratic charges Q λ 0 2 and the momentum P.
Once again, this theorem applies only for finite regular zeros λ 0 . Following a similar argument to the one given in the previous subsection, it also holds for the charges Q ∞ n if we consider Poisson brackets only weakly.

Zero curvature equations between the M λ 0 n 's
In subsection 6.1, we showed that the dynamics of the Lax matrix L(λ, x) under the Hamiltonian flow of the local charge Q λ 0 n takes the form of a zero curvature equation with a matrix M λ 0 n (λ, x). We thus exhibited a hierarchy of integrable equations, corresponding to the commutativity of the connections ∇ λ 0 n with ∇ x . This can be seen as the compatibility condition of the two auxiliary linear problems ∇ x Ψ = 0 and ∇ λ 0 n Ψ = 0, with Ψ a function on the phase space, valued in the connected and simply connected Lie group with Lie algebra g. In this subsection, we prove that the connections ∇ λ 0 n and ∇ µ 0 m also commute with one another (except when λ 0 is finite and µ 0 = ∞). This can be seen as the simultaneous compatibility of all auxiliary linear problems ∇ λ 0 n Ψ = 0 and it takes the form of zero curvature equations: Theorem 6.3. Let λ 0 , µ 0 ∈ Z, n ∈ E λ 0 and m ∈ E µ 0 . We have the zero curvature equation This subsection is entirely devoted to the proof of Theorem 6.3. After stating some general results, we will treat separately the cases λ 0 = µ 0 and λ 0 = µ 0 . Note that for the latter, we only have a complete proof for an algebra g of type B, C and D. For g of type A, we verified Theorem 6.3 for the first degrees n and m and conjecture that it holds more generally for any n and m. To improve the clarity of the subsection, some technical details of the proof are presented in appendix B.
Here also the theorem concerns the finite regular zeros λ 0 and µ 0 . The method presented in this subsection also applies for λ 0 = µ 0 = ∞ as L ∞ (λ, x) also satisfies an r/s-system with twist function (Theorem 2.1). However, the theorem does not hold when λ 0 is finite and µ 0 = ∞, even if Poisson brackets are considered only weakly.
The currents J k are extracted from T k . But in general, the charges are constructed from currents K k which are extracted from W k , where the definition of W k depends on g and σ (see subsections 3.5 and 4.7). In particular, we have K λ 0 k (x) = W k (λ 0 , x) for a non-cyclotomic regular zero λ 0 . For the origin, which is cyclotomic, Using the expression of W k and M k in terms of T k and N k , we see that Z λµ nm (ρ, x, y) contains several types of terms: and r k < T − 1, hence Ξ λµ 0 kl (ρ, x, y) λ r k = 0. Similarly Ξ µ 0 λ lk (ρ, x, y) λ r k = 0. Thus, the coefficient of λ rn in Z λµ 0 nm (ρ, x, y) vanishes, as required. This ends the proof of Theorem 6.3 for different regular zeros λ 0 and µ 0 .

Zero curvature equations at a non-cyclotomic regular zero
Let us now prove Theorem 6.3 for λ 0 = µ 0 . We start with the case where λ 0 is a non-cyclotomic point. We then want to show that dy Z λ 0 λ 0 nm (ρ, x, y) = 0.
As in section 3, we treat separately the Lie algebras of type B, C and D and the Lie algebras of type A. Suppose first that g is of type B, C or D. In this case, the currents K λ 0 2n are equal to the currents J λ 0 2n (see subsections 3.2 and 3.5) and the corresponding Lax matrices M λ 0 2n are equal to the matrices N λ 0 2n . Thus, Z λ 0 λ 0 2n 2m (ρ, x, y) is simply equal to Y λ 0 λ 0 2n 2m (ρ, x, y) (see paragraph 6.3.1). According to equation (6.15), we have where Ξ was defined in equation (6.14). To avoid cluttering the argument in the present paragraph with too many technicalities, we postpone the details of the computation of Ξ λ 0 λ 0 2n 2m (ρ, x, y) in appendix B.1. We find where the function f λ 0 2n 2m satisfies f λ 0 2n 2m = f λ 0 2m 2n (cf. appendix B.1). It then follows that from which we deduce that dy Y λ 0 λ 0 2n 2m (ρ, x, y) = 0, as required.
Suppose now that g is of type A. In this case, the currents K λ 0 n are different from the currents J λ 0 n and we therefore have to consider Z λ 0 λ 0 nm rather than simply Y λ 0 λ 0 nm . According to the discussion at the end of paragraph 6.3.1, it contains polynomials in the J λ 0 p 's and N λ 0 p 's, multiplied by either Ξ λ 0 λ 0 kl (ρ, x, y) or J λ 0 k (x), J λ 0 l (y) . This last Poisson bracket is given by equation (3.10) and is expressed in terms of the J λ 0 p 's. As for type B, C and D, we compute the expression of Ξ λ 0 λ 0 kl in appendix B.1. We find for some function f λ 0 kl such that f λ 0 kl = f λ 0 lk . Hence Z λ 0 λ 0 nm can be expressed in terms of the J λ 0 p 's and N λ 0 p 's, up to terms involving f λ 0 kl . The latter are always of the form with α a constant. Moreover, one can check that for any such term, there is also a similar one but with an opposite sign and k and l interchanged. Using the symmetry property f λ 0 kl = f λ 0 lk , one can then conclude that these terms always vanish.
Therefore Z λ 0 λ 0 nm can be expressed in terms of the J λ 0 p 's and N λ 0 p 's. Using the first explicit expressions (3.12) for the current K λ 0 n and the corresponding expressions for the matrices M λ 0 n , it can be check directly that dy Z λ 0 λ 0 nm (ρ, x, y) = 0 for the first few degrees n, m. Specifically, we have checked this for degrees n and m up to 7. In particular, we observed that we could not have chosen different coefficients in equation (3.12) for these zero curvature equations to hold (in the same way that these coefficients were uniquely fixed by requiring the involution of Q λ 0 n and Q λ 0 m ). Based on these strong observations, we conjecture that it holds for any n, m ∈ E λ 0 .

Zero curvature equations at a cyclotomic regular zero
Finally, let us prove Theorem 6.3 for λ 0 = µ 0 = 0, which is a cyclotomic point. Remember that K 0 n (x) and M 0 n (ρ, x) are extracted as the coefficient of λ rn in the power series expansion of W n (λ, x) and M n (λ ; ρ, x) where r n is the smallest power appearing in these expansions. That is, Theorem 6.3 for λ 0 = µ 0 = 0 is equivalent to the statement that dy Z λλ nm (ρ, x, y) λ rn+rm = 0.
The computation of Ξ λλ 2n 2m λ rn+rm is performed in appendix B.2. The final result is 2n 2m a function symmetric under the exchange of n and m. By virtue of this symmetry we find that the terms involving f disappear in Y λλ 2n 2m (ρ, x, y) λ rn+rm , while the other terms vanish when integrated over y, as required.
Consider now g = sl(d, C) of type A. The construction of the currents K 0 k depends on σ being inner or not (see subsections 4.4, 4.5 and 4.7). If σ is inner, then the currents K 0 k and W k are equal to the currents J 0 k and T k . In this case, we have The expression for Ξ λλ nm (ρ, x, y) λ rn+rm is given by equation (B.7) of appendix B.2. It has the same structure as in the case of types B, C and D: the same arguments then apply and we conclude that the integration of Y λλ nm (ρ, x, y) λ rn+rm over y vanishes.
Finally, consider g = sl(d, C) of type A with σ not inner. In this case, the currents K 0 k (x) and W k (λ, x) are constructed as polynomials of respectively J 0 k (x) and T k (λ, x). The corresponding structure of Z λµ nm is discussed at the end of paragraph 6.3.1. In particular, Z λλ nm (ρ, x, y) is composed of two types of terms: • Ξ λλ kl (ρ, x, y), multiplied by polynomials in the T j (λ, ·)'s, • T k (λ, x), T l (λ, y) , multiplied by polynomials in the T j (λ, ·)'s and N j (λ ; ρ, ·)'s.
We want to extract the coefficient of λ rn+rm in Z λλ nm (ρ, x, y). Recall that the powers of λ appearing in T j (λ, ·) and N j (λ ; ρ, ·) are of the form r j + aT , with a ∈ Z ≥0 , and that the coefficients corresponding to a = 0 are J 0 j (·) and N 0 j (ρ, ·). In the two types of terms mentioned above, one can check that the terms with a > 0 will not contribute to the coefficient of λ rn+rm . More precisely, Z λλ nm (ρ, x, y) λ rn+rm is composed of two types of terms: • Ξ λλ kl (ρ, x, y) λ r k +r l , multiplied by polynomials in the J 0 j (·)'s, • J 0 k (x), J 0 l (y) , multiplied by polynomials in the J . j (·)'s and N 0 j (ρ, ·)'s.
The Poisson brackets J 0 k (x), J 0 l (y) are given by equation (4.16). The expression for Ξ λλ kl (ρ, x, y) λ r k +r l is worked out in appendix B.2 and reads Ξ λλ kl (ρ, x, y) kl is a function invariant under the interchange of k and l. The rest of the argument follows closely that given in the non-cyclotomic case. Specifically, the terms containing f (0) kl are seen to vanish by virtue of this symmetry property. We can thus express Z λλ nm (ρ, x, y) λ rn+rm in terms of the J 0 k 's and N 0 k 's only. One then can check explicitly that this expression vanishes when integrated over y, as required. We verified this for the first few degrees n and m (up to 8) and different values of T (from 2 to 6). We therefore conjecture that this is also true for any n, m ∈ E 0 and any T .

Applications
In paragraph 2.1.4, we gave a list of integrable σ-models which fit the framework of the present article. In this section, we apply the methods developed in the previous sections to these particular examples, analyse the results and compare them to some existing work in the literature. These models were recently re-interpreted as particular examples of so-called dihedral affine Gaudin models [8]. We explain in the last part of this section how the framework of dihedral affine Gaudin models is particularly suited to apply the methods of the present article.

Principal chiral model and its deformations
Let us start with the simplest integrable σ-model, the Principal Chiral Model (PCM). The study of local charges of the PCM is already well known and was treated in the reference [44]: these results were the principal motivation and guideline for the present article. In particular, one of the aims was to generalise the construction of [44] to a wider class of models, among which are the integrable two-parameters deformations of the PCM (dPCM). We shall discuss the latter in this subsection.
The integrable structure of the dPCM was discussed in subsection 2.1.4. In this case, σ = Id so that T = 1. Their Lax matrix and twist function are given by equations (2.13) and (2.15) respectively. In the language of this article, the regular zeros of these deformed models are +1 and −1. The evaluation of ϕ(λ)L(λ, x) at these zeros gives the fields J ± (x) of equation (2.17).
The local charges Q ±1 n constructed in the present article are related to the traces of powers of these fields. In the undeformed case, i.e. for the PCM, these fields coincide with the fields j ± = −g −1 ∂ ± g used in reference [44] to construct the local charges. Thus, the method presented in this article gives back the results of [44] for the PCM, as expected. In the deformed case, it generalises these results, while keeping a similar structure in the construction: in particular, we obtain two towers of local charges in involution, corresponding to the two chiralities of the model, and the spin of these charges is still related to the exponents of the affine Kac-Moody algebra g, as it was in the PCM case [44].
The Hamiltonian and momentum of the deformed PCM are given by equations (2.16). One can check that these are related to the quadratic charges Q ±1 2 as follows , .
In particular, the Hamiltonian belongs to the algebra of local charges in involution so that these charges are conserved (see also the discussion at the end of subsection 5.1).
This observation also allows one to recover the temporal component M(λ, x) of the Lax pair of the model (see the paragraph below equation (6.6)). More precisely, the equation of motion of the dPCM can be recast as the Lax equation (6.1), where This zero curvature equation (6.1) is the first among a whole hierarchy of integrable equations generated by the local charges Q ±1 n (cf. subsection 6.1).
In particular, this result was used in subsection 6.2 to show that the local charges Q ±1 n are in involution with the non-local charges extracted from the monodromy of the Lax matrix L(λ, x) (see Theorem 6.2). In [44], it was shown that the local charges of the undeformed PCM Poisson commute with the non-local charges generating the classical Yangian symmetry of the model. In the framework of this article, if we consider the model on the real line R, we expect these non-local charges to be extracted from the expansion of (a gauge transformation of) the monodromy around the pole λ = 0 of the twist function of the PCM. For Yang-Baxter deformations (k = 0 and η = 0, see paragraph 2.1.4), this Yangian symmetry gets deformed to a quantum affine symmetry [42,43]. In particular, studying the monodromy around the poles ±iη of the twist function of the Yang-Baxter model, one can extract a q-deformed affine Poisson-Hopf algebra U q ( g). We have therefore proved that this algebra of non-local charges is in involution with the algebra of local charges consisting of the Q ±1 n 's.
As mentioned in the paragraph 2.1.4, the PCM and its deformation are defined on a real Lie group G 0 , whose Lie algebra g 0 is a real form of g. This real form is characterised by a semi-linear involutive automorphism τ . The Lax matrix (2.13) of these models satisfies the reality condition (5.4). Moreover, the twist function (2.15) verifies the reality condition (5.5) and the regular zeros of the model (+1 and −1) are real. Thus, the discussion of the subsection 5.3 applies to these models and the charges Q ±1 n are real (possibly up to a redefinition of some Q ±1 n by a factor of i, depending on τ ).

Bi-Yang-Baxter σ-model
There exists another two-parameter deformation of the PCM, the so-called bi-Yang-Baxter (bYB) σ-model [17,29]. Its Hamiltonian integrability was established in [30], by viewing it as a deformation of the Z 2 -coset model G 0 × G 0 /G 0,diag . In this formulation, the bYB model falls into the framework of the present article but has to be considered as a model with a gauge constraint (this is also the approach of the reference [8]). This can be treated with the methods developed here, using the results of subsections 2.2 and 5.2 to take into account the gauge constraint. The analysis of the local charges for this model would then be close to the one for the Z 2 -coset model that we perform in the next subsection.
In this subsection, we choose to treat the bYB model in its gauge fixed formulation, which is more natural if we see it as a deformation of the PCM. As we will see, this formulation does not fit exactly within the framework of this article, but we will explain how one can overcome this difficulty by further relaxing the general assumptions we made. Even though this generalisation could have been done throughout the entire article, we chose here to work in a more restricting but more common framework for clarity and simplicity. In this regard, the present subsection is also used to illustrate, a posteriori, how the methods and results we found apply under the generalised conditions.
The Hamiltonian integrability of the bYB model in its gauge fixed formulation was studied in the last section of the article [30]. In order to see this model as a deformation of the PCM, one should first perform a change of variables from the conventions of [30]: more precisely, we relate the spectral parameter z of [30] with the spectral parameter λ of this article by the involutive transformation z = 1 − λ 1 + λ .
In particular, after this change of spectral parameter, the Lax matrix of the gauge fixed-model takes the form for some g-valued fields j 0 , j 1 and j ∞ . In particular, the field j ∞ vanishes when the deformation parameters η andη go to zero. We then recover the Lax matrix (2.13) of the undeformed PCM. In this sense, the gauge fixed model with spectral parameter λ describes an integrable deformation of the PCM.
The equation (3.8) of [30] gives the twist function of the bYB model in terms of the spectral parameter z, which we shall denote χ bYB (z) here. We define a corresponding function ϕ bYB (λ) of the spectral parameter λ by the relation It is found to be of the form where a, b, c and d depend on the deformation parameters η andη. This is to be compared with the twist function (2.15) (where k = A = 0) of the PCM. The presence of a global factor in the twist function is directly related to the global factor K in the definition of the action of the model: it can be reabsorbed by setting K to a particular value.
In the undeformed case η =η = 0, one has b = c = 0. Thus, up to the global factor, the function ϕ bYB (λ) coincides with the twist function of the PCM in this case. The bi-Yang-Baxter deformation thus has for effect to deform the poles of the twist function of the PCM: the two double poles at 0 and ∞ in ϕ PCM (λ) dλ split into four simple poles on the imaginary axis. The zeros of the twist function stay undeformed in this procedure: indeed, the zeros of ϕ bYB are +1 and −1, as for the PCM and its deformations. Considering the expression (7.2) of the Lax matrix of the bYB model, we see that these zeros are regular.
It was shown in [30] that the Lax matrix of the gauge fixed model satisfies a non-ultralocal Poisson bracket of the form (2.1). Using the parameters of [30] and after performing the change of spectral parameter described above, we find that the R-matrix describing this Poisson bracket takes the form In the above equation, R is the non-split solution of the modified CYBE considered in [30] and we use the standard R 0 -matrix In the undeformed case,η and c vanishes and we recover the integrable structure of the PCM, as described in paragraph 2.1.4. However, in general, the matrix R bYB is different from R 0 . Thus, the bYB model does not fit exactly within the framework described in subsection 2.1. Let us note that R bYB is still a solution of the classical Yang-Baxter equation (CYBE) (2.2), using the fact that R is a non-split solution of the modified classical Yang-Baxter equation (see [30]) and the identity C ij , X i + X j = 0, true for any X ∈ g. Note that the coefficients of C 12 and R 12 in equation (7.3) are strongly constrained by the requirement that R bYB fulfils the CYBE.
As explained above, R bYB is different from R 0 and therefore we cannot directly apply the results of this article. However, going through the details of the proofs of these results for a non-constrained model with T = 1, we see that the only properties of the matrix R 0 that we used are the CYBE (for the zero curvature equations), the fact that R 0 12 (λ, µ) is holomorphic at pairs (λ 0 , µ 0 ) of distinct regular zeros and the asymptotic property (3.2) near a regular zero µ = λ 0 . The matrix R bYB also satisfies the CYBE, as explained above. Moreover, one easily checks that it also verifies the holomorphy condition and the asymptotic property mentioned above. Thus, the results we found in this article also apply to the bYB model. This is a general observation: we can also treat the models where the matrix R 0 is replaced by a matrix R ′ satisfying some similar properties. More precisely, we require that R ′ • obeys the CYBE (2.2), • is holomorphic at (λ 0 , µ 0 ) with λ 0 and µ 0 different regular zeros in Z, • verifies the asymptotic property (3.2) around non-cyclotomic regular zeros, • satisfies the equation (4.8) for U 12 (λ, λ) around a cyclotomic regular zero, up to a term O(λ 2T −3 ) (which would not contribute to some (r n + r m ) th power of λ in (4.9)).
In particular, let us consider a matrix R ′ of the form like the matrix R bYB . Then R ′ 12 (λ, µ) is holomorphic for λ and µ going to different regular zeros if D is holomorphic at any regular zero (this is for example the case for the bYB model). This condition also ensures that the asymptotic property (3.2) is satisfied by R ′ . In the same way the condition on Let us note, however, that these conditions do not allow to treat the case where infinity is a regular zero in the same way that we did in this article (subsections 2.2 and 4). This would require, among other conditions, that the asymptotic properties (2.23) at infinity are also satisfied by the matrix R ′ . One can check that a matrix R ′ of the form (7.4) can never satisfy the second property of equation (2.23).
As explained above, we can apply the construction of local charges in involution to the bYB model. These local charges will be very similar to the ones of the PCM and its deformations, described in the previous subsection, so we shall not enter into much details here. Let us note that these charges are related to traces of powers of j 0 (x) ± j 1 (x), where j 0 and j 1 are the fields appearing in the Lax matrix (7.2). As in the case of the PCM (see previous subsection), the Hamiltonian and the momentum of the bYB model are related to the quadratic charges Q ±1 2 by the relation .
In particular, the local charges constructed above are all conserved.

Z T -coset models and their deformations
In this subsection, we discuss the construction of local charges in involution for Z T -coset models (and the deformations of Z 2 -coset models). These models were described in paragraph 2.1.4. The order T of σ is strictly greater than one. Their twist function and Lax matrix are given by equations (2.20) and (2.19). Such local charges were constructed for symmetric spaces, i.e. Z 2 -cosets, in references [46] and [47]: we shall compare these results with the ones of this article.
As already mentioned in the paragraph 2.1.4, the regular zeros of these models are the origin and infinity, which are both cyclotomic points. We shall therefore apply here the construction of section 4. Moreover, all these models possess a gauge constraint Π (0) (see paragraph 2.1.4), which is identified with the field at infinity C(x) described in subsection 2.2. The results of subsection 5.2 ensure that the densities of the local charges that we construct here are gauge invariant. Indeed, by Theorem 5.3, these densities Poisson commute with the constraint C.
As in the case of the PCM (see subsection 7.1), the degrees of the local charges are related to the exponents of the affine Kac-Moody algebra g plus one (here also, we do not consider the exponents corresponding to the Pfaffian for type D). However, as explained in section 4, the fact that the regular zeros of the model are cyclotomic makes some of the exponents 'drop out', in the sense that we cannot construct a charge of the corresponding degree. Recall that a degree n (corresponding to an exponent n − 1) drops out if r n is equal to T − 1 (where r n was defined in subsection 4.1).
Let us study this in more detail for the case of Z 2 -cosets. In particular, we shall compare this phenomenon of exponents dropping out with some results of reference [46]. Indeed, in this reference, some local charges in involution were constructed for symmetric spaces (i.e. Z 2cosets). These symmetric spaces correspond to quotients G 0 /G σ 0 of the real Lie group G 0 by the subgroup of fixed points under the involutive automorphism σ (see paragraph 2.1.4). Such spaces were classified, up to isomorphism, for classical compact groups G 0 .
In particular, the possible exponents (i.e. the degrees minus one) of the local charges for each symmetric space of this classification were listed in Table 1 of [46]: they form a (potentially proper) subset of the exponents of g. A simple case by case computation of the integers r n for these symmetric spaces, and thus these automorphisms σ, shows that the exponents of g which do not appear in this list are exactly the exponents that drop out in the formalism of the present article. We therefore recover the structure of the degrees of local charges found in [46] (except for the integer h of [46], which we could not interpret in the present formalism).
An explicit computation of the traces of powers of ϕ Z 2 (λ)L Z 2 (λ, x) around the origin λ = 0 shows that the charges constructed in this article coincide, up to some factors, with the ones constructed in reference [46]. The two regular zeros 0 and ∞ correspond to the two chiralities of the model. The article [46] focused on one particular chirality. Here, we also have the Poisson brackets between the two towers of local charges constructed in this way. Indeed, according to Theorem 5.2, we show that these two towers of charges Poisson commute weakly.
This article also generalises the results of [46] in different directions. First of all, the present formalism also allows to treat the integrable deformations of the Z 2 -coset model. Indeed, as explained in paragraph 2.1.4, the regular zeros of these models are still 0 and ∞ and so the methods developed here still apply. The main generalisation is that this article does not restrict to (compact) symmetric spaces and also generalises the construction to any Z T -coset model. Finally, in this article we have also studied the hierarchy of equations induced by the flow of these local charges.
The reference [47] deals with the local charges in involution for the supersymmetric models on symmetric spaces, working with a superalgebra g. Although we did not consider such models in this article, we expect the construction to extend to these theories, by working with the Grassmann envelope of g and replacing all traces by supertraces. One should however be careful about how the automorphism σ is extented to the whole matrix algebra (see appendix A) in these supersymmetric cases. Such considerations could allow the construction of local charges in involution for supersymmetric σ-models whose target space includes AdS manifolds, with possible applications to the hybrid formulations of string theory.
We end this subsection by observing that the Hamiltonian of the Z T -coset model is related to the quadratic charge Q 0 2 at the origin and the momentum P Z T of the theory by where ζ was defined in equation (2.10) (note that this expression also holds for the deformed Z 2 model). Thus, we conclude that the local charges constructed above are conserved, as they commute (at least weakly) with the Hamiltonian.
As the Z T -coset models are constrained models, their Hamiltonian is defined up to a term Tr µ(x)C(x) , where µ is a g-valued Lagrange multiplier. In this sense, H Z T defined above is a particular choice of such a Hamiltonian, which generates a strong zero curvature equation (6.1) on the Lax matrix L. Another choice of Hamiltonian involves the quadratic charge Q ∞ 2 extracted at infinity, namely where ζ ∞ is defined in the same way than ζ by ζ ∞ (α T ) = αψ(α). This Hamiltonian is weakly equal to H Z T and generates a strong curvature equation on the Lax matrix L ∞ .

Dihedral affine Gaudin models
Let us end this section by discussing briefly dihedral affine Gaudin models (DAGM) and their relation to the present article. These models were defined and studied recently in reference [8], extending the notion of cyclotomic Gaudin models of [50]. In particular, it was shown that all the integrable σ-models mentioned in the previous subsections and paragraph 2.1.4 are examples of such DAGM. In fact, we will argue in this subsection that the construction of local charges of the present article automatically applies to a broader class of DAGM than just these integrable σ-models. We use the notations of [8] and refer the interested reader to the original article for the details.
Reference [8] defines a DAGM as a Hamiltonian field theory with a Poisson algebra of local observables S ℓ ( g D C ), depending on the following data: • a semi-simple Lie algebra g ; • an automorphism σ of g of order T and a semi-linear involutive automorphism τ , forming the dihedral group Π = D 2T ; • a set of points z ⊂ P 1 containing infinity (not to be confused with the set of regular zeros Z of this article), associated with positive integers n x (x ∈ z), encoded as a divisor D = x∈z n x x ; • complex numbers ℓ x p (p = 0, . . . , n x − 1) for all x ∈ z \{∞} and ℓ ∞ q (q = 1, . . . , n ∞ − 1), called the levels of the DAGM.
The dihedral group Π acts on P 1 via the multiplication by ω and complex conjugation. The set z is chosen such that any two points of z are in disjoints orbits under the action of Π. We define the Gaudin poles Πz as the set of all images of points in z by Π. For this article, we shall restrict the discussion to the case where the highest levels ℓ x nx−1 are all non-zero (this is the main case treated in article [8]).
The twist function ϕ and Lax matrix L of the DAGM are introduced in the form of a "natural" connection ∇ = ϕ(λ)∂ x + ϕ(λ)L(λ, x). In particular, the twist function is uniquely defined by the levels ℓ x p of the model (see equation (4.44) of [8] and below). The Lax matrix L(λ, x) then naturally satisfies an r/s-type Poisson bracket (2.1) with twist function ϕ(λ). Moreover, the Lax matrix and the twist function of the DAGM are shown in [8] to verify the equivariance properties (2.8) and (2.9) and the reality conditions (5.4) and (5.5). Hence the DAGM automatically fits in the framework of the present article, described in subsection 2.1.
Let us now look at the zeros of the twist function. In the construction of [8], the connection ∇ defined above (and so in particular the twist function) has poles exactly at the Gaudin poles Πz. The zeros of the twist function then do not belong to Πz. Yet, the matrix component S(λ, x) = ϕ(λ)L(λ, x) of ∇ has poles only at the points Πz. Thus, all zeros of the twist function are regular zeros. The methods developed in this article thus apply naturally to the DAGM.
Finally, let us relate the discussion of the point at infinity in the present article, based on the results of subsection 2.2, to its treatment in a DAGM given in [8]. We will argue that the notions studied in that subsection can naturally be transposed to the DAGM setting. For that, we restrict attention to the case where n ∞ = 1. As one can check from [8], this hypothesis is equivalent to the condition that P (α, x), defined in equation (2.21), is regular at α = 0. With this condition fulfilled, we can then define the field C(x) as in equation (2.22). Throughout the present article, based on the examples of Z T -coset models, this field C(x) was seen as a gauge constraint (with Theorem 5.3 stating that the local charges Q λ 0 n are gauge invariant). This interpretation generalises to any DAGM, as the field C(x) computed above can be seen to coincide with the constraint defined in the paragraph 4.5.3 of the reference [8].
In order to define a constrained DAGM, it is required in [8] that the levels ℓ x 0 satisfy some additional relation (4.64). In the context of this article, this relation ensures the regularity of the twist function ψ(α), defined in equation (2.11), at α = 0. If we suppose that the DAGM is such that T > 1, then this regularity, combined with the equivariance property (2.31), implies that ψ(0) = 0. In this case, the point at infinity is then a regular zero of the model and the methods developed in the present article apply.

Outlook
In the present paper we showed how, in any integrable field theory described by an r/s-system with twist function ϕ(λ) and underlying finite-dimensional Lie algebra g of classical type, one can assign an infinite tower of local conserved charges Q λ 0 n , n ∈ E λ 0 , in involution to each regular zero λ 0 of ϕ(λ). If the regular zero λ 0 is non-cyclotomic then the set of degrees E λ 0 of these charges consists of one plus the exponents of the untwisted affine Kac-Moody algebra g associated with g (excluding the Pfaffian in the case of type D). This generalises the results established in [44,45] for the principal chiral model, with or without a Wess-Zumino term, on a real Lie group G 0 of classical type. On the other hand, if λ 0 is cyclotomic then some exponents may 'drop' and the set of degrees E λ 0 forms a subset of the set of exponents of g plus one. In the case where the order of the cyclotomy is T = 2, a simple computation reveals that the 'drops' occur precisely when g is of type A and the automorphism σ is inner. In particular, by performing a direct comparison with the results of [46] we find perfect agreement with the set of exponents naturally associated with a compact symmetric space G 0 /G σ 0 listed in Table 1 of that paper (see also [51]) describing the degrees of the local charges in involution in a symmetric space σ-model. We were, however, unable to interpret the Coxeter number of the compact symmetric space G 0 /G σ 0 in our present formalism. It would be interesting to understand whether there is a more fundamental connection between these two descriptions of the degrees of local charges in involution for symmetric space σ-models. This would also provide insight into the extension of the definition of the exponents to the case of Z T -cosets.
As noted in subsection 7.4, the regularity assumption on the zeros of the twist function is intimately related to a condition imposed on the levels of a dihedral affine Gaudin model in [8], so that the results of the present paper automatically apply to r/s-systems arising in the latter framework. In light of this fact, a noteworthy feature of the densities of the local charges Q λ 0 n is that they were all obtained from the particular combination ϕ(λ)L(λ, x) of the twist function and Lax matrix of the model. Indeed, the significance of this expression is that it forms the matrix component of the connection ∇ = ϕ(λ)∂ x + ϕ(λ)L(λ, x) characterising the dihedral affine Gaudin model. As argued in [51], the invariant tensors used to build local charges in involution in the principal chiral model and symmetric space σ-models can be obtained directly from the Weyl-invariant tensors appearing in the leading derivative-free terms of the local charges in involution of g-mKdV theory, given by the Drinfel'd-Sokolov construction. On the other hand, it follows from the results of [8] that g-mKdV theory is itself a cyclotomic affine Gaudin model whose corresponding connection ∇ is the starting point in the construction of Drinfel'd-Sokolov. It is therefore natural to speculate that a more general tower of local charges in involution, involving derivatives of the fields appearing in the Lax matrix L(λ, x), may be constructed starting from the dihedral affine Gaudin model connection ∇ itself.
With the application to dihedral affine Gaudin models in mind, it is also interesting to note that the results of the present paper bear a strong resemblance to the structure of commuting integrals of motion in a (quantum) Gaudin model for a finite-dimensional Lie algebra g. Indeed, it is well known from the work of Feigin, Frenkel and Reshetikhin [52] (see also [53,54]) that the algebra of commuting integrals of motion in a Gaudin model with N ∈ Z ≥1 sites, known as the Gaudin algebra, is generated by elements of the algebra of quantum observables U(g) ⊗N whose degrees are given by exponents of g plus one. An important open problem is to obtain an analogue of this result in the case of quantum Gaudin models associated with affine Kac-Moody algebras. Although this is still currently out of reach, the results of the present article may furnish useful hints in this direction by providing a similar description of local Poisson commuting integrals of motion in classical (dihedral) affine Gaudin models.
A Extension of the automorphism to the whole space of matrices We consider a Lie algebra g of classical type A, B, C or D, in its defining matricial representation. We therefore regard elements of g as acting linearly on a vector space V , i.e. as element of the space F of endomorphisms of V . Note that we can consider some connected matrix group G ⊂ F whose Lie algebra is g. Let σ be an automorphism of g, of finite order T . In this article, we are considering powers of elements of g, which do not belong to g in general but are elements of F . Thus, we want to extend the automorphism σ to the whole space of matrices F , in a "natural way".

A.1 The conjugacy case
Let us begin with the case where σ is inner, i.e. when σ : X ∈ g → QXQ −1 for some Q ∈ G.
Then the extension of σ to F , which by a slight abuse of notation we still denote as σ, can be naturally defined as This covers the case of types B and C, as they do not have any non trivial diagram automorphism.
Let us now consider the algebra D n , i.e. g = so(2n, C), for n ≥ 5. In this case, there always exists one non trivial diagram automorphism. However, this automorphism can be realised on the defining representation as an external conjugation: σ : X ∈ g → QXQ −1 where Q is not in the group SO(2n, C) but belongs to O(2n, C). In this case, the endomorphism σ as defined in equation (A.1) still naturally extends σ on F .
Let us say a few words on the algebra D 4 = so(8, C). It is known to have 6 diagram automorphisms, forming the triality, isomorphic to the symmetric group S 3 . One of them, of order 2, can also be realised as conjugation by a matrix Q ∈ O(8, C) and so extends to F = M 8 (C) by equation (A.1). The other non-trivial ones cannot be realised in any "natural way" on the defining representation and thus cannot easily be extended to the whole space M 8 (C). We shall not consider them in this article.
We now describe the properties of the extension σ defined in (A.1). The fact that the automorphism σ of g is of order T is equivalent to the fact that Q T belongs to the centraliser of g in F , By Schur's lemma, this implies that Q T = λ Id for some λ ∈ C. Therefore σ on F defined by (A.1) is also of order T . We shall make extensive use of the following five obvious properties of σ as defined in (A.1): σ(XY ) = σ(X)σ(Y ), σ(X n ) = σ(X) n , σ(Id) = Id, Tr σ(X) = Tr(X), Tr σ(X)σ(Y ) = Tr(XY ), for any X, Y in F .

A.2 The transpose case
The last case that we have to treat is the one of a type A algebra, with σ being not inner. We thus consider the defining representation g = sl(d, C). The action of σ on g can then always be expressed as σ : X ∈ g → −QX T Q −1 , where X T is the transpose of X and Q is a matrix in SL(d, C). Here also we can naturally extend σ to an endomorphism of F = M d (C), which we still denote σ, by letting σ : F −→ F X −→ −Q X T Q −1 . (A.2) Once again, let us investigate the properties of σ. As the automorphism σ of g is not inner, its order T must be even, and we shall write T = 2S. We note that σ 2 acts as conjugation by R = Q(Q T ) −1 . The fact that σ T = (σ 2 ) S = Id| g is thus equivalent to the fact that R S belongs to the centraliser Z F (g). Thus R S = λ Id for some λ ∈ C and so σ defined in (A.2) is also of order T . We end the subsection by noting the following five properties of σ: σ(X n ) = (−1) n−1 σ(X) n , σ(Id) = −Id, Tr σ(X) = −Tr(X), Tr σ(X)σ(Y ) = Tr(XY ), for any X, Y in F .

B Computation of Ξ
In this appendix, we give the details of the computation in some particular cases of the term Ξ λµ nm (ρ, x, y), defined by (6.14).

B.1 At a non-cyclotomic regular zero
We first suppose that λ 0 is a non-cyclotomic regular zero and that g is of type B, C or D.
In particular, note that f λ 0 nm = f λ 0 mn .

B.2 Around a cyclotomic regular zero
This subsection is devoted to the computation of Ξ λλ nm (ρ, x, y) around the origin λ = 0 and more precisely to the computation of the coefficient of λ rn+rm in its series expansion. Our starting point is the definition (6.14) of Ξ λµ nm . To evaluate this equation at µ = λ, we will need the expression of U 12 (λ, λ). We saw in section 4 that this is given by equation (4.8).
The presence of the partial Casimir C (0) 12 in this equation will gives rise to projections of S m−1 (λ, y) onto the grading F (0) = {Z ∈ F | σ(Z) = Z} of the matrix algebra. More precisely, the calculations will involve g λ mn (ρ, x, y) = nmλ −2 ζ(λ T )Tr 2 R 0 12 (ρ, λ) where S (0) p denotes the projection of S p on F (0) . As we are computing the coefficient of λ rn+rm in Ξ λλ nm , we will consider the λ rn+rm -term of g λ nm . Let us show that this term is actually always zero. Using the conventions and results of the subsections 4.1 and 4.2, in particular the integers α and q m , we see that the smallest power of λ appearing in g λ nm is a = αT − 2 + q m , as the S p (λ, x)'s are regular at λ = 0. Recall that r n and r m are both strictly less than T − 1 when n, m ∈ E 0 and that in this case, we have q m = r m + 1 (see subsection 4.2). Thus a = αT −1 + r m and hence a > r n + r m since α ≥ 1 and T − 1 > r n . We can then conclude that the coefficient of λ rn+rm in g λ nm vanishes, as announced.
We will also need the function Tr 2 R 0 12 (ρ, λ)S k+l (λ, x) 2 ∂ x S(λ, x) 2 S n+m−4−k−l (λ, x) 2 , similar to the function f λ 0 nm defined in the non-cyclotomic case (see previous subsection) and which possesses the same symmetry property f λ nm = f λ mn . As for g λ nm , we will use more precisely the function f (0) nm = f λ nm λ rn+rm , which is also symmetric under the exchange of n and m.
We now compute the coefficient of λ rn+rm in Ξ λλ nm . As explained at the beginning of this subsection, g λ nm does not contribute to this term and the contribution of f λ nm is defined as f nm . The contribution from the third term is calculated as in the case of types B, C and D. In particular, it vanishes when r n + r m is strictly less than T − 2.
Finally, let us discuss the contribution of the last term. First of all, we note that if a = 1, λ T ζ ′ (λ T ) − aζ(λ T ) = O(λ 2T ). Thus the powers of λ in this term are greater than 2T − 2. Yet, we have r n + r m < 2T − 2 for n, m ∈ E 0 , so this term does not contribute to the λ rn+rm -term in this case. Hence for σ inner, we have Ξ λλ nm (ρ, x, y) as in the case of types B, C and D. Suppose now that σ is not inner, so that a = 0. Then the smallest power of λ in the last term of equation (B.6) is greater than or equal to T − 2 + r n−1 + r m−1 . In subsection 4.5, we have shown that this is equal to r n + r m if both r n and r m are greater than S − 1 and that it is strictly greater than r n + r m otherwise. In conclusion, we find Ξ λλ nm (ρ, x, y)