Integrable deformations of coupled sigma-models

We construct integrability-preserving deformations of the integrable $\sigma$-model coupling together $N$ copies of the Principal Chiral Model. These deformed theories are obtained using the formalism of affine Gaudin models, by applying various combinations of Yang-Baxter and $\lambda$-deformations to the different copies of the undeformed model. We describe these models both in the Hamiltonian and Lagrangian formulation and give explicit expressions of their action and Lax pair. In particular, we recover through this construction various integrable $\lambda$-deformed models previously introduced in the literature. Finally, we discuss the relation of the present work with the semi-homolomorphic four-dimensional Chern-Simons theory.


Introduction and summary of the results
Integrable non-linear σ-models form an important class of two-dimensional classical integrable field theories. Their study was initiated more than 40 years ago and has found applications in various domains of physics, such as the AdS/CFT correspondence (see for instance the review [1]) and condensed matter theory [2]. A prototypical example of integrable σ-model is given by the Principal Chiral Model on a real semi-simple Lie group G 0 , with or without Wess-Zumino term. It describes the dynamics of a G 0 -valued field g(x + , x − ), where x ± = (t ± x)/2 denote the two-dimensional light-cone coordinates. Let g 0 be the Lie algebra of G 0 and κ the opposite of its Killing form. The action of this model is then given by where ρ and are constant parameters, ∂ ± denote the derivatives with respect to x ± and I WZ [g] is the Wess-Zumino term of g. The integrability of this model relies on the fact that its equation of motion can be recast in the form of a zero curvature equation ∂ + L − (z) − ∂ − L + (z) + L + (z), L − (z) = 0 on a Lax pair L ± (z). This Lax pair is valued in the complexification of g 0 and depends on an auxiliary complex parameter z, called the spectral parameter. It was shown by Klimčík in [3,4] that the Principal Chiral Model (without Wess-Zumino term) admits a continuous integrable deformation, called the Yang-Baxter model, which generalises to an arbitrary group G 0 a model constructed in [5] for the group SU (2). This deformed model depends on the choice of a skew-symmetric R-matrix on g 0 , i.e. a linear operator R : g 0 → g 0 satisfying the modified classical Yang-Baxter equation [RX, RY ] − R[RX, Y ] − R[X, RY ] = −c 2 [X, Y ] for every X, Y ∈ g 0 , with c equal to 1 or i. The action of the Yang-Baxter model is given by where η is the deformation parameter and R g = Ad −1 g •R•Ad g . This construction was later extended in various directions. For instance, one can construct integrable Yang-Baxter deformations of symmetric space σ-models [6], of superstrings on semi-symmetric spaces [7,8] and of the Principal Chiral Model with Wess-Zumino term [9]. Alternatively, one can also consider deformations based on homogeneous R-matrices [10], satisfying the non-modified (c = 0) classical Yang-Baxter equation.
Another type of integrable deformed σ-model, called the λ-model, was constructed by Sfetsos in [11]. It corresponds to a deformation of the non-abelian T-dual of the Principal Chiral Model (without Wess-Zumino term) and generalises a result obtained in [12] for the group SU (2). Its action is defined as where and λ are constant parameters and S WZW, [g] is the action of the conformal Wess-Zumino-Witten model at level (i.e. the action (1.1) with ρ = /2). Similarly to the Yang-Baxter deformation, the λ-deformation can be generalised to symmetric-space σ-models [13] and superstrings on semi-symmetric spaces [14].
The existence of a Lax pair ensures that the model admits an infinite number of conserved charges, extracted from the monodromy of the Lax matrix L(z) = 1 2 L + (z) − L − (z) . In order to show the integrability of this model, one has to prove that these conserved charges are in involution. For integrable σ-models, this is done by showing that the Poisson bracket of the Lax matrix takes the form of a non-ultralocal Maillet bracket [15,16]. This was proved for the Principal Chiral Model in [17], for the Yang-Baxter model, with and without Wess-Zumino term, in [6] and [9] and for the λ-model in [13]. For all these cases, the Maillet bracket takes a particular form, which is encoded in a rational function ϕ(z) of the spectral parameter, called the twist function [17][18][19][20] (see also [21]). These results shed light on the common algebraic structure underlying the integrability of this family of models and led to their reinterpretation as part of a larger class of integrable field theories, called (realisations of) affine Gaudin models [22]. In this formalism, the twist function and the Lax matrix of the model arise naturally from representations of untwisted affine Kac-Moody algebras.
Recently, the formalism of affine Gaudin models has been applied to generate an infinite family of new integrable classical σ-models [23][24][25]. More precisely, these models are obtained by coupling in a non-trivial way an arbitrary number of Principal Chiral Models with Wess-Zumino terms on the same Lie group G 0 . The fact that these σ-models are constructed as realisations of affine Gaudin models ensures that they are integrable (more precisely, they possess a Lax pair, whose spatial component satisfies a Maillet bracket with twist function). Let us briefly describe the coupled model with N copies. In order to write its action in a compact form, we gather the N G 0 -valued fields of the model into a unique field f = g (1) , · · · , g (N ) ∈ G N 0 and consider the block diagonal bilinear form κ N on g N 0 , which restricts on each factor g 0 to κ. The action of the coupled model then takes the form . (1.2) In this equation, denotes a linear operator on g N 0 , which can be seen as a N × N matrix of operators on g 0 with entries rs = ρ rs Id proportional to the identity. The action (1.2) is then characterised by the N (N + 1) coefficients ρ rs and r . For generic values of these coefficients, the model is not integrable. The particular model of [23,24], which is integrable since it is constructed as a realisation of affine Gaudin model, corresponds to a specific choice of these coefficients. More precisely, they are expressed in terms of 3N − 1 free parameters in a way which, for brevity, we will not describe in this introduction. The Lax connection of the model takes the form where (f −1 ∂ ± f ) r is the r-th component g (r) −1 ∂ ± g (r) of f −1 ∂ ± f ∈ g N 0 and the α r (z)'s are rational functions of the spectral parameter, whose expressions in terms of the 3N − 1 defining parameters of the model are also known explicitly [23,24].
It is natural at this point to search for integrable deformations of this coupled σ-model. It was explained in [24] that such deformations exist and that they can also be defined as realisations of affine Gaudin models. For instance, one can apply a Yang-Baxter deformation to any of the N copies of the model. Moreover, if one of the copy has no Wess-Zumino term, it is possible to consider a corresponding λ-deformation, which would then be more precisely a deformation of the model where this copy of the Principal Chiral Model has been replaced by its non-abelian T-dual. In general, one can then consider any combinations of these deformations on the different copies, leading to a whole panorama of different models.
Although these integrable deformed coupled σ-models are known to exist, they have not been constructed explicitly so far and have yet to be fully understood. In particular, since they are defined as realisations of affine Gaudin models, they are inherently formulated in the Hamiltonian framework. It is then an important aspect in the understanding of these models and of their properties to formulate them in the Lagrangian framework and to find an explicit expression of their action. In particular, this would give us access to the geometry underlying these theories, i.e. the deformed metric and B-field of the target space G N 0 which define these σ-models. In addition to clarifying the structure of the models at the classical level, describing their Lagrangian formulation can also benefit the understanding of their quantum properties, as for example the one-loop renormalisation of σ-models is controlled by the curvature of their underlying geometry. It is also an important problem to express the Lax pair of the model in terms of the Lagrangian fields, in order to understand how the integrable structure of the model manifests itself in the Lagrangian formulation. The explicit construction of the action and Lagrangian Lax pair of these integrable deformed coupled σ-models is the main subject of this article.
Several examples of integrable σ-models coupling together λ-models were proposed by Georgiou and Sfetsos in [26][27][28][29], using a different approach than the one considered in this article. Moreover, it was shown very recently that these models satisfy a Maillet bracket and possess a twist function [30]. As an application of the general construction developed in this article, we will show that these models can be obtained as limits of the ones obtained using affine Gaudin models.
Before sketching the methods used in this article to construct and study integrable deformations of coupled σ-models, let us illustrate briefly some of its main results. Let us first consider the model with N copies of the Principal Chiral Model with Wess-Zumino term, each subject to a Yang-Baxter deformation. It is defined by 4N − 1 parameters, which can be thought of as the 3N − 1 parameters of the undeformed model and N deformation parameters, and by the choice of N R-matrices R r on g 0 1 . The action of the model then takes the form In this expression, ± and U ± are operators on g N 0 and t denotes the transpose of operators with respect to the bilinear form κ N . Similarly to in the undeformed model, the entries of the operators ± are of the form ( ± ) rs = ρ ± rs Id, with the coefficients ρ ± rs expressed explicitly in terms of the 4N − 1 1 The R-matrix Rr is assumed to satisfy the additional property R 3 r = c 2 r Rr, except if the r-th copy does not possess a Wess-Zumino term, i.e. if r = 0.
defining parameters of the model. The entries U ± rs in the r-th row of the operators U ± are expressed as polynomials of degree two (or degree one if r = 0) in Ad −1 g (r) •R r •Ad g (r) , with coefficients also explicitly known as functions of the defining parameters. In the undeformed limit (where the N deformation parameters are taken to 0), the coefficients ρ + rs and ρ − rs converge respectively to ρ sr − δ rs r /2 and ρ rs + δ rs r /2 and the operators U ± simply become the identity operator on g N 0 . The action (1.3) then reduces to the action (1.1) of the undeformed model. The operators U ± entering the action of the model also control the Lax pair of the model. Indeed, the latter is given by with α r (z) explicit deformations of the rational functions α r (z) considered in the undeformed model. Let us now consider the model coupling together N copies of the λ-model. This is a deformation of the model coupling N copies of the non-abelian T-dual of the Principal Chiral Model, which is equivalent to the model (1.2) with no Wess-Zumino terms. This undeformed model then possesses 2N − 1 free parameters and its N -fold λ-deformation is described by 3N − 1 parameters. The action of this model takes the form where K and M are operators on g N 0 , which can be seen as N × N matrices with scalar entries K rs = δ rs r Id and M rs = µ rs Id. The model is then characterised by the coefficients r and µ rs , which are expressed explicitly in terms of the 3N −1 defining parameters of the models. Actions of this form were already considered in the article [29]. In particular, it was argued in this reference that the truncation of this model where all the coefficients µ rs vanish except for the coefficients µ 11 , · · · , µ (N −1)1 and µ N 2 , · · · , µ N N defines an integrable model with 3N − 2 parameters. This truncation can be seen as a particular limit of the model constructed above, with one deformation parameter less. Although the model considered here extends this truncation by introducing only one additional parameter, this extension has a non-trivial effect on the structure of the model, as all the coefficients µ rs become generically non-zero in this model.
In the main text of this article, we also construct explicitly the integrable model coupling together N 1 Yang-Baxter models and N 2 λ-models, whose action takes a form which mixes the structures of the above actions (1.3) and (1.4). For brevity, we will not describe this action in the introduction. All these deformed models involve the inverse of operators on g N 0 . These operators can be seen as N × N matrices whose entries are operators on g 0 . In particular, the non-commutativity of these entries makes the explicit inversion of these operators a non-straightforward problem. In the case of models with two copies only, we show how to perform this inversion explicitly. More precisely, we find an expression of these inverse operators which involves inversions of operators on one copy of g 0 only. Using this result, we give more explicit expressions of the models coupling together two Yang-Baxter models or two λ-models.
Let us now briefly sketch the methods used in this article, which are based on the formalism of affine Gaudin models. To illustrate these methods, it is useful to come back to the deformed models with only one copy and describe their structure as realisations of affine Gaudin models. An important object characterising affine Gaudin models is their twist function, which is the rational function of the spectral parameter controlling the Poisson bracket of their Lax matrix. For the Yang-Baxter model and the λ-model, this twist function possesses two simple poles z ± in the complex plane. Each of these poles z ± corresponds to a so-called site of the underlying affine Gaudin model and is associated with a Kac-Moody current J ± . This current belongs to the Poisson algebra of Hamiltonian observables of the model, which for both the Yang-Baxter model and the λ-model is the algebra A G 0 of canonical fields on the cotangent bundle T * G 0 . The currents J ± satisfy the standard Poisson brackets of Kac-Moody currents and Poisson commute one with another.
It is a standard result in the literature that the Hamiltonian integrable structure of the Yang-Baxter and λ-models is characterised by two commuting Kac-Moody currents [6,9,13,31]. This is what motivated their reinterpretation as affine Gaudin models in [22]. An important remark to make here is that although the Yang-Baxter model and the λ-model both possess commuting Kac-Moody currents in the same Poisson algebra A G 0 , the expression of these currents in terms of the fundamental fields of this algebra is different. It is this expression which characterises the model one considers and in particular differentiates the Yang-Baxter model and the λ-model. In the terminology of [24], the datum of N commuting Kac-Moody currents in a Poisson algebra is called a Kac-Moody realisation with N sites. In particular, the Yang-Baxter and λ-models define two different Kac-Moody realisations with two sites, in the same algebra A G 0 .
The integrable coupled deformed models considered in this article are constructed as realisations of affine Gaudin models with 2N sites. Their twist function is thus a rational function of the spectral parameter z with 2N simple poles, that we gather in pairs z ± r , r ∈ {1, · · · , N }. We attach to these pairs N independent copies of either the Yang-Baxter realisation or the λ-realisation. The algebra of observables of these models is then the N -fold tensor product A ⊗N G 0 , which can be seen as the algebra of canonical fields on T * G N 0 . The models are then defined in the Hamiltonian framework: in particular, their Hamiltonian is constructed as the spatial integral of a particular quadratic combination of the 2N Kac-Moody currents attached to the 2N sites, following the general formalism of affine Gaudin models [22,24]. This definition of the Hamiltonian ensures that these models are integrable: their equation of motion can be recast as a zero curvature equation on a Lax pair and the corresponding Lax matrix satisfies a Maillet bracket, controlled by the choice of twist function made above. As these models are defined in the Hamiltonian framework, one then has to perform an inverse Legendre transform to obtain their Lagrangian formulation, and in particular their action and Lagrangian Lax pair. In this article, we do this using interpolation methods, which generalise to the deformed case the techniques used in [24] to treat the undeformed coupled model. The algebra of Hamiltonian observables of the model being the algebra of canonical fields on T * G N 0 , the models are formulated in the Lagrangian framework in terms of N G 0 -valued fields g (r) , which are the fields introduced earlier.
The most important building blocks for the construction of these models are the Yang-Baxter and λ-realisations, which are Kac-Moody realisations in A G 0 . In this article, we treat these two realisations in a uniform way, by introducing a general ansatz for the form of the corresponding Kac-Moody currents in A G 0 , which includes these two examples. Using the fact that this ansatz should describe Kac-Moody currents, we identify certain key properties that it should satisfy in general. These properties then allow us to obtain a general expression for the action and the Lax pair of the models based on the combinations of any number of Kac-Moody realisations obeying an ansatz of this form. We then apply these results to the case of a model constructed from Yang-Baxter realisations and/or λ-realisations. In this case, the particular form of these realisations allows for further simplifications of the action, which for instance lead to the action introduced above for the cases with N Yang-Baxter realisations or N λ-realisations.
As a side result, we comment in this article on the relation of the models constructed here with the 4d semi-holomorphic Chern-Simons theory. This theory was introduced in [32] and was related to integrable systems and in particular integrable lattice models in [33][34][35][36]. More recently, it was shown in [37] how to generate integrable two-dimensional field theories from this four-dimensional theory (see also [38][39][40] for further developments). The reference [37] treated two different classes of models, corresponding to so-called order and disorder defects. In particular, the Principal Chiral Model with Wess-Zumino term (1.1) and its coupled version (1.2) were obtained in this formalism as models with disorder defects. The canonical analysis of the general models with disorder defects was performed in a subsequent article [39], which showed in particular that all these integrable field theories are realisations of affine Gaudin models. Finally, it was shown in [40] how to obtain the Yang-Baxter model and the λ-model in this framework. It is thus natural to search for a construction of the deformed coupled σ-models considered here from the 4d semi-holomorphic Chern-Simons theory. In this article, we present this construction explicitly and relate it to the affine Gaudin model approach.
The plan of this article is the following. In Section 2, we explain the construction of the models in the Hamiltonian framework. More precisely, we first describe in details in Subsection 2.1 the Kac-Moody realisations in A G 0 that serve as building blocks for this construction. We then proceed to construct the models as realisations of affine Gaudin models in Subsection 2.2. We go on to perform the inverse Legendre transform of these models in Section 3, constructing in particular their action and their Lagrangian Lax pair. The results of Sections 2 and 3 are obtained using the general ansatz for the Kac-Moody realisations mentioned above in this introduction. We then study the models obtained from combinations of Yang-Baxter realisations and λ-realisations in Section 4: in particular, we find a simple expression of the action of these field theories and show that the σ-models constructed in [26][27][28][29] can be obtained as particular limits of the ones constructed in this section. Finally, in Section 5, we explain the relation of this work with the 4d semi-holomorphic Chern-Simons theory. Some technical results are gathered in Appendices A and B.

Hamiltonian formulation
In this section, we define the integrable field theories that we will consider in this article. These theories are constructed as realisations of affine Gaudin models (AGM), following the general terminology of [22,24], and as such are then naturally defined in the Hamiltonian formalism. As explained in [22,24], the basic building blocks for the construction of realisations of AGM are the so-called Takiff realisations. In this article, we will be interested in a particular class of such realisations, which are given concretely by a pair of Kac-Moody currents in a certain Poisson algebra. As these particular Kac-Moody realisations are the basic building blocks of the models we will consider, we will start by describing them in details in Subsection 2.1, before proceeding to the construction of the models themselves in Subsection 2.2. For conciseness, we will not reintroduce here the general formalism of AGM and Takiff realisations and refer to [22,24] for the details.

Kac-Moody realisations in
All the Kac-Moody realisations that we shall consider in this article are defined in the same Poisson algebra A G 0 , which is the algebra of canonical fields on the cotangent bundle of a Lie group G 0 . Let us then begin by describing this algebra.
Conventions and notation. Let us consider a finite-dimensional semi-simple real Lie algebra g 0 . Let us also introduce the opposite of its Killing form κ, which is a non-degenerate bilinear form on g 0 . We will denote a basis for g 0 by (I a ) a∈{1,··· ,n} and its dual basis with respect to κ by (I a ) a∈{1,··· ,n} . It is then possible to define the split quadratic Casimir of g 0 as the following element: in g 0 ⊗ g 0 , which is independent of the choice of basis (here and in the following, we use the standard tensorial notations i).
The Lie algebra g 0 can be seen as the real form of a complex Lie algebra g, or, in other words, as the subalgebra of fixed points of an antilinear involutive automorphism τ of g. A basis for g over C is then given by (I a ) a∈{1,··· ,n} . We note that the split quadratic Casimir (2.1) of the algebra is real, in the sense that it satisfies τ 1 C 12 = τ 2 C 12 = C 12 .
To conclude, let us also mention the fact that by choosing g 0 to be the compact form of g, the bilinear form κ becomes a positive scalar product on g 0 .
Canonical fields on T * G 0 . Let G 0 be a connected real Lie group with Lie algebra g 0 . We will now consider fields taking values in the cotangent bundle T * G 0 and depending on a space coordinate x in a one dimensional space D, which for us will be either the real line R or the circle S 1 . Let us also consider the algebra of functionals generated by these fields 2 , which we will denote by A G 0 . It can be conveniently described in the following way.
Firstly, we note that acting by a translation on the base G 0 of T * G 0 , it is always possible to send the cotangent space at a point p ∈ G 0 to the one at the identity Id ∈ G 0 , which is just the dual g * 0 of the Lie algebra g 0 . As we supposed g 0 to be semi-simple, we then have a canonical isomorphism between g * 0 and g 0 through the bilinear form κ. This further implies that also T * G 0 and G 0 × g 0 are isomorphic to each other. Hence, it is possible to describe a field on T * G 0 by a pair of fields g(x) in G 0 and X(x) in g 0 . Now, as T * G 0 is a cotangent bundle, it possesses a canonical symplectic structure. This means that A G 0 comes naturally equipped with a Poisson bracket, which, in terms of the fields g and X, can be written as where C 12 is the split quadratic Casimir of g 0 and δ xy = δ(x − y) is the Dirac delta-distribution.
The current j(x) and the momentum. Let us now define the following g 0 -valued current: which, from (2.2), satisfies the Poisson brackets Let us also consider the quantity From (2.2) and (2.3), one can check that its Hamiltonian flow generates the spatial derivatives on both g(x) and X(x): and Hence, this is nothing but the momentum of the algebra of observables A G 0 .
The current W (x) and the Wess-Zumino term. As shown for example in [24], it is also possible to define another g 0 -valued current W (x) in A G 0 with Poisson brackets and However, in this article, we shall not need the precise definition of W (x) and thus we refer to [24] for details. A further property of this current is that it satisfies the following orthogonality relation: As a final remark, we note that through this current it is possible to define the Wess-Zumino term of g [41][42][43]. Indeed, briefly considering the field g to be dependent on a time coordinate t ∈ R (in the Hamiltonian formulation, this time dependence is implicitly defined by the choice of a Hamiltonian), the Wess-Zumino term of g is given by (see for instance [24]) (2.8)

Kac-Moody currents
Commuting Kac-Moody currents. We are now in a position to introduce the Kac-Moody realisations that will serve as basic building blocks for the construction of the integrable models of Subsection 2.2. Such realisations are characterised by two commuting Kac-Moody currents in A G 0 , i.e. two g-valued fields J ± (x) satisfying the Poisson brackets where ± are constant numbers called the levels. Currents of such kind have already been found to play an important role in the study of integrable deformations of σ-models [6,9,13,31], leading to examples of Kac-Moody realisations such as the Yang-Baxter realisation (with or without Wess-Zumino term) and the λ-realisation [24]. These examples will be described more in detail in Subsection 2.1.3. For the time being, we focus on aspects which are common to all the realisations we shall describe, in order to keep the treatment as general and uniform as possible.
In particular, in all the examples we shall consider, the Kac-Moody currents J ± (x) are expressed as linear combinations of the g 0 -valued currents X(x), j(x) and W (x) introduced in Subsection 2.1.1. Moreover, the currents X(x) and W (x) always appear through the unique combination for some real constant which depends on the particular realisation. As one can see from (2.8), the current W is related to the Wess-Zumino term of the corresponding field g. Because of this relation, and as we will see more precisely in Subsection 3.2, the presence of the current W in the realisation, i.e. the non-vanishing of , will lead to the presence of a corresponding Wess-Zumino term in the action of the model.
From now on, we will suppose that the Kac-Moody currents J ± (x) take the form where B ± , C ± : g → g are linear operators on the Lie algebra g. We will allow these operators to be dynamical (and thus have non-trivial Poisson brackets with other quantities in A G 0 ), but will suppose them to depend only on the field g (that is, not on X or derivatives of g). As we shall see in Subsection 2.1.3, both the Yang-Baxter realisation and the λ-realisation can be retrieved in this formalism by making some specific choices for the operators B ± and C ± . Let us note that, in general, these operators cannot be arbitrary. Indeed, they should be chosen such that the currents (2.10) satisfy the brackets (2.9). We will not try to write here the most general conditions on B ± and C ± for these brackets to hold. However, as explained in details in Appendix A, one can already obtain some useful constraints on these operators by focusing on the non-ultralocal terms in the brackets (2.9), i.e. terms proportional to the derivative of the Dirac distribution. More precisely, one finds that B ± and C ± should satisfy the following identities: where we have introduced the transpose t O with respect to the form κ for an operator O on the Lie algebra g.
Reality conditions. In order for the models that we will construct from these realisations to be real, one has to impose some reality conditions on both the currents J ± and the levels ± . There are two possible types of conditions that we shall consider. In the first case, we suppose that the currents are invariant under the antilinear involutive automorphism τ (i.e. they are g 0 -valued) and the corresponding levels are real: τ (J ± (x)) = J ± (x) and ± = ± . (2.12) In the second case, one requires the currents to be conjugate with respect to τ and the levels to be complex conjugate to each other: Momentum and suitability of the realisation. To conclude this section, we will now prove that all the realisations that we are considering here are suitable (in the language of [24]). In particular, this will later allow a simple characterisation of Lorentz invariance for the integrable models that we will build from them.
To start with, it is simple to check that from the relations (2.11) obeyed by the operators B ± and C ± , one can derive the following additional identities: These, together with the definition of the currents (2.10) above and the identity (2.7), allow one to prove that the momentum (2.4) of the Poisson algebra A G 0 can be re-expressed as From [22,24], one recognises on the right-hand side the Segal-Sugawara integrals of the Kac-Moody realisation. This implies that the realisations described above are indeed suitable.

Examples of realisations
We will now review some relevant examples of Kac-Moody realisations.
Inhomogeneous Yang-Baxter realisation without Wess-Zumino term. Let us start by considering a solution R : g 0 → g 0 of the modified classical Yang-Baxter equation (mCYBE): with c = 1 (so-called split case) or c = i (non-split case), which we suppose to be skew-symmetric with respect to the non-degenerate form κ: The Kac-Moody currents for the inhomogeneous Yang-Baxter realisation without Wess-Zumino term then read [6,24,31] where γ is a real constant and The proof that these are Kac-Moody currents can be found in [6], where the levels are found to be . (2.17) Note in particular that the levels ± are opposite to one another. Moreover, the reality conditions discussed in Subsection 2.1.2 are satisfied. In particular, in the split case (c = 1) the currents J ± are g 0 -valued and the levels ± are real, hence (2.12) is satisfied. In the non-split case (c = i) instead, it is a simple check that the currents and the levels satisfy (2.13).
In the general language of Subsection 2.1.2, we see that the current W does not appear in the expression (2.16), which means that for this realisation we take the coefficient in (2.10) to be zero. According to what has been discussed in the previous subsection, this justifies the fact that the models constructed from this realisation will not contain the Wess-Zumino term of g. Finally, by comparing with (2.10) we read for the operators B ± and C ± : One easily checks that these operators satisfy the identities (2.11), as expected.
Inhomogeneous Yang-Baxter realisation with Wess-Zumino term. The inhomogeneous Yang-Baxter realisation defined in the previous paragraph has no Wess-Zumino term, i.e. does not contain the current W (x) (or equivalently has = 0). Following [9], one can generalise this construction to include the current W (x) and thus a non-zero coefficient , at least when the R-matrix underlying the realisation satisfies the additional condition R 3 = c 2 R, with c as in the right-hand side of the mCYBE (2.15) (note in particular that the standard Drinfeld-Jimbo R-matrix satisfies this condition). The levels of this generalised realisation are given by with γ a real constant. Comparing to the levels (2.17) of the realisation without Wess-Zumino term, one sees that turning on the coefficient corresponds to relaxing the fact that the levels ± are opposite one to another.
The Kac-Moody currents of the inhomogeneous Yang-Baxter realisation with Wess-Zumino term can be computed from the results of [9, Section 3], up to a few differences in the conventions 3 . In the present notations, they read where we recall that Y = X − W and where we have defined the quantities Similarly to the case without Wess-Zumino term, it is simple to check that the reality conditions are satisfied for both the choices c = 1 and c = i. From the form of the currents, one can make the following identifications comparing to Equation (2.10): Let us note that, as expected, the identities (2.11) are again satisfied by these operators B ± and C ± (using the fact that we restrict here to R-matrices satisfying R 3 = c 2 R).
λ-realisation. For the λ-realisation, the Kac-Moody currents are given by [13,24,31]: Note that, similarly to the inhomogeneous Yang-Baxter realisation without Wess-Zumino term, these levels ± are opposite one to another. In this case, the reality condition (2.12) is satisfied, since the currents J ± are g 0 -valued and the levels ± are real.
To conclude, comparing to Equation (2.10), one sees that for the λ-realisation the operators B ± and C ± have the following form: 20) and again one can check that the identities (2.11) are satisfied.

Definition as realisations of affine Gaudin models
Sites, levels and twist function. In this subsection, we proceed to constructing the integrable field theories that we will consider in this article as realisations of AGM, following the general formalism and terminology of [22,24]. As AGM, the models that we will consider possess 2N sites of multiplicity one, which we gather in pairs (r, +) and (r, −) with r ∈ {1, · · · , N }. The position of the site (r, ±) in the complex plane C will be denoted by z ± r . Since each site (r, ±) is of multiplicity one, it is associated with one level, which is a non-zero constant number and which we will denote by ± r . Following [22,24], let us also fix a non-zero real number ∞ . Altogether, this data specifies the so-called twist function of the AGM, which in the present case reads where z ∈ C is an auxiliary complex parameter, called the spectral parameter.
Kac-Moody currents and the algebra A. To each site (r, ±) we attach a g-valued field J (r) ± (x) in the algebra of observables of the model A, which we now describe. As explained in [22,24], the Poisson brackets of these fields are specified by the choice of levels ± r made above. More precisely, we have the following: Thus, the models that we consider are constructed from N independent pairs of commuting Kac-Moody currents (J (r) , · · · , N }. We have described in detail in Subsection 2.1 how such a pair can be realised in the algebra A G 0 . A natural way to realise the 2N currents J (r) ± is then to consider N independent realisations in A G 0 of the type described in Subsection 2.1. Concretely, this means that we choose the algebra of observables of the models to be the tensor product A = A ⊗N G 0 , with the currents J (r) ± belonging to the r th -tensor factor 1 ⊗ · · · ⊗ A G 0 ⊗ · · · ⊗ 1 in A. This r th -tensor factor is generated by a pair of canonical fields g (r) (x) and X (r) (x), valued respectively in the group G 0 and the Lie algebra g 0 , which are the equivalent of the fields g(x) and X(x) introduced in Subsection 2.1.1 in order to describe one copy of A G 0 . Similarly, one can define from these canonical fields the equivalent in the r th -tensor factor of the currents j(x) and W (x), which we shall denote by j (r) (x) and W (r) (x). Following the discussion above, we then also define the currents J (r) ± as the analogues in the r th -tensor factor of the Kac-Moody currents J ± described in Subsection 2.1.2. Therefore, they take the form and r is a real constant number depending on the choice of representation in the r th -tensor factor. The B (r) ± 's are linear operators on the Lie algebra g, which are the equivalent in the r-th tensor factor of the operators B ± and C ± introduced in Subsection 2.1.2. In particular, they depend only on g (r) and satisfy analogous identities to the ones of equation (2.11).
Gaudin Lax matrix. We are now in a position to define the remaining building block that, in the next section, will allow us to write down an Hamiltonian for the model. This is the so-called Gaudin Lax matrix of the model, which we define as the following g-valued field [22,24]: (2.25) Reality conditions. As we discussed in Subsection 2.1.2, in order for the models which we construct in this article to be real, we have to impose some reality conditions. For each pair of sites (r, ±), there are two cases. In the first one, we suppose the positions of the two sites z ± r to be real and that the condition (2.12) on the currents J (r) ± and the levels ± r holds. In the second case, we assume instead that the the positions of the sites are complex conjugate to each other and that the currents and levels satisfy the condition (2.13).
These conditions can be summarised in terms of the twist function and the Gaudin Lax matrix as the following equivariance relations: and ϕ(z) = ϕ(z).

Hamiltonian and momentum
Hamiltonian. In order to construct the Hamiltonian of the model, we start by rewriting the twist function in terms of its zeroes ζ i (i ∈ {1, · · · , 2N }), which, for future convenience, we will suppose to be real and distinct. As we assumed ∞ to be non-zero, we can thus rewrite the twist function as .

(2.26)
Let us consider the spectral parameter dependent local charge and define, for i = 1, · · · , 2N , or, more explicitly, These are local charges quadratic in the currents J (r) ± which, as proven in [24], satisfy {Q i , Q j } = 0 for all i and j. We define the Hamiltonian of the model to be the linear combination for some real numbers i . This then generates the time evolution of the model through the Hamiltonian flow Note that, as a consequence of the reality conditions we introduced, H is real [24].
Momentum and relativistic invariance. Recall that, in Subsection 2.1.2, we proved that the Kac-Moody realisations in A G 0 that we are considering are suitable. According to [24], this gives some additional information on the space-time properties of the model. Firstly, the momentum of the algebra A (i.e. the generator of spatial translations with respect to x) is given by the following expression: Secondly, requiring relativistic invariance of the model restricts the choice of the coefficients i in the definition of H to for every i ∈ {1, · · · , 2N }. We then see that there is a natural division of the indices i ∈ {1, · · · , 2N } labelling the zeroes ζ i into the sets In the rest of this article, we will suppose that there are as many i 's equal to +1 as i 's equal to −1 (i.e. that the sets I ± are both of size

Integrability
Lax pair and zero curvature equation. We define the Lax matrix of the model to be the following g-valued field [22,24]: By construction, it has poles at the zeroes ζ i of the twist function. More precisely, it can be rewritten as [24] L(z, From here, one can check that the time evolution of L(z, x) takes the form of a zero curvature equation, Therefore, the model admits a Lax pair representation with Lax pair (L, M).
Maillet bracket and integrability. The fact that the equations of motion of the model take the form of a zero curvature equation allows one to extract an infinite number of charges from the monodromy of the Lax matrix L(z, x). The integrability of the model then follows from the fact that these are in involution, which is a consequence of the Lax matrix satisfying the following Poisson bracket: where the R-matrix is defined to be The bracket (2.34) is an example of a Maillet non-ultralocal bracket [15,16]. One can check that it satisfies the Jacobi identity due to the fact that the R-matrix is a solution of the classical Yang-Baxter equation Lax pair in light-cone coordinates. As we will need this in Section 3, let us briefly discuss the reparametrisation of the Lax pair in light-cone components. Let us firstly introduce the light-cone coordinates x ± = (t ± x)/2 and the corresponding derivatives ∂ ± = ∂ t ± ∂ x . The zero curvature equation (2.32) can then be rewritten as where we have introduced the light-cone Lax pair Finally, from Equations (2.31) and (2.33), one finds the following expression for L ± (z, x): in terms of the split of the zeroes into the two sets I ± introduced in Subsection 2.2.2.

Exploring the "space of models"
Gaudin parameters. Let us describe the "space of models" that we are considering in this article by summarising what are the defining parameters of the integrable field theories that we have constructed so far. As affine Gaudin models, these theories are characterised by the following quantities, that we shall refer to as Gaudin parameters: • the positions z ± r ; • the levels ± r ; • the constant term ∞ in the twist function ; • the Kac-Moody realisations with levels ± r attached to each pair of sites (r, ±). As explained in [24,Subsection 1.4.2], there exists a redundancy between the Gaudin parameters of the model, corresponding to the freedom of translating and dilating the spectral parameter. Indeed, the model with parameters z ± r , ± r and ∞ as above is invariant under the transformation where a and b are real numbers with a = 0 and where we keep the levels ± r and the Kac-Moody realisations fixed. Note that one can fix the dilation redundancy (corresponding to the parameter a in the transformation above) by setting the constant term ∞ to a specific value. Similarly, one can fix the translation redundancy (corresponding to the parameter b) by setting one of the positions z ± r to a specific point.
Note that the Gaudin parameters introduced above are in general not all real but should satisfy the reality conditions described in Subsections 2.1.2 and 2.2.1. Let us then discuss what are the real parameters of the models. Note first that the constant term ∞ is always assumed to be real. Moreover, recall that for each pair of sites (r, ±), there are two possible reality conditions: either the positions z ± r and the levels ± r are real or they form pairs of complex conjugate numbers. We will encode the choice of reality condition for the sites (r, ±) by introducing a number c r , which is defined to be 1 in the first case and i in the second one. In particular, z ± r and ± r can then be written using the following parametrisation: where the parameters z r , η r , r and [1] r are real. As we shall see, this particular choice of parametrisation will also be convenient for the interpretation of the models as deformations in the next subsection. Note that it is equivalent to defining ± are given by Equation (2.19) (replacing g by g (r) and c by the number c r ∈ {1, i} defined in the previous paragraph, which encodes the choice of reality conditions for the sites (r, ±)).
Similarly, if one chooses the λ-realisation, the operators B (r) ± and C (r) ± are given by Equation (2.20) (with g replaced by g (r) ), while r is given by − + r /2. Note however that one can choose the λrealisation only if the levels ± r are real (i.e. c r = 1 in the notations of the previous paragraph) and are such that This is in contrast with the case of the inhomogeneous Yang-Baxter realisation with Wess-Zumino term considered above, where the levels ± r are not subject to any constraints (other than the reality conditions).
Note that the choice of a Yang-Baxter realisation at the sites (r, ±) comes with the additional freedom of choosing a skew-symmetric R-matrix R r , solution of the mCYBE (2.15). As explained in Subsection 2.1.3, this operator should in general satisfy the additional property R 3 r = c 2 r R r . However, if the levels ± r satisfy the constraint (2.39), i.e. if one considers a Yang-Baxter realisation without Wess-Zumino term, one does not need to require this additional condition on R r .
The space of models. The discussion above concerns the choice of realisation for one pair of sites (r, ±). One can then construct different models by considering different combinations of realisations for the N pairs (1, ±), · · · , (N, ±) describing the models. In particular, one can consider a model with N 1 copies of the Yang-Baxter realisation and N 2 copies of the λ-realisation, where N 1 + N 2 = N . Let us discuss what are the free parameters of this theory. As explained in the previous paragraphs, the model is described by the 4N + 1 Gaudin parameters z ± r , ± r and ∞ , or equivalently by the 4N + 1 real parameters z r , η r , r , [1] r and ∞ . Taking into account the translation and dilation redundancy (2.36) and the fact that the levels corresponding to the λ-realisations should satisfy the constraints (2.39), we arrive at the conclusion that this model is described by 3N + N 1 − 1 free parameters. Note that in addition to these parameters, which specify its structure as an AGM, the model is also determined by the choice of N 1 R-matrices for the Yang-Baxter realisations (which do not need to be identical).
As was explained in [22], see also [24], the models with only one realisation, i.e. with N = 1, correspond to well-known integrable σ-models, which served as basis for defining the Yang-Baxter and λ-realisations. Indeed, the inhomogeneous Yang-Baxter realisation (without or with Wess-Zumino term) is defined in such a way that the AGM with one copy of this realisation, corresponding in the above paragraph to N 1 = 1 and N 2 = 0, coincides with the so-called Yang-Baxter σ-model, without [3,4] or with [9] Wess-Zumino term. Similarly, the AGM with one copy of the λ-realisation, i.e. with N 1 = 0 and N 2 = 1, yields the so-called λ-model [11]. The model defined above with arbitrary numbers N 1 and N 2 is thus a generalisation of these models. According to the general coupling procedure described in [24,Subsection 2.3.3], it corresponds to coupling together N 1 copies of the Yang-Baxter model and N 2 copies of the λ-model in a non-trivial way which however ensures the integrability of this interacting model (as, by construction, it is a realisation of AGM).
Zeroes versus levels. Let us end this subsection with some remarks about a possible more convenient reparametrisation of the models that we are considering. Recall from Subsections 2.2.2 and 2.2.3 that in order to define the Hamiltonian and express the Lax pair of the models, one uses the zeroes ζ i , i ∈ {1, · · · , 2N }, of the twist function. These zeroes are related implicitly to the Gaudin parameters z ± r , ± r and ∞ through the equation ϕ(ζ i ) = 0, with the twist function ϕ(z) defined in terms of the Gaudin parameters as in (2.21). This equation is equivalent to a polynomial equation of degree 2N in ζ i . Thus, it is in general impossible to give an explicit expression of the zeroes ζ i in terms of the Gaudin parameters.
One way of bypassing this difficulty is to consider as defining parameters of the models the positions z ± r , the zeroes ζ i and the constant term ∞ . One then defines the twist function of the model by Equation (2.26) instead of Equation (2.21) and the levels ± r as the corresponding residues: .
The main advantage of this re-parametrisation is that all the relevant quantities that are used to describe the models, in particular the levels ± r and the Hamiltonian H, can be written as rational expressions of the parameters z ± r , ζ i and ∞ . Note however that this parametrisation has a disadvantage when one wants to consider λ-realisations and/or Yang-Baxter realisations without Wess-Zumino terms. Indeed, for these realisations, the levels should satisfy the additional constraint (2.39), which translates in a rather complicated algebraic condition on the parameters z ± r and ζ i , using the above expressions for the levels. Finally, let us note that the translation and dilation redundancy (2.36) among the Gaudin parameters can be re-expressed in terms of this new parametrisation as the invariance of the model under the transformation

Recovering undeformed models
In this subsection, following the results of [24], we discuss how the model defined above by taking N 1 Yang-Baxter realisations and N 2 λ-realisations can be interpreted as a deformation of a simpler model. This result generalises the well known facts that the Yang-Baxter model (with or without Wess-Zumino term) is a deformation of the Principal Chiral Model (PCM, with or without Wess-Zumino term) and the λ-model is a deformation of the non-abelian T-dual of the PCM. In the present language, these correspond respectively to the cases (N 1 = 1, N 2 = 0) and (N 1 = 0, N 2 = 1). The undeformed limit of the model with arbitrary N 1 and N 2 corresponds to a theory coupling together N 1 copies of the PCM (with Wess-Zumino terms) and N 2 copies of its non-abelian T-dual. This undeformed model is also defined as a realisation of AGM but possesses a sligthly different sites structure. Indeed, in the language of [22,24], instead of the 2N sites (r, ±) of multiplicity one, it possesses N sites (r) of multiplicity two. These sites correspond to double poles in the twist function and the Gaudin Lax matrix of the model and are associated with so-called Takiff realisations of multiplicity two, which generalise the notion of Kac-Moody realisations for sites of multiplicity greater than one. As we shall now explain, the site (r) of multiplicity two is obtained from the pair of sites (r, ±) in the deformed model by making their positions z + r and z − r collide, while controlling the behaviour of the corresponding levels ± r .
Colliding two simple poles into a double pole. Let us focus here on one pair of sites (r, ±). In order to isolate the parts of the twist function and the Gaudin Lax matrix of the model corresponding to this pair, let us rewrite them as where ϕ and Γ contain all the information related to the other sites. Using the parameters c r , z r , η r , [0] r and [1] r introduced in the previous Subsection (see Equation (2.38)), one can rewrite the twist function as As mentioned above, the undeformed limit corresponds to making the two positions z + r and z − r collide at the point z r and thus to taking η r → 0. In particular, this leads us to interpret η r as a deformation parameter. We aim here to recover, in the limit η r → 0, a model with a site of multiplicity two, i.e. with a double pole in its twist function. It is then clear from Equation (2.40) that this is the case if one supposes that the quantities r and [1] r stay finite when η r goes to 0. From now on, we will thus define the undeformed limit as taking η r → 0 while keeping [0] r and [1] r finite (let us note that the levels ± r of the sites (r, ±) then diverge, as one can see from Equation (2.37)). In this limit, the twist function becomes Following the terminology of [22,24], this corresponds to the twist function of an AGM with a site (r) of multiplicity two, with position z r and Takiff levels r and [1] r (and with the other sites, contained in ϕ(z), as in the deformed model).
A similar argument applies to the Gaudin Lax matrix of the model. Let us suppose that the Kac-Moody currents J (r) ± are such that the limits are finite. Then the Gaudin Lax matrix becomes in the undeformed limit: Thus, J [0] and J (r) [1] are the Takiff currents attached to the site (r) of the undeformed model 5 . Let us now discuss this undeformed limit for the Yang-Baxter realisation and the λ-realisation. 5 Starting from the Kac-Moody Poisson brackets (2.22) of the currents J (r) ± , one can indeed show that in the undeformed limit, the currents J (r) [0] and J (r) [1] satisfy the brackets of Takiff currents with levels r and [1] r .
From the Yang-Baxter to the PCM realisation. Let us suppose that the sites (r, ±) are associated with a Yang-Baxter realisation with Wess-Zumino term, as described in Subsection 2.1.3. Let us first note that for this realisation, the Wess-Zumino coefficient is given by r = − [0] r /2. In particular, the undeformed limit defined in the previous paragraph can then be seen as taking η r to 0 while keeping r and [1] r finite. Let us denote by R r the R-matrix associated with this Yang-Baxter limit and introduce R (r) = Ad −1 g (r) • R r • Ad g (r) and Π (r) = Id − R (r) 2 /c 2 r . The Kac-Moody currents of the realisation are then given by Let us now consider the undeformed limit, i.e. taking η r to 0 while keeping r and [1] r finite. One first observes that in this limit, the coefficient δ r tends to 0. Using this, one finds that the limits J and J (r) [1] defined in Equation (2.41) are indeed finite and simply read and J (r) Thus, the undeformed limit described in the previous paragraph is well defined. Moreover, one recognises in the above equation the Takiff currents of the PCM+WZ realisation (with levels From the λ-realisation to the non-abelian T-dual realisation. A similar mechanism to the one described above for the Yang-Baxter realisation provides the underformed limit of the λ-realisation, yielding the so-called non-abelian T-dual realisation, as defined in [24, Subsection 4.3.1]. This limit requires however a more subtle treatment. Indeed, if one were to consider the currents J (r) ± of the λrealisation in terms of the fields g (r) and X (r) and take the limits (2.41) "naively", one would encounter divergent expressions, making the undeformed limit procedure ill-defined. In order to obtain a well defined limit, one has to consider the fields g (r) and X (r) as depending on the deformation parameter η r and suppose that they obey a well-chosen asymptotic expansion when η r goes to 0. In particular, one of the consequences of this more subtle limit is that it changes the algebra of observables of the realisation: from the algebra A G 0 of canonical fields on T * G 0 (generated by g (r) and X (r) ), one goes in the limit to the algebra A g 0 of canonical fields on T * g 0 , which is the algebra supporting the nonabelian T-dual realisation. For brevity, we will not re-explain this procedure in the present article and refer to [24,Subsection 4.4

.3] for details.
Undeformed limits of the coupled models. Let us consider the model defined in the previous subsection by coupling together N 1 copies of the Yang-Baxter model and N 2 copies of the λ-model. For each pair of sites (r, ±), one can consider the corresponding undeformed limit η r → 0. One would then obtain a model where the r-th copy reduces to either an undeformed PCM with Wess-Zumino term or a non-abelian T-dual of the PCM (depending on whether we started with a Yang-Baxter realisation or a λ-realisation at the sites (r, ±)), still interacting non-trivially with the other N − 1 copies. One can then consider different combinations of these undeformed limits on any number of copies, yielding various limits of the model. All these limits can be seen as deformations of a completely undeformed model, obtained by taking the limit where all the deformation parameters η 1 , · · · , η N are sent to 0. This undeformed model is the coupling of N 1 copies of the PCM with Wess-Zumino terms and N 2 copies of the non-abelian T-dual of the PCM. In particular, if one considers N 2 = 0, one obtains the model coupling together N copies of the PCM with Wess-Zumino term: this is the integrable coupled σ-model first introduced in [23] and whose detailed construction was presented in [24,Subsection 3.3].
Although it is defined in a different way, let us note also that the undeformed model with N 2 = 0 copies of the non-abelian T-dual of the PCM is in fact canonically equivalent to the model with N = N 1 + N 2 copies of the PCM, where N 2 of these copies have no Wess-Zumino term. This is because the non-abelian T-dual realisation is related to the PCM realisation without Wess-Zumino term by a canonical transformation [44]. Thus, the general model with N 1 Yang-Baxter realisations and N 2 λ-realisations can be seen as a deformation of the model coupling N 1 PCM with Wess-Zumino term and N 2 PCM without Wess-Zumino term (which is a particular case of the model introduced in [23]) after having first T-dualised the N 2 copies without Wess-Zumino term.
Homogeneous Yang-Baxter limit. For completeness, let us end this subsection by mentioning briefly another possible limit of the models considered here, which corresponds to going from an inhomogeneous Yang-Baxter realisation to a homogeneous Yang-Baxter realisation 6 . Let us consider an inhomogeneous Yang-Baxter realisation without Wess-Zumino term and with R-matrix R, which satisfies the mCYBE (2.15). So far, we considered the coefficient c appearing in the mCYBE as being either 1 or i, depending on the type of reality conditions imposed on the realisation. However, one easily checks that the construction of the Yang-Baxter realisation as recalled in Subsection 2.1.3 holds without changes for any c = 0 (the realisation is then equivalent to the one with c = 1 or c = i by rescaling the matrix R). The homogeneous limit consists in taking the limit c → 0 of this realisation while also making the corresponding simple poles in the twist function collide (see for example [25]). Similarly to what happens for the undeformed limit described in this subsection, this yields a model with a site of multiplicity two, to which is attached the so-called homogeneous Yang-Baxter realisation, as defined in [24,Subsection 4.1.1]. This realisation corresponds to a deformation of the PCM realisation without Wess-Zumino term by a homogeneous R-matrix, i.e. a solution of the (non-modified) CYBE: which corresponds to the limit c → 0 of the mCYBE.
Summary. Although we introduced them as limits, the PCM, non-abelian T-dual and homogeneous Yang-Baxter realisations can be constructed independently, as was done for example in [24]. One can then consider AGM containing these realisations. In general, one can construct a model coupling together any combination of PCMs, non-abelian T-dual models, homogeneous and inhomogeneous Yang-Baxter models and λ-models. Up to taking appropriate limits, the present article then covers all these possibilities. In particular, one can obtain a model with N − 1 copies of the PCM and one homogeneous Yang-Baxter realisation: one then recovers the model studied in [24, Appendix D] as the simplest illustration of the various possible integrable deformations of coupled integrable σ-models.

Lagrangian formulation
In this section, our aim will be to describe the models introduced in Section 2.2 in the Lagrangian formulation. Recall that in the Hamiltonian formulation, the degrees of freedom of these models are the fields g (r) (x) and X (r) (x), describing canonical fields valued in N independent copies of the cotangent bundle T * G 0 . The fields g (r) (x) are the "coordinate fields" of the models, valued in the base G 0 of the bundle T * G 0 . The momentum fields conjugate to these coordinate fields, which take values in the fibers of the bundle, are then encoded in the fields X (r) (x) (see for instance [24, Subsection 3.1.1] for details). In order to pass to the Lagrangian formulation, one has to consider the coordinate fields g (r) (x, t) as depending explicitly on the time variable t ∈ R, defined by the Hamiltonian of the model, and express the momentum fields of the theory, encoded in X (r) , in terms of these Lagrangian fields g (r) (x, t) and their derivatives ∂ t g (r) (x, t) and ∂ x g (r) (x, t). Finally, one obtains the action of the model as a functional of g (r) (x, t) by performing an inverse Legendre transform on their Hamiltonian.
In the present case, we will obtain the Lagrangian expression of the fields X (r) in a rather indirect way. Indeed, as we shall see, these fields can be expressed naturally in terms of the Lax pair of the model. For this reason, we will start by determining the Lagrangian expression of the latter.

Lax pair in the Lagrangian formulation
Maurer-Cartan currents in terms of the Lax pair. Let us begin by considering the time evolution of the fields g (r) . In the Hamiltonian formulation, this is given by their Poisson bracket with the Hamiltonian. More explicitly, recalling the definition (2.28) of the latter, one expresses the temporal Maurer-Cartan current j (r) From the expression (2.27) of the charges Q i , we then have The Poisson bracket in the integrand is calculated by inserting the definition (2.25) of Γ(z, x), yielding: where we have also used the fact that J (s) ± is in the s th -tensor factor in A N G 0 and thus Poisson commutes with g (r) if r = s. In order to calculate the Poisson brackets on the right hand side we then use the definition (2.23) of the currents J (r) ± in terms of Y (r) and j (r) . Note that, firstly, the Poisson brackets of g (r) with j (r) vanish. Moreover, the brackets of g (r) with the operators B ± also give no contribution as we assumed that these operators depend only on g (r) . Thus, we have to take into account only the terms coming from the Poisson bracket of g (r) with Y (r) , so that g (r) where we have used the fact that for any operator O on g, one has O 2 C 12 = t O 1 C 12 . Putting everything together, we conclude that Using the expression of the temporal component of the Lax pair (2.33), this can be re-expressed in the following way: j (r) . Moreover, by repeating this argument replacing the Hamiltonian by the momentum P, expressed in terms of the charges Q i as in (2.29), and using the expression (2.31) for the spatial component of the Lax pair, one finds a similar relation for the currents j (r) : Therefore, using light-cone coordinates, we find that the Maurer-Cartan currents j (r) take the following rather simple form in terms of the Lax pair: Lagrangian Lax pair from interpolation. Our goal in this subsection is to find a Lagrangian expression of the Lax pair, i.e. an expression of L ± (z) in terms of the Maurer-Cartan currents j (r) ± . We note that Equation (3.1) relates these currents to the evaluations L ± (z + r ) and L ± (z − r ) of the Lax pair at the positions z + r and z − r . As we shall now explain, this relation is enough to reconstruct the expression of L ± (z) in terms of j (r) ± for all values of the spectral parameter z. Let us define J (r) for r = 1, · · · , N . From Equation (2.35), one sees that L ± (z) is a rational function of z with N simple poles, situated at the zeroes of the twist function ζ i , for i ∈ I ± (recall that we have supposed that the subsets I ± are both of size N ). It is a standard result that such a function is completely determined by its evaluation at N pairwise distinct points. In particular, L ± (z) can be expressed in terms of its evaluations at the positions z ± r , i.e. the currents J (r) ± introduced above. More precisely, one has the following interpolation formula (see also Lemma B.2 of [24]) .

(3.4)
We are now in a position to rewrite the Lax pair in terms of the currents j where we have defined In the following, we will see the operators U ± rs as the entries of some matrix operators U ± , so that U ± rs = (U ± ) rs . Note that U ± are then N × N matrices with non-commutative entries. To conclude, we rewrite the currents J where (U −1 ± ) rs denote the entries of the inverse of the matrix operators U ± . Reinserting now (3.7) in (3.3) then gives an expression of L ± (z) in terms of the currents j (r) ± : Note that this is a formal relation, as it involves the inverse of the matrix operators U ± . Performing explicitly this inversion is in general not straightforward because of the non-commutativity of the entries of U ± (for example, one cannot use the general expression for the inverse of a matrix in terms of its comatrix). We will explain in Subsection 3.2 how this is done explicitly in the case of two copies.
Different interpolations and factorisations of the twist function. We conclude this section by making an important remark about Equations (3.2) and (3.3). In these equations, we decided to express the component L + (z), resp. L − (z), of the Lax pair in terms of its evaluations at the positions z + r , resp. z − r . Let us stress here that this choice is arbitrary, as one could have chosen for example to interpolate L + (z) and L − (z) through their evaluations at the positions z − r and z + r respectively 7 . More generally, one could have considered the evaluations J (r) where the σ r 's take values in the set {+1, −1} for every r. The interpolation equation (3.3) would then become . (3.9) Following the method developed in the previous paragraph, one would then express the currents J  ± by a relation similar to Equation (3.7), with the operators U ± replaced by some different operators U ± . Re-inserting this expression in Equation (3.9) would then give L ± (z) in terms of j (r) ± , similarly to Equation (3.8). This expression can be shown to coincide with Equation (3.8) as one should expect, considering that they correspond to two ways of expressing the same object L ± (z). Similarly, all the methods and computations developed in the rest of this subsection can be applied starting from an arbitrary choice of interpolation, i.e. from an arbitrary choice of σ r 's: the end results (in particular the expression of the action of the model in terms of the Maurer-Cartan currents that will be obtained in the next subsection) can then be shown to be independent of this choice. For this reason, and to avoid unnecessary cumbersome notations, we will use in the rest of this article a particular choice of σ r 's, namely σ r = +1 for every r, corresponding to the choice made originally in the previous paragraph.
To conclude this paragraph, let us discuss a reinterpretation of the functions ϕ ±,r (z) appearing in the interpolation formula (3.9) and of the freedom encoded in the choice of σ r 's in terms of the twist function (2.26) of the model. Let us rewrite the latter in the following factorised form: . (3.10) The functions ϕ ±,r (z) can then be re-expressed as ϕ ±,r (z) = (z − z ±σr r ) ϕ ± (z). Moreover, we observe that the freedom in the choice of the σ r 's gets now reinterpreted as the existence of different ways of factorising the twist function. Indeed, redistributing the pairs of factors (z−z + r ) and (z−z − r ) associated to the paired sites (r, ±) into the definition (3.10) of ϕ ± (z) amounts to changing the values of the σ r 's 8 . In the rest of this article and in agreement with the notations of the previous paragraph, we will denote by ϕ ± (z) the functions ϕ ± (z) corresponding to the choice σ r = +1 for every r ∈ {1, · · · , N }.

Inverse Legendre transform and action of the models
Lagrangian expression of the momentum. We are now in a position to perform the first step towards writing down the inverse Legendre transform of the model, i.e. re-expressing the fields X (r) , which encode the momentum fields of the theory, in terms of Lagrangian fields. Let us first note from 7 Note in particular that the indices ± of L±(z) are conceptually totally unrelated to the labels ± of the positions z ± r . Indeed, the former are space-time indices corresponding to the light-cone directions in R × D while the latter are abstract labels distinguishing the two sites (r, +) and (r, −). 8 Note that contrarily to the poles z ± r , the zeroes ζi of the twist function cannot be redistributed differently between the functions ϕ+(z) and ϕ−(z), as they are naturally associated with one or the other depending on whether the index i belongs to the set I+ or I−. Equation (2.24) that the fields Y (r) and X (r) are related through the current W (r) . As explained in Subsection 2.1.1, this current W (r) is expressed in terms of the field g (r) and its spatial derivative (and not the momentum fields) and has thus a direct Lagrangian expression. Thus, finding the Lagrangian expression of X (r) is equivalent to finding the Lagrangian expression of Y (r) . As we shall now see, the latter is easier to find, using the Lagrangian expression of the Lax pair obtained in the previous paragraph. From the definition (2.30) of the Lax matrix L(z), one can prove that (see also [24,Equation (2.22) where to obtain the second equality we have used the definition (2.23) of the currents J (r) ± . Then, using the identities (2.14) satisfied by the operators B (r) ± and C (r) ± , we find the following expression for Y (r) : Using the light-cone components of the Lax pair, this can be rewritten as From the Lagrangian expression (3.3) of L ± (z), one then finds that where we have defined Similarly to the operators U ± rs in the previous subsection, we will see the operators V ± rs as the entries of some N × N matrix of operators V ± , so that V ± rs = (V ± ) rs .
Action in terms of j ± 's. The action of the models is obtained as the following inverse Legendre transform of the Hamiltonian (see for instance [24]): where both X (r) and H should be replaced by their expressions in terms of Lagrangian fields. Recalling the definitions (2.24) and (2.8), one can rewrite the action in terms of the fields Y (r) making the Wess-Zumino terms of g (r) appear: From here, reinserting the expression (3.12) of Y (r) in terms of the currents J (r) ± , we find: (1) , · · · , g (N ) ] = 1 2 (3.14) We note that the terms in the second line are not Lorentz invariant. However, one shows that these are cancelled by the term containing the Hamiltonian (for brevity, we give the proof of this result in Appendix B), so that we eventually get S[g (1) , · · · , g (N ) ] = 1 2 N r,s Action in terms of Maurer-Cartan currents. To conclude this subsection, we proceed to compute the expression of the action in terms of the j (r) ± 's only. This is done through the formal inversion relation (3.7). As a final result we obtain where we have defined O rs as the entries of the following matrix operator: Finally, using the identities (2.11), one proves that the second term in this definition is equal to the first one, so that we get: Model with two copies. In this paragraph, we give an explicit expression for the inversion of the operator matrices U ± and consequently for the coupling operator O in the case of a model with two copies only, i.e. with N = 2. In order to do so, one has to make a further assumption about the operators B ± . More precisely, we will suppose that they satisfy the following commutation relation Let us note that, crucially, this additional condition is satisfied by the Yang-Baxter realisation (with or without Wess-Zumino term) and the λ-realisation, as can be checked easily from Equations (2.18), (2.19) and (2.20) 9 .
As we have noted in Subsection 3.1, the fact that it is not straightforward to invert the operator matrices U ± is due to the non-commutativity of their entries. However, using the additional assumption (3.19) made on the operators B (r) ± , one shows that: Thus, even if the entries of U ± are not all commutative, this shows that the ones on a same line commute with one another. This fact will allow us to find an explicit expression of the inverse of U ± .
Let us introduce the operators If the entries U ± rs of U ± were commutative, the objects ∆ ± 1 and ∆ ± 2 would be equal and would correspond to the inverse of the determinant of the 2 × 2 matrix U ± . In the present case, these operators 9 It is not obvious whether this condition is an accidental property of these particular realisations or if it can be derived more generally as a consequence of the fact that J (r) ± are Kac-moody currents, as was for example the case for the identities (2.11) (see Appendix A). ∆ ± r are the inverse of non-commutative versions of the determinant. In terms of these, the inverse of the operator U ± is then given by Indeed, one checks explicitly that The property (3.20) then ensures that the off-diagonal terms vanish, while the definition (3.21) of the operators ∆ ± r is such that the diagonal terms are the identity operator, thus proving that the matrix (3.22) is the inverse of U ± 10 . The expression (3.22) is a non-commutative generalisation of the standard comatrix formula for the inverse of a 2 × 2 matrix, where in particular one takes into account the non-commutativity of the entries by considering different "inverse determinants" ∆ ± r in the different columns.
To give a more compact expression of the entries of U −1 ± , let us introduce the notationr, defined for every r ∈ {1, 2} byr ∈ {1, 2} \ r (i.e.1 = 2 and2 = 1). Then, one has Reinserting the above results into the expression (3.18) of the operator O, one can compute its entries O rs , which appear in the action (3.16) of the model, yielding

Parameters of the models
In Subsection 2.2.4, we have discussed what are the defining parameters of the models, from their construction as realisations of affine Gaudin models. Let us briefly give some additional comments on the subject in the light of the Lagrangian formulation of the models.
Functions ϕ ± (z). Recall the functions ϕ ± (z) and ϕ ±,r (z) = (z−z ± r )ϕ ± (z) introduced in Subsection 3.1. It is clear from the results of this section that these functions play an important role in describing the Lagrangian formulation of the models. For example, they are used to obtain the Lagrangian expression (3.3) of the Lax pair. Similarly, they enter the definitions (3.6) and (3.13) of the operators U ± and V ± , which are then used to express the operator O appearing in the action (3.16) of the model. Note that the definition of the operators U ± and V ± also involves the operators B (r) ± and C (r) ± , which characterise the choice of Kac-Moody realisations of the model. In particular, these realisations depend on the levels ± r . For completeness, let us thus note that the latter can also be expressed quite easily using the functions ϕ ± (z) and ϕ ±,r (z): indeed, we have Finally, let us note that these levels also determine the coefficients r of the Wess-Zumino terms in the action (3.16) of the model. Thus, the datum of the functions ϕ ± (z) is enough to describe completely the model in its Lagrangian formulation.
Parameters (z ± r , ν ± r , ∞ ). Recall that in Subsection 2.2.4, we discussed two possible sets of parameters for the model: the "Gaudin parameters" (z ± r , ± r , ∞ ) and the parameters (z ± r , ζ i , ∞ ), where the datum of the levels ± r has been replaced by the datum of the zeroes ζ i of the twist function. In particular, recall that the second parametrisation is in general more convenient as the zeroes play an important role in the description of the model and as they cannot be expressed explicitly in terms of the levels, whereas the levels can be expressed rationally in terms of the zeroes. Recall also that choosing the parametrisation using the zeroes is however less convenient to describe models with λrealisations and/or Yang-Baxter realisations without Wess-Zumino term. Indeed, these realisations require that the levels ± r satisfy the additional constraint (2.39), which translates into a complicated algebraic condition on the parameters (z ± r , ζ i ). These observations motivate the introduction of a third possible set of parameters (z ± r , ν ± r , ∞ ), which is in some sense intermediate between the two sets described above and which circumvents the various issues related to solving algebraic equations. In this parametrisation, the datum of the levels ± r or of the zeroes ζ i is replaced by the datum of the coefficients ν ± r = ϕ ±,r (z ± r ). Note that these coefficients characterise the partial fraction decomposition of the functions ϕ ± (z): In particular, the levels ± r can then be expressed in terms of these parameters as Thus, the condition (2.39), which the levels ± r should satisfy in order to attach a λ-realisation or a Yang-Baxter realisation without Wess-Zumino term to the sites (r, ±), becomes If one considers a model with N 2 λ-realisations, one has to impose N 2 relations as the one above, which form a linear system on the corresponding set of coefficients ν + r (or equivalently on the corresponding ν − r ). This is the advantage of this parametrisation, as one then has to solve linear constraints on the parameters instead of algebraic ones when using the zeroes. In particular, the solutions of these constraints are rational expressions of the remaining free parameters (however potentially quite complicated). This will be useful later in Subsection 4.1.2 when we will study the model with N coupled λ-models.
Coefficients ρ ± rs . Let us end this subsection by introducing some coefficients which will be useful to study the undeformed limit of the models in Subsection 3.4 and specific examples of models in Section 4. We define Using the expression (3.24) of the levels ± r , one shows that these coefficients can be rewritten as . (3.27) Using this expression, the operators U ± rs and V ± rs introduced in (3.6) and (3.13) can be re-expressed as (3.28)

Undeformed limit
As explained in Subsection 2.2.5, one can see the model constructed from N 1 inhomogeneous Yang-Baxter realisations and N 2 λ-realisations as a deformation of a simpler model, coupling together N 1 copies of the PCM with Wess-Zumino term and N 2 copies of the non-abelian T-dual of the PCM. This was understood by means of the so-called undeformed limit, in which the positions z + r and z − r collide for every pair of sites (r, ±), or equivalently by letting the parameters η r go to 0, while keeping r and [1] r finite (see Subsection 2.2.5 for details). The goal of this subsection will be to complete this discussion by studying this limit in the Lagrangian formulation of the model, focusing mostly on Yang-Baxter realisations (as explained in Subsection 2.2.5, the undeformed limit of λ-realisations requires a more subtle treatment which we will not detail here for conciseness). In particular, this will allow us to compare the methods and results presented in the previous subsections for deformed models to the ones presented in [23,24] for undeformed ones.
Interpolation formula. Let us focus for the moment on a pair of sites (r, ±), which we suppose to be associated with a Yang-Baxter realisation with Wess-Zumino term. The corresponding operators B (r) ± are then given by B (r) where δ r is defined in Equation (2.43), R (r) = Ad −1 g (r) • R r • Ad g (r) and Π (r) = Id − R (r) 2 /c 2 r . Recall that in Subsection 3.1, we found the Lagrangian expression of the Lax pair by interpolation methods, using the fact that one can express the Maurer-Cartan currents j (r) ± in terms of the evaluation of the Lax pair at the positions z ± r by Equation (3.1). In the present case, this equation can be rewritten as As recalled above, the undeformed limit corresponds to making the positions z ± r collide to the same point z r . It is then clear that in the undeformed limit, the above formula simply becomes j (r) ± = L ± (z r ). (3.30) This is precisely the interpolation formula obtained in [24,Equation (3.33)] for the model coupling N PCM with Wess-Zumino terms. In this reference, this formula plays a key role in obtaining the Lagrangian expression of the Lax pair of this model. The method developed in Subsection 3.1 of this article is thus a generalisation of the one of [24] to include deformed realisations. Recall from Equation (3.2) that in the deformed model, the currents J (r) ± are defined as the evaluations L ± (z ± r ). It is then clear from the above equation that in the undeformed limit, these currents J ± (see Equation (3.5)), becomes the identity in the undeformed limit, or equivalently, in components: For completeness, let us comment briefly on the homogeneous Yang-Baxter limit considered at the end of Subsection 2.2.5 (note that we considered the homogeneous limit only for realisations without Wess-Zumino term, in which case r = δ r = 0). Recall that this limit corresponds to taking the coefficient c r to 0. Recall also that the positions z ± r are given by z r ± c r η r . Thus, in the limit c r → 0, the equation (3.29) becomes j (r) ± = L ± (z r ) + η r R (r) L ± (z r ), (3.32) where L ± (z) denotes the derivative of L ± (z) with respect to the spectral parameter z. This is the equivalent of the equation (D.7) of [24], which was obtained when studying a model with N − 1 PCM realisations and one homogeneous Yang-Baxter realisation. It is interesting to compare the equations (3.29), (3.30) and (3.32): the undeformed interpolation formula (3.30) is corrected by a derivative term L ± (z r ) for an homogeneous Yang-Baxter deformation and by a finite difference term L ± (z r + c r η r ) − L ± (z r − c r η r ) /2c r for an inhomogeneous Yang-Baxter realisation.
Lagrangian expression of Y (r) . Recall from Subsection 3.2 that after the derivation of the Lagrangian expression of the Lax pair, the next step for performing the inverse Legendre transform of the model is to find the Lagrangian expression of the field Y (r) , which encodes the momentum fields of the model. This was done using Equation (3.11), which expresses Y (r) in terms of the Lax pair, through the operators C (r) ± . For a Yang-Baxter realisation, it can be re-written, after a few manipulations, as The undeformed limit correspond to taking η r to 0 while keeping [1] r and r = − r /2 finite. Recalling that z ± r = z r ± c r η r , the above equation then becomes in this limit This then coincides with the equation (3.36) of [24].
Recall from Subsection 3.2 that Equation (3.11) allows us to rewrite Y (r) in terms of the currents J (r) ± and the operators V ± rs , in Equation (3.12). In the undeformed limit, the currents J ± . Moreover, one can study the behaviour of the undeformed limit of the operators V ± rs using their expression (3.28). In particular, the coefficients ρ ± rs in this expression, defined by Equation (3.26), can be shown to converge in the undeformed limit to: with the coefficients ρ rs as defined in [24,Equation (3.40)]. Note that in this limit, the expression of the coefficient r also coincides with its expression in [24,Equation (3.38)]. Using the above limit of the coefficients ρ ± rs , as well as the expression (2.37) of the levels ± r in terms of the coefficients [0] r = −2 r and [1] r which stay finite in the undeformed limit, one can compute the limit of the operators V ± rs starting from their expression (3.28): In particular, Equation (3.12) agrees with [24,Equation (3.39)] in the undeformed limit: Action. Finally, we are now in a position to calculate the undeformed limit of the action of the model with N copies of the Yang-Baxter realisation. By reinserting the limits (3.31) and (3.34) in the expression (3.17) for the operator O, we find: Comparing to Equation (3.49) of [24], one sees that the action (3.16) then reduces to the one of N coupled copies of the PCM with Wess-Zumino term: S g (1) , · · · , g (N ) = dt dx Undeformed and q-deformed symmetries. The undeformed model (3.35) possesses N global symmetries acting by left translation on the fields g (r) : where h 1 , · · · , h N are constant elements of G 0 . Indeed, these transformations leave the Maurer-Cartan currents j (r) ± = g (r) −1 ∂ ± g (r) and the Wess-Zumino terms I WZ g (r) invariant. These global symmetries are broken by the introduction of deformations. Indeed, let us consider the model with N copies of the Yang-Baxter model studied in this subsection. The entries of the operators U ± and V ± are expressed in terms of the operators R (r) = Ad −1 It is a well-known result [6] (see also [45][46][47]) that in the Yang-Baxter model (with one copy and without Wess-Zumino term), this broken symmetry is in fact deformed into a q-deformed Poisson-Lie symmetry. Based on this result, it was explained in [24] that this is in general the case for every affine Gaudin model with a Yang-Baxter realisation (without Wess-Zumino term). In particular, the model coupling N copies of the Yang-Baxter models without Wess-Zumino term then possesses N q-deformed symmetries, which replace the translation symmetries (3.36). Their action on the fields of the model can be computed using the results of [48]: in particular, let us note that this action is non-local.
As the bilinear form κ is invariant under conjugacy, the undeformed model also possesses a global symmetry acting by simultaneous right translation on all the fields g (r) : with h a constant element of G 0 . As explained in [24], it corresponds to the diagonal symmetry of the underlying affine Gaudin model. As such, it is not broken by applying Yang-Baxter deformations to the various copies of the model. ± and the Wess-Zumino terms I WZ g (r) are invariant under this transformation, it is thus a symmetry of the deformed action (3.16). Note that a similar result holds for models involving λ-realisations: in this case, the corresponding fields g (r) should not transform by right multiplication but by conjugacy g (r) → h −1 g (r) h, while the fields corresponding to Yang-Baxter realisations still transform by right multiplication by h.

Yang-Baxter and λ-deformed coupled models
The action (3.16) presented in the previous section was obtained using the general ansatz introduced in Subsection 2.1.2 for the form of the Kac-Moody realisations defining the model. In this section, we specialise these results to the model constructed from N 1 copies of the Yang-Baxter realisation and N 2 copies of the λ-realisation. As we shall see, the particular form of these realisations will allow us to rewrite the action of this model in a simpler form. In particular, we will show that the integrable σ-model introduced in [29] corresponds to a particular limit of the model constructed from N copies of the λ-realisation. We will then focus on models with two copies and will rewrite their action in a more explicit form, using the expressions (3.22) and (3.23) of the inverse of U ± and of the operators O rs obtained in this case.

Deformed model with N 1 Yang-Baxter realisations and N 2 λ-realisations
Let us consider a model made up of N 1 copies of the Yang-Baxter realisation with Wess-Zumino term and N 2 copies of the λ-realisation. Let us now associate the former to the first N 1 pairs of sites (r, ±) and the latter to the last N 2 pairs. Then, from (2.19) and (2.20), one obtains, for the operators B . We observe that the relations C (r) ± respectively hold in the first and in the second case. Thus, from (3.6) and (3.13) and after a few manipulations, we obtain for the entries of the operator V ± : where the coefficients ρ ± rs have been defined in (3.26). From the expressions (3.17) and (3.18) of the operator O found in the previous section, we are now in a position to write the action of the model. We choose to express the entries O rs of this operator as in (3.17) for 1 ≤ r ≤ N 1 and as in the second equality in (3.18) for N 1 < r ≤ N 2 11 . Reinserting (4.1) in the form of the action (3.16), we obtain with α ± rs = ρ ± rs , 1 ≤ r ≤ N 1 , α + rs = − r δ rs , N 1 < r ≤ N and where S WZW, r [g (r) ] denotes the Wess-Zumino-Witten action of g (r) with level r :

Model with N Yang-Baxter realisations
Let us now briefly discuss the model with copies of the Yang-Baxter realisation only. In order to write its action in a compact form, let us introduce the G N 0 -valued field f = (g (1) , · · · , g (N ) ) and the block diagonal bilinear form κ N on g N 0 which restricts to κ on each factor g 0 . In this case, the action (4.2) gets simplified to where ± are operators on g N 0 which can be seen as N × N matrices with entries ( ± ) rs = ρ ± rs Id.
Let us describe more explicitly the operators U ± appearing in the action (4.3). From the expressions of the operators B (r) ± and C (r) ± for a Yang-Baxter realisation, one finds that and Π (r) = Id − R (r) 2 /c 2 r , and Let us end this subsection by presenting an alternative form of the action of the model. Let us introduce the operator c, with entries c rs = c r δ rs Id. Then, one can further rewrite the operator U ± as Finally, introducing ± = ± (Id − cϑ ± ) −1 , one can rewrite the action of the model in the form This way of writing the action of the model is quite similar to the way the action of the Yang-Baxter model with one copy is expressed and thus seems rather natural. Let us note however that it has some downsides compared to the expression (4.3). Indeed, the entries ρ ± rs and θ ± rs of the operators ± and ϑ ± appearing in the expression above are not straightforwardly expressed in terms of the parameters of the models (contrarily to the coefficients ρ ± rs and θ ± rs which were used in the previous formulation) as their definition involves the inversion of the operator Id − cϑ ± .
From the expression of the action above, one can simply check that its undeformed limit yields the action of N coupled PCMs with Wess-Zumino terms presented in [24]. Indeed, in this limit, the parameters θ ± rs and thus also the operators θ ± , go to zero. In particular, the coefficients ρ ± rs and ρ ± rs have the same limit. From Equation (3.33), we then see that in this limit, ρ + sr + ρ − rs → 2 ρ rs , with ρ rs as defined in [24].

Model with N λ-realisations
Action. Let us now discuss the case where we take λ-realisations only. For this model, using the notations introduced in the previous subsection, the action reads where K is an operator on g N 0 with entries K rs = r δ rs Id. From the expression of the operators B (r) ± of the λ-realisation, one can rewrite the operator U − as where the coefficients µ rs are defined as . (4.8) The action of the model then takes the simple form Parameters. Let us discuss what are the defining parameters of the model. We will use the parameterisation (z ± r , ν ± r , ∞ ) introduced in Subsection 3.3. As explained in this subsection, these parameters are convenient to take into account the fact that the levels ± r of the models should satisfy the constraints + r + − r = 0, that one has to impose to consider λ-realisations. Indeed, these constraints translate into the conditions (3.25) on the parameters z ± r and ν ± r . One can solve this condition by expressing the parameters ν + r in terms of z ± r and ν − r : (4.10) The remaining 3N + 1 parameters (z ± r , ν − r , ∞ ) are unconstrained: taking into account the translation and dilation redundancy among these parameters (see Subsection 2.2.4), the model is thus defined by 3N − 1 free parameters (for concreteness, one can for example fix this redundancy by fixing the values of ∞ and of one of the positions z ± r ). The coefficients µ rs defined in Equation (4.8) can be expressed in terms of this parametrisation as Similarly, the coefficient r appearing in the action (4.9) is given by where ν + r is replaced by its expression (4.10).
Comparison with [29]. Actions of the form (4.9) have been considered in [29] (and in [26][27][28] for the case N = 2, see Subsection 4.2.2). More precisely, the action (4.9) is identical to the action (2.13) of [29], with the matrix λ −1 in this reference identified in the present language with √ KM √ K −1 , or in components λ −1 rs = r / s µ rs . It was shown in [29] that the model defined by taking all entries of λ −1 to be zero except for λ −1 11 , · · · , λ −1 (N −1)1 and λ −1 N 2 , · · · , λ −1 N N is integrable. Let us now explain how this model can be obtained as a limit of the one constructed above by coupling together N λ-realisations. We introduce the following reparametrisation of the positions z ± r of the model: 12) in terms of new parameters y 1 , · · · , y N −1 , y 2 , · · · , y N and γ. We used here the translation redundancy on the parameters z ± r to fix the value of z − 1 to 0. Recall that one can also use the dilation redundancy to fix the value of ∞ : for future convenience, we choose here to fix it to Using this parametrisation, the model is then described by the 3N − 1 free parameters y 1 , · · · , y N −1 , y 2 , · · · , y N , ν − 1 , · · · , ν − N and γ. The limit we shall consider in this paragraph is γ → 0, while keeping the remaining parameters fixed.
Let us comment on the limit γ → 0 considered above. This limit consists in singling out two sets of positions Z 1 = {z − 1 , z + 1 , · · · , z + N −1 } and Z 2 = {z − 2 , · · · , z − N , z + N } and sending the distance between these two sets to infinity. It is thus quite similar to the decoupling procedure considered in [24,Subsection 2.3.3] 12 . According to this procedure, the sites (1, −), (1, +), · · · , (N − 1, +) corresponding to the positions Z 1 thus cease to interact with the sites (2, −), · · · , (N, −), (N, +) corresponding to the positions Z 2 in the limit γ → 0. This explains the structure of the model considered in [29], where the fields g (2) , · · · , g (N −1) have no interactions one with another. The theory before taking the limit γ → 0 then defines a non-trivial integrable generalisation of this model: indeed, although it corresponds to adding only one parameter, this introduces non-trivial interactions between all the different fields g (r) , as the coefficients µ rs then become generically all non-zero.
Following the decoupling procedure of [24], one describes the integrability of the model in the limit γ → 0 using two independent Lax pairs, which are obtained as two different limits of the initial Lax pair L ± (z). More precisely, let us consider: and L (2) It is clear that, before taking the limit γ → 0, both L ± (z) and L ± (z + γ −1 ) satisfy a zero curvature equation (as L ± (z) does) and thus still do after taking the limit. The reason behind the necessity of considering these two Lax pairs is that, loosely speaking, the Lax pair L ± (z) loses the information about the positions Z 2 in the limit γ → 0: the Lax pair L (1) ± (z) then only "corresponds to" the positions Z 1 (see [24] for a more precise treatment). Considering the shift of the spectral parameter by γ −1 , as done in the definition of L (2) ± (z), exchanges the roles of the sets Z 1 and Z 2 , so that the second Lax pair L (2) ± (z) contains the information about the positions Z 2 . This is coherent with [29], where the integrable truncation was described using two Lax pairs.
The Hamiltonian analysis of the corresponding Lax matrices was performed recently in [30], where it was shown that their Poisson brackets are described by twist functions. In the language of affine Gaudin models used above, these twist functions are obtained from the twist function ϕ(z) of the original model by a limit similar to the one of Equation (4.13) (see [24]): and One then finds that the twist function ϕ (1) has poles at the points {y 1 /ν − 1 , · · · , y N −1 /ν − 1 , 0} while the twist function ϕ (2) (z) has poles at the points {0, y 2 , · · · , y N }. Up to dilation and translation, these poles coincide with the ones obtained in [30].

Deformed models with two copies
Recall from Subsection 3.2 that in the case of a model with two copies, one can rewrite the operators U −1 ± and O rs more explicitly as in Equations (3.22) and (3.23). Using these results, we study in this subsection the models with two Yang-Baxter realisations and two λ-realisations.

Model with two Yang-Baxter realisations
Let us consider first the model with two Yang-Baxter realisations. In this case, we will use the first expression of the operators O rs in Equation (3.23). The entries of the operator U + can be read from (4.4) while the entries of V + are related to the ones of U + by equation (4.1a). Using the notationr introduced in Subsection 3.2, one then obtains the following expression for the operators O rs : where R (r) = c r Id + R (r) + c r δ r Π (r) , det(θ + ) = θ + 11 θ + 22 − θ + 12 θ + 21 and θ + rs is given by Equation (4.5).

Model with two λ-realisations
Let us now consider the model with two λ-realisations. Its action is given by Equation (4.6) with N = 2. Reinserting the explicit form (4.7) of the operator U − and calculating its inverse through (3.22), we find that in this case the operator KU −1 − appearing in the action is explicitly given by , with det(µ) = µ 11 µ 22 − µ 12 µ 21 and µ rs given by Equation (4.8).
Let us end this subsection by comparing this result with the ones of [26][27][28]. Indeed, the integrable σ-models introduced in these references can be obtained from the model above by taking limits similar to the one considered in Subsection 4.1.2 (which allowed us to compare the model with N copies of the λ-realisation with the integrable model of [29]).

Relation with 4d semi-holomorphic Chern-Simons theory
In this section, we explain how the models considered in this article can be obtained using the approach proposed recently by Costello and Yamazaki to generate integrable 2d field theories from 4d semiholomorphic Chern-Simons theory [37] (see [32][33][34][35][36][38][39][40] for additional references on this variant of Chern-Simons theory and its relation to integrable systems). Note that, in the terminology of [37], we restrict our attention here to 4d Chern-Simons theory with disorder defects. It was shown in [37] that the PCM with Wess-Zumino term and the integrable σ-model coupling N of its copies can be obtained from this approach. It was subsequently shown in [39] that the integrable 2d field theories obtained from 4d Chern-Simons theory with disorder defects are realisations of AGM. Moreover, it was explained in [40] how the Yang-Baxter model and the λ-model can also be derived following this approach. It is thus natural to search for a generalisation of these results for the AGM coupling together N 1 copies of the Yang-Baxter model and N 2 copies of the λ-model, which is the integrable field theory constructed in the present article.

4d semi-holomorphic Chern-Simons theory and integrable field theories
In this subsection, we will briefly sketch the method developed in [37] to generate integrable 2d field theories from 4d semi-holomorphic Chern-Simons theory. We will not explain this method in details here and mainly focus on the aspects that will be concretely relevant for the purpose of this article (we then refer to [37,40] for details). We will follow here the conventions of [40], which are in agreement with the ones used in the rest of this article.
4d Chern-Simons theory. The semi-holomorphic Chern-Simons theory is defined on the 4d manifold R × D × P 1 : the R × D part of this manifold corresponds to the 2d space-time with coordinates (t, x) of the resulting integrable field theory (here the spatial manifold D can be either the real line R or the circle S 1 , as in the rest of this article), while the Riemann sphere P 1 gives rise to the spectral parameter z of this integrable model. The 4d Chern-Simons theory is partly characterised by the choice of a meromorphic 1-form ω = ϕ(z)dz on P 1 : as shown in [39], the corresponding rational function ϕ(z) is the twist function of the resulting integrable model. The dynamical fields of the four-dimensional theory are the components A + , A − and Az of a g-valued gauge field A along the light-cone directions x ± of R×D and the anti-holomorphic directionz of P 1 (note that the component of A in the z-direction decouples from the theory and is not a physical degree of freedom). In addition to the choice of ω made above, the theory is then fully determined by specifying appropriate boundary conditions on the field A at the poles Z ⊂ P 1 of ω, i.e. at the poles of the twist function (see [37,40] and the next subsection for details). The action of the semi-holomorphic Chern-Simons theory is defined as [32] S where CS(A) is the standard Chern-Simons 3-form of A.
Parametrisation of the gauge field. In order to relate the 4d Chern-Simons theory to an integrable 2d model, one parametrises the gauge field components in the following form where g and L ± are fields respectively valued in the group G and the Lie algebra g. The equation of motion obtained by varying the action (5.1) with respect to Az then ensures that the fields L ± depend meromorphicaly on z. Moreover, the equations of motion obtained by varying A ± show that they also satisfy a zero curvature equation These two properties make the field L ± a good candidate for being the Lax pair of a 2d integrable model on R × D.
The fields of the 2d theory. Let us now explain how this integrable 2d field theory is constructed. For z in the Riemann sphere P 1 and a field φ on R × D × P 1 , we will denote by φ| z the field on R × D obtained by evaluating φ at the point z on the Riemann sphere. It is explained in [37,40] that for a point z ∈ P 1 \ Z which is not a pole of ω, the 2d field g| z can be set to a constant field equal to the identity of G by an appropriate gauge transformation on the gauge field A. The fact that we restrict here to points z which are not poles of ω is due to the fact that this gauge transformation on A should preserve the boundary conditions imposed on A at these poles and mentioned above (see [37,40] for details). Thus, the 2d fields g| z , z ∈ P 1 \ Z, are not physical degrees of freedom of the model. The dynamical fields of the 2d model we aim to construct are then defined to be the remaining degrees of freedom contained in g, i.e. its evaluations { g| z 0 } z 0 ∈Z at the poles of ω. Let us mention that in general, one should also consider the fields ∂ p z g| z 0 obtained by evaluating derivatives of g at the points z 0 ∈ Z: however, as explained in [37,40], for the boundary conditions considered in these references and that we shall consider in this article, these degrees of freedom can also be eliminated by gauge transformations. So far, we have considered only the degrees of freedom contained in the field g, which, as we see from Equation (5.2), encodes the component Az of the gauge field. Let us now consider the component A ± and thus the field L ± . As explained above, the equation of motion of Az ensures that L ± is meromorphic in z. In fact, it also implies that L ± can have poles in P 1 only at the zeroes of ω. This constrains quite strongly the dependence of L ± in terms of the variable z ∈ P 1 . Let us be more precise. As ω will ultimately be given by the twist function of the resulting 2d theory, let us denote its zeroes {ζ i } i∈{1,··· ,M } , in agreement with what was done in the rest of this article. These zeroes can be separated into two sets {ζ i } i∈I ± , labelled by subsets I + and I − of {1, · · · , M }, depending on which of the component L + or L − has a pole at ζ i (see [40] for details). This fixes the z-dependence of the fields L ± : more precisely, they are of the form for some 2d g-valued fields U i , U ∞ + and U ∞ − on R × D. In this equation, we have written the Lax pair as L ± (z) to stress its dependence on the spectral parameter z: note however that it also depends on the coordinates (t, x) ∈ R × D, through the 2d fields U i and U ∞ ± . Recall that the gauge field A obeys some specific boundary conditions at the poles z 0 ∈ Z of ω, which translate into conditions on the evaluations {L ± | z 0 } z 0 ∈Z and { g| z 0 } z 0 ∈Z . As observed in [37,40] and as we shall see in this article, these boundary conditions, combined with the z-dependence (5.3) of L ± , specify completely L ± in terms of the 2d fields { g| z 0 } z 0 ∈Z . The field L ± then does not contain any additional degrees of freedom and is interpreted as the Lax pair of the resulting 2d field theory on { g| z 0 } z 0 ∈Z (indeed, recall also from the previous paragraph that, on-shell, it satisfies a zero curvature equation on R × D).
Let us end this paragraph by the following remark. As argued above, the fields { g| z 0 } z 0 ∈Z describe all the degrees of freedom of the resulting 2d model. However, in general, these degrees of freedom are not all physical: there are some residual gauge symmetries acting on these fields, which depend on the type of boundary conditions considered. Moreover, there always exists an additional redundancy on these fields which consists on multiplying all of them on the right by an arbitrary G-valued field h on R × D (see [37,40]). This redundancy can be used to fix one of the fields { g| z 0 } z 0 ∈Z to the identity.
The effective 2d action. To complete the description of the 2d field theory obtained through this method, one has to describe its action. This is done by performing the integration over P 1 in the 4d action (5.1), resulting on an effective 2d action on R × D depending on the 2d fields { g| z 0 } z 0 ∈Z . However, we will not need the details of this procedure in the following and thus refer to [37,40] for details. In particular, it was shown in [40] that, for the type of boundary conditions that we shall consider in this article, this 2d action simply reads: where I WZ g| z 0 is the Wess-Zumino term of g| z 0 and j {z 0 } ± is defined as the Maurer-Cartan current

The models
Our aim in this section is to show explicitly that a certain class of 2d integrable field theories obtained using the Chern-Simons approach described in the previous subsection can be identified with the affine Gaudin models coupling together an arbitrary number of copies of inhomogeneous Yang-Baxter realisations and λ-realisations, as considered in the rest of this article. Let us then start by defining the particular class of 4d Chern-Simons theories that we shall consider here. As explained in [37,40] and recalled in the previous subsection, the 4d semi-holomorphic Chern-Simons theory is characterised by the choice of the meromorphic 1-form ω and of the boundary conditions on A at the poles Z of ω. Let us then define the 1-form and boundary conditions that we shall consider here.

1-form ω
Following [39] (see also the summary in the previous subsection), the meromorphic 1-form ω characterising the models obtained from the 4d Chern-Simons approach should coincide with ϕ(z)dz, where ϕ(z) is the twist function of these models when seen as realisations of AGM. As we aim to recover the models constructed in this article, we will then choose ω to be given by the twist function (2.26) considered in the previous sections, i.e.
This 1-form has 2N simple poles at the points z ± r and a double pole at ∞. In the language of the previous subsection, one then has Z = {z + 1 , z − 1 , · · · , z + N , z − N , ∞}. Following the notations of this article, let us define ± r as the residues of ω at the poles z ± r , which coincide with the levels of the model when seen as a realisation of AGM.
λ-boundary condition. Let us describe the second type of boundary condition at the pair of simple poles z ± r that we shall consider, which we call the λ-boundary condition. It can be imposed only if the poles z ± r and the residues ± r are real and satisfy the additional condition + r + − r = 0 (note that this is identical to the condition (2.39) that one should impose to consider a λ-realisation in an affine Gaudin model). The λ-boundary condition is then simply given by For a λ-boundary condition, we define the parameter r = − r /2 = − + r /2, which is equal to the Wess-Zumino coefficent r defined for a λ-realisation (see Subsection 2.1.3).

Fields of the model
Let us consider a 4d Chern-Simons theory with ω as in Equation (5.5) and with N 1 Yang-Baxter boundary conditions and N 2 λ-boundary conditions. Let us describe what are the dynamical fields of this model. As recalled in the previous subsection, these fields are given by the evaluations g| z 0 of the field g at the poles z 0 ∈ Z of ω and thus by the 2N + 1 fields g| ∞ , g| z + r and g| z − r . However, as mentionned in the previous subsection and explained in [37,40], we can eliminate some of these degrees of freedom. In particular, recall from the previous subsection that we can fix one of the fields g| z 0 to the identity: here, we will choose to fix the field at infinity g| ∞ . Moreover, as explained in [40], if one considers a Yang-Baxter boundary condition or a λ-boundary condition at the pair of poles z ± r , there exists a residual gauge symmetry on the fields g| z + r and g| z − r . In the case of a Yang-Baxter boundary condition, this gauge symmetry can be fixed by imposing g| z + r = g| z − r : we then define g (r) as their common value. For the λ-boundary condition, it can be fixed instead by imposing g| z − r = Id: we then define g (r) = g| z + r . To summarise, the fields of the model are the N group-valued fields g (1) , · · · , g (N ) and we have g| ∞ = Id, YB-BC: g| z + r = g| z − r = g (r) , λ-BC: g| z + r = g (r) , g| z − r = Id.

Identification of the two approaches
Let us consider the 2d integrable field theory defined in the previous subsection with N 1 Yang-Baxter boundary conditions and N 2 λ-boundary conditions. We will prove in this subsection that it can be identified with the AGM with N 1 Yang-Baxter realisations and N 2 λ-realisations studied in the previous sections. In order to do so, we shall show that the two approaches lead to the same Lax pair as well as the same action.

Identification of the Lax pairs
Let us consider the Lax pair of the model coming from 4d Chern-Simons theory as given by Equation (5.3). Let us now express it in terms of the fields g (r) of the model, using the boundary conditions that are imposed on the gauge field A at the poles z 0 ∈ Z of ω.
Pole at infinity. Let us start with the pole at z 0 = ∞, for which the boundary condition is simply defined by Equation (5.6). From the fact that g| ∞ = Id (see Subsection 5.2.3) and the expression (5.3) of the Lax pair L ± (z), it is clear that the evaluation of the gauge field (5.2) at z = ∞ gives Combining this with the boundary condition (5.6), we then get that U ∞ ± = 0.
Pair of poles with Yang-Baxter boundary condition. Let us know consider a pair of simple poles z ± r and let us suppose that we imposed on this pair a Yang-Baxter boundary condition (5.7). As explained in Subsection 5.2.3, in this case, we have g| z + r = g| z − r = g (r) . Thus, the evaluation of the gauge field (5.2) at z ε r , for ε ∈ {+1, −1}, is given by where j (r) are the Maurer-Cartan currents of the field g (r) . After a few manipulations, the Yang-Baxter boundary condition (5.7) then becomes j (r) Noting that R r is skew-symmetric and Π r is symmetric, this can be rewritten as j (r) . The operators B (r) ± found here coincide exactly with the operators, denoted in the same way in the rest of this article, coming from a Yang-Baxter realisation (see Subsection 2.1.3). The above equation is then equivalent to the equation (3.1) obtained in the context of affine Gaudin models.
Pair of poles with λ-boundary condition. Let us now consider a pair of simple poles z ± r with the λ-boundary condition (5.8). In this case, we have g| z + r = g (r) and g| z − r = Id (see Subsection 5.2.3). Thus, the evaluations of the gauge field (5.2) at z + r and z − r read Similarly to what was done in the previous paragraph for the Yang-Baxter boundary condition, the λ-boundary condition (5.8) can then be rewritten j (r) .

The operators B
(r) ± coincide with the ones introduced in the previous sections for a λ-realisation (see Subsection 2.1.3). As for the Yang-Baxter boundary condition, we then recover the equation (3.1) obtained through the affine Gaudin model approach.
Summary. Let us summarise the results of this subsection. We have proved from the boundary condition at z = ∞ that the fields U ∞ ± vanish. The component L ± (z) of the Lax pair (5.3) has then no constant term and has simple poles at the zeroes {ζ i } i∈I ± . Thus, it has the same meromorphic z-dependence as the Lax pair (2.35) of the corresponding affine Gaudin model. Moreover, we showed that the boundary conditions imposed at the pairs of simple poles z ± r in the Chern-Simons approach coincide exactly with the Equation (3.1) obtained in the affine Gaudin model approach. Recall from Subsection 3.1 that this equation, combined with the meromorphic z-dependence mentioned above, allowed us to express the Lax pair L ± (z) in terms of the Maurer-Cartan currents j (r) ± by means of interpolation techniques. This ensures that the Lax pairs obtained from the two approaches can be identified.

Identification of the actions
Let us end this section by showing that the action obtained by the Chern-Simons approach for the model with N 1 Yang-Baxter and N 2 λ-boundary conditions coincides with the one of the AGM with N 1 Yang-Baxter and N 2 λ-realisations, computed in Section 3. As recalled in Subsection 5.1, the former is given by Equation (5.4). Since we proved in the previous subsection that the Lax pair L ± (z) of the two models coincide, one can then re-insert in this equation the expression (3.3) of L ± (z) obtained in the AGM approach using interpolation techniques. As the twist function has simple poles at z ± r with residues ± r , we then get Moreover, recall that the field g| ∞ has been set to the identity. The action (5.4) then becomes Recall from Subsection 5.2.3 that the fields g| z ± r are related to the fundamental fields g (r) of the model, in a way which depends on the type of boundary conditions considered at the poles z ± r . Equation (5.9) then expresses the action of the model in terms of the Maurer-Cartan currents j (r) ± , the currents J (r) ± and the Wess-Zumino terms of the fields g (r) . In the AGM approach, we obtained a similar expression for the action in Equation (3.15). In the rest of this subsection, we shall show that these two expressions coincide, thus proving that the models obtained from the 4d Chern-Simons and the AGM approaches are identical. For that, we will prove that for every r ∈ {1, · · · , N }, we have In order to show these identities, one needs to distinguish the cases where the pair of poles z ± r is associated with a Yang-Baxter boundary condition or a λ-boundary condition in the Chern-Simons approach and, correspondingly, with a Yang-Baxter realisation or a λ-realisation in the AGM approach.
Yang-Baxter boundary condition. Let us start with the Yang-Baxter boundary condition. In this case, recall that g| z + r = g| z − r = g (r) and that we defined the Wess-Zumino coefficient to be r = −( + r + − r )/2. The Wess-Zumino terms in Equation (5.9) corresponding to these poles thus satisfy Equation (5.10). Let us now focus on the term Υ r . Note first that j ± . This implies that the operators U ± rs and V ± rs defined in Equations (3.6) and (3.13) satisfy V ± rs = ρ ± rs Id ± r 2 U ± rs .
Using the expression (3.7) of J Re-inserting this identity in the above expression for Υ r then shows that it satisfies Equation (5.11), as required.
λ-boundary condition. Let us consider now a pair of poles z ± r associated with a λ-boundary condition. One then has g| z + r = g (r) and g| z − r = Id (see Subsection 5.2.3). Recall moreover from Subsection 5.2.2 that the Wess-Zumino coefficient is defined for λ-boundary conditions as r = − + r /2. Thus, the Wess-Zumino terms corresponding to the poles z ± r in the action (5.9) are given by Equation Reinserting Equations (5.13) and (5.14) in the expression (5.12) of Υ r , one sees that Υ r satisfies Equation (5.11), as required.

Conclusion and perspectives
In this article, we constructed integrable deformations of the coupled σ-model introduced in [23], using the formalism of affine Gaudin models. In particular, we obtained explicit expressions of the action and Lax pair of the deformed models corresponding to arbitrary combinations of Yang-Baxter and λ-deformations. Moreover, we showed that the integrable coupled λ-models introduced recently in [26][27][28][29] can be seen as particular limits of the models constructed here. Let us now conclude by discussing some perspectives of the present work.
As explained in Subsection 3.4, the deformed models constructed in this article which involve a Yang-Baxter realisation without Wess-Zumino term possess a corresponding q-deformed Poisson-Lie symmetry, which replaces the left translation symmetry of the undeformed model. It is well known that the Yang-Baxter model (with one copy and without Wess-Zumino term) in fact possesses a larger (infinite) symmetry algebra, satisfying the relations of an affine q-deformed Poisson algebra [49] (see also [46,47,50]), which replaces the Yangian symmetry of the undeformed Principal Chiral Model [51].
It would be interesting to understand whether such infinite extensions of the q-deformed symmetries also exist for the deformed coupled models and what would be their underlying algebraic structure.
The integrable deformed models constructed in this article still possess an undeformed symmetry, corresponding to the diagonal symmetry of the underlying affine Gaudin model, which acts on the fields g (r) by right multiplication (g (r) → g (r) h) or conjugacy (g (r) → h −1 g (r) h), depending on whether the realisations at sites (r, ±) are Yang-Baxter realisations or λ-realisations. It was explained in [25] that for a general realisation of affine Gaudin model of the type considered in [24], one can construct an integrable Yang-Baxter deformation which breaks its diagonal symmetry. Thus, one can introduce a further integrable deformation of the deformed coupled σ-models constructed in this article. As explained in [25], this deformation procedure involves gauging the model and thus requires treating Hamiltonian first-class constraints. For brevity, we chose not to treat these deformations in the present article: however, we expect that they can be studied using similar methods to the ones developed here. For the case with one copy only, it was conjectured in [25] that these further deformed theories coincide with already known models, namely the bi-Yang-Baxter model (see [4,52,53] for the case without Wess-Zumino term and [54] for the case with Wess-Zumino term) and the generalised λ-model [55].
It is known that the Yang-Baxter and λ-models are Poisson-Lie T-dual [56][57][58] to one another [31,55,59,60], while the Yang-Baxter model with Wess-Zumino term is Poisson-Lie T-dual to itself with different parameters [61]. It would be interesting to investigate the various possible dualities between the coupled models constructed in this article and how they would manifest themselves in the underlying geometry of their target space G N 0 . The study of Poisson-Lie T-dualities between deformed σ-models led to their reformulation as E-models [60,62,63], making their duality properties manifest. A natural direction to explore is thus to search for a similar reformulation of the coupled models constructed here.
The results of Section 5 illustrate once again the deep relation existing between the approaches to two-dimensional integrable field theories from affine Gaudin models [22] and from four-dimensional semi-holomorphic Chern-Simons theory [37], first established in [39] and further supported in [40]. In particular, the analysis conducted in this section strengthens the apparent correspondence between the choice of realisations in the first approach and the choice of boundary conditions in the second one. It would be interesting to understand in more details this correspondence. Let us also note that the construction of the Yang-Baxter and λ-boundary conditions introduced in [40], which uses isotropic subalgebras of the complex or real double of g 0 , is reminiscent of the structure underlying Poisson-Lie T-duality and E-models and could thus provide interesting directions for investigating the questions raised in the previous paragraph.
It would also be interesting to explore the quantum properties of these classically integrable deformed σ-models. For example, a natural question is whether these models are one-loop renormalisable and if there exist conformal fixed points in this space of models. The results obtained in [27][28][29] about the renormalisation of the coupled λ-models introduced in these references (which are limits of the models considered here), already show a rich structure in their renormalisation group flow. As a further possible step, it would be interesting to investigate the higher-loops renormalisability of these models and, if needed, the corresponding quantum corrections of their underlying geometry, as recently studied in [64][65][66] for non-coupled models.

Appendices
A Proof of the identities (2.11) In this appendix we will present the calculation of the non-ultralocal terms (i.e. terms containing derivatives of the delta distribution) in the bracket (2.9), using the ansatz (2.10) for the currents J ± in terms of the operators B ± and C ± . In particular, we will show that this computation implies that these operators satisfy the identities (2.11). Let us start by noting that in order to perform this computation, we need the Poisson brackets between the following objects: B ± , Y , C ± and j. However, let us recall that we have assumed the operators B ± and C ± to depend only on the field g (and not on its derivative ∂ x g). Thus, the non-ultralocal terms in the brackets of J ± can only come from the brackets between the fields Y and j. More precisely, for , σ ∈ {±}, we have Comparing with Equation (A.1), we then see that the operators B ± and C ± should satisfy the identities (2.11).
B Simplification of the action (3.14) In this appendix, we show that the non-Lorentz invariant terms appearing in the second line of the action (3.14) cancel with the term in the first line containing the Hamiltonian. For that, let us start by computing the expression of the Hamiltonian in terms of Lagrangian fields.
Hamiltonian in terms of Lagrangian fields. We will proceed here in a similar fashion to what has been done in [24]. Let us start by noting that, combining the equations (2.27) and (2.28), the Hamiltonian can be rewritten as where we have used the fact that i = ±1 for i ∈ I ± . We then need to look for the Lagrangian expression of the quantities Γ(ζ i ). This is done by relating them to residues of the Lax pair. More precisely, let us fix i ∈ I ± : from (2.35), we have Γ(ζ i ) = ± 1 2 ϕ (ζ i ) res z=ζ i L ± (z).