3-dimensional mixed BF theory and Hitchin's integrable system

The affine Gaudin model, associated with an untwisted affine Kac-Moody algebra, is known to arise from a certain gauge fixing of 4-dimensional mixed topological-holomorphic Chern-Simons theory in the Hamiltonian framework. We show that the finite Gaudin model, associated with a finite-dimensional semisimple Lie algebra, or more generally the tamely ramified Hitchin system on an arbitrary Riemann surface, can likewise be obtained from a similar gauge fixing of 3-dimensional mixed BF theory in the Hamiltonian framework.


Introduction
The Heisenberg spin chain can be obtained from 4-dimensional mixed topologicalholomorphic Chern-Simons theory on R 2 ×C by introducing certain line defects along the topological plane R 2 for each site of the spin chain [C1,C2,CWY1,CWY2]. This elegant description of the Heisenberg spin chain is ultimately possible because the integrable structure of the latter is underpinned by the quantum Yang-Baxter equation R 12 (z 1 , z 2 )R 13 (z 1 , z 3 )R 23 (z 2 , z 3 ) = R 23 (z 2 , z 3 )R 13 (z 1 , z 3 )R 12 (z 1 , z 2 ) z 2 z 3 z 1 = z 2 z 3 z 1 which is neatly encoded in the semi-topological nature of the 4-dimensional Chern-Simons theory. By contrast, the integrability of the (classical and quantum) Gaudin model [G], or more generally the Hitchin system [H], is underpinned by the classical Yang-Baxter equation [r 12 (z 1 , z 2 ), r 13 (z 1 , z 3 )] = [r 23 (z 2 , z 3 ), r 12 (z 1 , z 2 )] − [r 32 (z 3 , z 2 ), r 13 (z 1 , z 3 )], whose topological origin is less clear. On the other hand, affine Gaudin models, whose integrability is also underpinned by the classical Yang-Baxter equation (1), see [V1], can be obtained [V2] from the same 4-dimensional mixed topological-holomorphic Chern-Simons theory [CY] on R 2 × CP 1 , this time by introducing surface defects along R 2 placed at the marked points z i ∈ C of the affine Gaudin model. A natural question is therefore whether the ordinary Gaudin model, associated with a finite-dimensional semisimple Lie algebra g rather than an affine Kac-Moody algebra, can also be described in a similar way.
The purpose of this paper is to show that the tamely ramified Hitchin system on a Riemann surface C, in which the Higgs field has simple poles at certain marked points on C, can be described using a collection of various line defects in 3-dimensional mixed (topological-holomorphic) BF theory on R × C. In particular, the Gaudin model is obtained as the special case when C = CP 1 . More precisely, we perform a Hamiltonian analysis of the 3d mixed BF theory, whose fields are a partial connection 1-form A and a (1, 0)-form B, with suitably chosen line defects. Using the condition Az = 0 to fix the gauge invariance, we find that the dynamics on the reduced phase space coincides with that of the tamely ramified Hitchin system. In particular, the (1, 0)-form B becomes meromorphic and gets identified with the Higgs field. This is completely analogous to the relationship found in [V2] between 4d mixed topologicalholomorphic Chern-Simons theory on Σ×CP 1 , with the cylinder Σ = R×S 1 , and the affine Gaudin model. In other words, our analysis shows that 3d mixed BF theory is to the Gaudin model what 4d Chern-Simons theory is to the affine Gaudin model. The plan of the paper is as follows.
In §2 we recall the definition of 3-dimensional mixed topological-holomorphic BF theory on R × C for some Riemann surface C, depending on a g-valued (1, 0)-form B and a g-valued connection 1-form A. In §2.3 we introduce two types of line defects, which we refer to as type A and B, respectively. Type A defects ensure that in the gauge Az = 0 the field B is meromorphic on-shell, with poles at the location of the defects. This field then gets identified with the Lax matrix L of the Gaudin model, or the Higgs field Φ of the tamely ramified Hitchin system. Adding the type B defect ensures that (the negative of) the time component of the connection 1-form A gets correctly identified with the matrix M in the Lax equation ∂ t L = [M, L]. We also derive the 1d action for the Gaudin model along the lines of [DLMV3,CSV].
In §3 we perform a Hamiltonian analysis of 3-dimensional mixed topologicalholomorphic BF theory with the two types of defects introduced in §2.3. There is a first class constraintμ generating the gauge symmetry of the theory. In the absence of type A defects,μ coincides with the moment map µ =∂ A Φ of the Hitchin system, after identifying the Higgs field Φ with the (1, 0)-form B of the 3-dimensional mixed topological-holomorphic BF theory. Upon imposing the condition Az = 0 to fix the gauge symmetry generated byμ, the component B z of the (1, 0)-form B becomes meromorphic with poles at the location z i of the type A line defects and is shown to satisfy the Lax algebra with respect to the Dirac bracket. Morever, the Hamiltonian on the reduced phrase space coincides with that of the Hitchin system.
We end the paper with a brief discussion of possible future directions in §4.
2 3d mixed BF theory on R × C Let G be a semisimple Lie group over C, with Lie algebra g and fix a non-degenerate invariant symmetric bilinear form ⟨⋅, ⋅⟩ ∶ g × g → C on g. Let C be a Riemann surface. We shall consider the 3-dimensional classical mixed topological-holomorphic BF theory on R×C, or 3d mixed BF theory for short -see e.g. [GW, GRW, Ra] where the theory is discussed using the BV formalism. The field content of this theory consists of a g-valued (1, 0)-form B on R × C, where the bigrading (p, q) corresponds to the one induced by the complex structure on C, together with a g-valued connection 1-form A on R × C. We denote the curvature of the latter as F (A) = dA + 1 2 [A, A]. The action of 3d mixed BF theory is then given by

Gauge invariance
The 3d mixed BF action (2) is trivially invariant under gauge transformations of the form A → A + χ for any g-valued (1, 0)-form χ on R × C. Indeed, χ drops out from the action since B is a (1, 0)-form. We can fix this invariance by requiring that A has no (1, 0)-component so that it locally takes the form A = Azdz + A t dt for some local coordinate t on R and a local holomorphic coordinate z on C. From now on we will always take A to be a partial connection of this form. More interestingly, the action (2) is invariant under the action of any G-valued function g on R × C acting by gauge transformations on the connection 1-form A and by conjugation on the field B, namely where∂ is the Dolbeault differential on C, i.e. the (0, 1)-component of the de Rham differential d C = ∂ +∂ on C, and d R denotes the de Rham differential on R. Indeed, under such a transformation the curvature 2-form F (A) transforms by conjugation F ( g A) = gF (A)g −1 and so the invariance of the action follows from the adjoint G-invariance of the bilinear form ⟨⋅, ⋅⟩.

Equations of motion
To derive the equations of motion we consider the variations B → B +ǫ and A → A+η by an arbitrary (1, 0)-form ǫ and 1-form η on R × C. Varying the action we find where in the third equality we have used the fact that A and η are both g-valued 1-forms, so that [η, A] = [A, η]. In the last equality we used Stokes's theorem, noting that ⟨B, dη⟩ = ⟨dB, η⟩ − d⟨B, η⟩, and the adjoint invariance of the bilinear form. The equation of motion for B is therefore F (A) = 0, or explicitlȳ while the equation of motion for A reads Letting z be a local holomorphic coordinate on C and t a global coordinate on R, and writing the two fields in components as B = B z dz and A = Azdz + A t dt, we can write the equations of motion (4) more explicitly in components as The first key observation to make here is that the equation of motion (5c) is very reminiscent of the Lax equation However, to make this superficial resemblance more precise we would need B z and −A t to both be holomorphic (or more generally meromorphic) in order to identify them with the Lax pair L and M of an integrable system. The second observation, based on the other two equations of motion (5a) and (5b), is that this can be achieved by working in the gauge where Az = 0. Indeed, in this gauge the two equations (5a) and (5b) reduce to ∂zA t = 0 and ∂zB z = 0, respectively, which express the fact that A t and B z are both holomorphic on C.

Introducing line defects
In the Lax equation (6) of an integrable system, however, L and M are more generally g-valued meromorphic functions with poles at certain marked points. This is, in fact, necessary if C has genus zero, i.e. when C = CP 1 . Moreover, as it stands there is no relation between B z and −A t in (5c), while in (6) the matrix M is typically built out of the Lax matrix L. We can fix both of these issues by introducing two different types of line defects in the action (2). We will refer to these as type A and type B line defects, since these will depend on the fields A and B, respectively.

Type A line defects
The Lax pair of the Gaudin model is formed of two g-valued meromorphic functions L and M on CP 1 with L having poles at certain marked points z i ∈ C for i = 1, . . . , N . In order to view B z and −A t as such a Lax pair, but working on a more general Riemann surface C, we would like them to be meromorphic instead of holomorphic, with B z having poles at certain marked points z i ∈ C. To this end, we pick and fix elements u i ∈ g and introduce G-valued fields h i on R for i = 1, . . . , N . Following [CY], see also [CSV], we add to the action (2) the following sum of line defects where ι z i ∶ R × {z i } ↪ R × C is the embedding of the line defect at z i into the total space. In particular, the pullback ι * z i A is just the evaluation of the component A t dt at the point z i ∈ C so that we can rewrite the defect action (7) more explicitly as In order to maintain the gauge invariance of the action (2) after adding (7) to it, we should require that the latter is itself gauge invariant. This can easily be achieved by supplementing the gauge transformations (3) of the fields A and B by the transformation for the G-valued fields h i , i = 1, . . . , N . Consider now the extended actioñ Since the defect action (7) does not depend on B, the equations of motion (4a) for B are unchanged. On the other hand, the equation of motion (5b) for A in a local holomorphic coordinate z on an open neighbourhood U of the point z i is now replaced by where we introducedû i ∶= h i u i h −1 i for each i = 1, . . . , N and δ zz i denotes the Dirac δ-distribution at the marked point z i with the property that for any function f ∶ U → C on the neighbourhood U ⊂ C of z i equipped with the local holomorphic coordinate z.
In the gauge Az = 0, the modified equation of motion (10) reads Using the fact that ∂z(z − z i ) −1 = −2πiδ zz i we may rewrite this equation as which tells us that B z has a simple pole at z i with residueû i there, i.e.
where O(1) denotes terms which are holomorphic at the point z i .
When C = CP 1 , corresponding to the Gaudin model, if we fix a global coordinate z on C ⊂ CP 1 and require B to have a simple pole also at infinity then we can explicitly write B as the g-valued meromorphic (1, 0)-form By varying the action (9) with respect to h i → e ǫ i h i for some g-valued function ǫ i on R we find N further equations of motion for i = 1, . . . , N . But given the meromorphic behaviour (12) of the (1, 0)-form B at each of the marked points z i , these are merely consequences of the equation of motion (5c) given by taking the residue at each z i , assuming that A t is regular at z i , as will be the case in §2.3.2.

Type B line defects
The type A line defects introduced in §2.3.1 ensured that B z is no longer holomorphic in the gauge Az = 0 but rather meromorphic with poles at certain marked points z i ∈ C. The type B line defects will have a similar effect on the field A t . However, since −A t is meant to play the role of M in the Lax pair (6), we want it to be built out of B z , which plays the role of the Lax matrix L. Let P ∶ g → C be a G-invariant polynomial on g and fix a point w ∈ C distinct from the marked points z i ∈ C for i = 1, . . . , N at which the type A line defects were inserted in §2.3.1. We consider the following line defect where z is a local holomorphic coordinate around the point w ∈ C and, writing B = B z dz in this coordinate, B z (w) denotes the evaluation of B z at the point w.
The G-invariance of the polynomial P ensures that the action (15) is gauge invariant. Therefore, adding it to the gauge invariant action (9) obtained so far, we obtain the full gauge invariant action Since the defect term (15) only depends on B, it will not modify the equations of motion for A. Only the equation of motion for B, namely (5a) which has so far remained unchanged, will be modified. To derive it we note that the variation of the defect action (15), under the variation B → B + ǫ considered in §2.2 with ǫ = ǫ z dz in the local holomorphic coordinate z, reads where in the third line we introduced the element P ′ (B z (w)) ∈ g such that the linear map ⟨P ′ (B z (w)), ⋅⟩ ∶ g → C is the derivative of P ∶ g → C at B z (w) and kept only the terms linear in ǫ z (w). It follows that (5a) is now replaced by In the gauge Az = 0 this simplifies to or in other words, In the case C = CP 1 this tells us that the expression in brackets is a constant. Taking this contant to be zero we therefore obtain which coincides with the usual expression for M = −A t in terms of L = B z , see for instance [BBT,(3.33)] in the case when g = gl r and the polynomial P ∶ gl r → C is given by X ↦ tr(X n ) for some n ∈ Z ≥1 . Indeed, in this case we have P ′ (X) = nX n−1 for any X ∈ gl r so that (19) becomes In connection with the Hamiltonian analysis to be performed in §3 below, where the classical r-matrix r 12 (z, w) = C12 w−z will be introduced in (43), note that we can rewrite (20) in the more recognisable form Substituting the expression (19) for A t into the equation of motion (5c) we obtain the desired Lax equation where we have explicitly written the dependence of B z on the spectral parameters. We thus expect from the general theory of integrable systems, see for instance the Proposition [BBT,p.47], that the time coordinate t along the topological direction of the 3-dimensional space R × C is identified, through the introduction of the type B defect (15), with the time induced by the Hamiltonian To confirm this we will move to the Hamiltonian formalism in §3 below.

Unifying 1-dimensional action
We have now shown that the gauge fixed equations of motion for 3d mixed BF theory in the presence of type A and B defects correspond exactly to the Lax equation (21) of the Gaudin model with Lax matrix L(z) = B z (z) given by (13), where the residueŝ u i = h i u i h −1 i are coadjoint orbits through the fixed elements u i ∈ g and parametrised by the dynamical G-valued variables h i ∈ G.
At this stage it is therefore natural to proceed along the lines of [DLMV3], where a unifying 2d action for integrable field theories of affine Gaudin type was derived from the 4d Chern-Simons action of [CY]. In a similar spirit, in the present context we would like to obtain a 1d action for the Gaudin model with Lax matrix (13) starting from the 3d mixed BF theory with both type A and type B defects. In fact, the procedure followed in [CSV] is closer in spirit to the present case since we do not have to deal with the presence of a meromorphic 1-form ω having zeroes, as in the 4d Chern-Simons action considered in [DLMV3].
Following [CSV], we will therefore substitute the solutions to the equations of motion (10) and (17) (but crucially not (5c)) in the gauge Az = 0, namely (13) and (19) respectively, into the full action (16). We will do this for the three pieces in the action separately. For the bulk action (2) we find where in the second equality we used the gauge Az = 0. In the last equality we used the fact that B z is identified with the Lax matrix L together with the identity (18), and then performed the integral over C using the presence of the δ-function.
For the type A defect action (7) we have where in the second equality we used (19) evaluated at z = z i and in the last line we recognised the sum over i in the second term as the expression for the Lax matrix L(w) = B z (w) in (13). Note that the second term on the right hand side exactly cancels the expression found above for the bulk action S 3d [A, B]. Finally, the type B defect action (15) is simply S B−def [B] = − ∫ R H P w dt using the expression (22) for the Hamiltonian alluded to in §2.3.2 and to be confirmed in §3. Putting all the above together, we deduce that the full action (16) reduces to the simple form where we have suppressed the dependence on the fields A and B since these have now been expressed in terms of the dynamical variables h i ∈ G and the fixed elements u i ∈ g for i = 1, . . . , N . We recognise (23) as the first order action associated with the Hamiltonian H P w in (22) but where the conjugate momentum X i ∈ g of h i ∈ G has been fixed to the constant element X i = −u i . This is consistent with the Hamiltonian analysis to be performed in the next section. Namely, we will find in §3.1.2 that there is a primary constraint X i + u i ≈ 0 on the conjugate momentum X i ∈ g of the dynamical variable h i ∈ G.
We can check directly that the equations of motion of the 1d action (23) are given by (14), with A t as in (19), by varying it with respect to h i → e ǫ i h i for some arbitrary g-valued variable ǫ i . Under this variation, the Lax matrix L(w) transforms to Hence, using the explicit expression H = P (L(w)) for the Hamiltonian, the variation of the action is given by where in the last equality we have used Stokes's theorem and the adjoint invariance of the bilinear form. The N equations of motion for the h i are therefore and using (19), we do indeed recover the equations of motion (14) found previously from adding in type A defects in §2.3.1.

Hamiltonian analysis
Throughout this section we shall work in some local coordinate z on some open subset of C. Our starting point is the Lagrangian density of the action (16) written in terms of the components of the g-valued bulk fields A = Azdz +A t dt and B = B z dz and in terms of the G-valued defect variables h i for i = 1, . . . , N which we introduced at the type A defects in §2.3.1, namely

Conjugate momenta and primary constraints
To move to the Hamiltonian formalism we first determine the conjugate momenta of the bulk fields Az, A t and B z and the defect variables h i . We shall find various primary constraints, some of which will be second class. We shall impose the latter strongly at this stage by introducing a corresponding Dirac bracket. To alleviate the notation, all Dirac brackets computed in this section will ultimately be renamed simply as {⋅, ⋅} before moving on to §3.2 where we work out the secondary constraints. We will begin in §3.1.1 by considering the conjugate momenta of the bulk fields Az, A t and B z , as the conjugate momenta to the G-valued defect variables h i will need to be handled with more care, as discussed in §3.1.2 below.

Bulk canonical fields
The conjugate momenta to the g-valued bulk fields Az, A t and B z are the g-valued fields given respectively by which satisfy the canonical Poisson bracket relations {Π t1 (z), A t2 (z ′ )} = C 12 δ zz ′ , We have three primary constraints associated with the bulk fields, namely The first is clearly first class and the latter two are second class with Poisson bracket We set these both strongly to zero immediately by introducing the corresponding Dirac bracket, with respect to which we still have the same relations between the remaining fields, i.e.
and hence, by an abuse of notation, we will continue to denote this Diract bracket as {⋅, ⋅}.

Defect canonical variables
We have yet to find the conjugate momenta to the G-valued variables h i , i = 1, . . . , N introduced at the type A defects. This can be done by working in local coordinates φ α on the group G where α ranges from 1 to dim G, the dimension of G. We refer the reader, for instance, to [L1, §3.1.2] for details. Each variable h i ∈ G can then be described locally in terms of the dim G variables φ α i ∶= φ α (h i ). The relevant part of the Lagrangian in finding the conjugate momenta is which in the second expression we have rewritten in terms of the local coordinates φ α i , where ∂ α denotes the partial derivative with respect to the coordinate φ α . The corresponding conjugate momenta are therefore given by and these have the usual canonical Poisson bracket relations To return to a coordinate free description of the phase space, we define a matrix L a α for some fixed basis {I a } of g such that This L a α is invertable and we denote the inverse as L α a following the conventions of [L1, §3.1.2]. We can then introduce a coordinate-free g-valued variable X i which encodes the conjugate momentum π i,α as where {I a } is the basis of g dual to {I a } with respect to the bilinear form ⟨⋅, ⋅⟩. We therefore have a coordinate free description of the phase space, with the canonical Poisson brackets in local coordinates being equivalent to for each i, j = 1, . . . , N . Using the definition of the matrix L a α in (28) we have It then follows from the expression (27) for π i,α above, derived from the Lagrangian, and the definition (29) of X i that we have a primary constraint of the form for each defect i = 1, . . . , N . These N primary constraints are not entirely first or second class. Indeed, their Poisson brackets are non-vanishing on the constraint surface (31) and are not generally invertible. Let {v i p } d i p=1 be a basis of the centraliser g u i ∶= ker(ad u i ) of the element u i ∈ g, with d i ∶= dim g u i for each i = 1, . . . , N . The first class part of each C i is given by the set of constraints C p i ∶= ⟨v i p , C i ⟩ for p = 1, . . . , d i . These satisfy the relations for every i, j = 1, . . . , N and p = 1, . . . , d i , where the last equality uses the fact that v i p ∈ g u i . In particular, we have {C p i , C q j } ≈ 0 for any q = 1, . . . , d j so that the set of constraints C p i for p = 1, . . . , d i , i = 1, . . . , N are indeed first class. It also follows from (30b) that the first class constraints C p i generate right multiplication of the h i by elements e ǫv i p of the centraliser G u i of u i in G -note that under such transformations the g-valued variablesû i are invariant.
Let us extend the basis We claim that the remaining constraintsC r i ∶= ⟨ṽ i r , C i ⟩ for r = 1, . . . , c i contained in C i are second class. We need to show that the matrix {C r i ,C s i } for r, s = 1, . . . , c i is invertible on the contraint surface C i ≈ 0. If this were not the case then we would have ∑ c i s=1 {C r i ,C s i }a s ≈ 0 for some a s ∈ C with s = 1, . . . , c i . On the other hand, we also know from (33) where in the third step we used (32). It follows that ∑ c i s=1 a sṽ i s ∈ g u i which contradicts the assumption that {ṽ i r } c i r=1 is the basis of some complement of g u i in g. We would like to impose suitable gauge fixing conditions D p i ≈ 0, for p = 1, . . . , d i , to fix the first class constraints C p i and move to a Dirac bracket {⋅, ⋅} * which fixes the constraints C i ≈ 0 strongly. In particular, we would like to compute the Dirac bracket {û i1 ,û j2 } * of the g-valued variablesû i = h i u i h −1 i for i = 1, . . . , N . It turns out that the result is independent of the choice of gauge fixing condition D p i ≈ 0. Indeed, consider the variablesX i ∶= h i X i h −1 i . One deduces from (30) that they have the Poisson brackets for each i, j = 1, . . . , N . In particular, it follows from (34b) that {C i1 ,X j2 } = 0 for any i, j = 1, . . . , N . Now the matrix of Poisson brackets between the set of all second class constraints C p i , D p i for p = 1, . . . , d i andC r i for r = 1, . . . , c i is of the block form where the first, second and third block rows and columns correspond to the set of constraints C p i , D p i andC r i , respectively. Each ' * ' denotes a possibly non-zero block matrix. The matrix (35) is invertible since the blocks in position (1, 2), (2, 1) and (3, 3) are all invertible by design. Its inverse is then schematically of the block form Since {C p i ,X j } = {C r i ,X j } = 0 for all p = 1, . . . , d i and r = 1, . . . , c i , the zero block in the middle of the right hand side of (36) implies that the Poisson brackets (34a) will remain unchanged when passing to the Dirac bracket, i.e. we have Finally, using the fact thatX i = −û i after imposing the constraint C i = 0 strongly, we deduce that the g-valued variablesû i for i = 1, . . . , N satisfy N commuting copies of the Kostant-Kirillov bracket To avoid overburdening the notation, and since we shall need to introduce a further Dirac bracket in §3.3, we will denote the Dirac bracket {⋅, ⋅} * introduced above simply as {⋅, ⋅} from now on.

Hamiltonian and secondary constraints
The Hamiltonian density is defined as the Legendre transform of the Lagrangian density (24). However, since the field A t is non-dynamical, i.e. there are no time derivatives of A t in the action, we shall perform the Legendre transform only with respect to the dynamical fields Az, B z and the dynamical variables h i . So we define where in the second line we have used the definition of the bulk constraint C z in (25) and of the defect constraints C i for i = 1, . . . , N in (31). Since we have already set these along with P z strongly to zero, we can drop the corresponding terms in the Hamiltonian density. The Hamiltonian is the integral of the Hamiltonian density over C, namely where in the first term of the right hand side we introduced the notation ⟪X, Y ⟫ ∶= C ⟨X, Y ⟩ dz ∧ dz for any g-valued fields X and Y on C.
Let µ denote the moment map of the Hitchin system [H] (we refer the reader to [BBT,§7.11] for a concise review of Hitchin systems) where in the second equality we have used the constraint C z = 0 in (25) which is now imposed strongly. Introducing also the g-valued field the Hamiltonian can be rewritten succinctly as

Gauge invariance
We need to ensure that the remaining primary constraint, Π t ≈ 0, is preserved under time evolution. That is, giving rise to the secondary constraintμ ≈ 0. We see from the canonical brackets (26) that −μ is the generator of gauge transformations (3) on the fields Az and B z since Note that the moment map µ satisfies the following Poisson bracket where in the second equality we used the relations (41), which also trivially hold withμ replaced by µ. In the third equality we have used the Jacobi identity and the fact that ∂zδ zz ′ + ∂z′δ zz ′ = 0, which follows using the identity ∂z(z − z ′ ) −1 = −2πiδ zz ′ . The Poisson bracket ofμ with itself is therefore where in the second equality we have used (37) for the second term. This vanishes on the constraint surface soμ is first class -we will set it strongly to zero with an appropriate gauge fixing condition in the following section. The time evolution ofμ is given by and therefore we have no tertiary constraints.

Gauge fixing and Lax formalism
We wish to fix the gauge invariance arising from the constraintμ ≈ 0. We will use the gauge fixing condition Az ≈ 0 and simultaneously impose this condition and the constraintμ ≈ 0 strongly by defining a new Dirac bracket. To this end, recall that where the first equality is (41a) and in the last step we have used the new constraint Az ≈ 0. This can be inverted, since Here the subscript (z ′ , 2) means that the pairing ⟨⋅, ⋅⟩ is taken in the second tensor space and the integration is with respect to z ′ . We therefore define the new Dirac bracket, denoted {⋅, ⋅} ⋆ for g-valued functions U and V on C, by .

Lax algebra
We will show that the Dirac bracket of B z with itself satisfies the Lax algebra where r 12 (z, z ′ ) is the standard classical r-matrix To compute this Dirac bracket, we begin by noting that (41b) implies {B z1 (z),μ 2 (z ′ )} = [C 12 , B z1 (z)]δ zz ′ .

Lax matrix
By definition ofμ in (39), it follows that setting this constraint and its gauge fixing condition to zero strongly, i.e.μ = 0 and Az = 0, leads to the equation or ∂zB z = −2πiû i δ zz i in a small neighbourhood of the point z i , which is equivalent to (11). This then leads to the local meromorphic behaviour (12)

Lax equation
At this point we have now fixed all the constraints strongly except for the primary constraint Π t ≈ 0. However, now thatμ = 0 is imposed strongly, the Hamiltonian (40) no longer involves the field A t and simply reduces to In particular, together with the Dirac bracket (42) this now implies the Lax equation (21) in the Hamiltonian formalism We deduce, as claimed at the end of §2, that the time flow ∂ t along the topological direction of the 3-dimensional space R × C is the one induced by the Hamiltonian H P w = P (B z (w)) with respect to the Dirac bracket, i.e. ∂ t f = {H P w , f } ⋆ for any function f of the Lax matrix B z . Focusing on such observables, we are also free to set Π t = 0 strongly since these all Poisson commute with Π t under the Dirac bracket {⋅, ⋅} ⋆ and so their bracket will remain unchanged after introducing a further Dirac bracket to fix the constraint Π t ≈ 0.

Involution
It is well know that the Lax algebra (42) implies the involution property for any pair of G-invariant polynomials P, Q ∶ g → C and distinct points w, z ∈ C. This can also be seen more directly from the above Hamiltonian analysis of 3d mixed BF theory as follows. Since H P w = P (B z (w)) only depends on the field B z we have the involution property with respect to the Poisson bracket (more precisely, recall that {⋅, ⋅} denotes the Dirac bracket introduced in §3.1), for any polynomials P, Q ∶ g → C and distinct points w, z ∈ C. But since H P w is gauge invariant for G-invariant polynomials P and −μ is the generator of gauge transformations, see (41), we have {μ(z), H P w } = 0. The involution property (47), for any polynomials P, Q ∶ g → C, therefore immediately implies the involution property (46), for any G-invariant polynomials P, Q ∶ g → C.

Discussion
In this article we showed that Gaudin models associated with a finite-dimensional semisimple Lie algebra, and more generally tamely ramified Hitchin systems, can be obtained from 3d mixed BF theory in the presence of certain line defects by moving to the Hamiltonian framework and fixing the gauge symmetry using the gauge fixing condition Az ≈ 0.
This 3-dimensional gauge theoretic origin of finite-dimensional Gaudin models is exactly analogous to that of affine Gaudin models in terms of 4-dimensional mixed topological-holomorphic Chern-Simons theory [V2].

Alternative realisations
The Lax matrix of the Gaudin model, or the Higgs field of the Hitchin system, arises from the (1, 0)-form B of the 3d mixed BF theory. In particular, after going to the gauge Az = 0 the latter becomes meromorphic with simple poles (12) at each z i , the location of the type A line defects. The specific choice of line defect (7) led to the residues of B at these simple poles being coadjoint orbitsû i = h i u i h −1 i of some fixed Lie algebra elements u i ∈ g. As is well known, and as we have rederived in the present setting in §3.1.2, such coadjoint orbits provide a realisation of the Kostant-Kirillov Poisson bracket (37).
It would be interesting to see if other realisations of the Kostant-Kirillov Poisson bracket can be obtained by making other choices of type A defects than (7). Indeed, since the field B z satisfies the Lax algebra (42) regardless of the choice of type A line defects we make, the residuesû i at each simple pole z i of B z will necessarily satisfy the Kostant-Kirillov bracket. As mentioned in the affine case in [V2, §4.1], it would be desirable to find the precise dictionary between the possible choices of type A line defects one can introduce in 3d mixed BF theory and the different types of possible representations of the Kostant-Kirillov bracket.

Generalised Gaudin models
We have focused in this paper on the case when the Lax matrix of the Gaudin model, or the Higgs field of the Hitchin system, has simple poles at the marked points z i .
It would be interesting to consider also type A line defects which would give rise to higher order poles in the Lax matrix in order to construct Gaudin models with irregular singularities [FFT,FFRy,VY2]. In the affine setting, generalised surface defects in 4-dimensional Chern-Simons theory leading to affine Gaudin models with irregular singularities were recently considered in [BSV, LV].
Other generalisations of the Gaudin model which one could try to relate to 3d mixed BF theory, or some generalisation thereof, include cyclotomic Gaudin models [Sk,VY1,VY2] or dihedral Gaudin models (see [V1] in the affine case), whose Lax matrices are equivariant under the action of cyclic or dihedral groups, respectively. In the affine case, such a generalisation was considered recently in [Sc2] where the symmetric space λ-model, which can be described as a Z 4 -cyclotomic affine Gaudin model, was obtained along the lines of [V2] starting from 4d Chern-Simons theory with a Z 4 -equivariance condition imposed on the gauge field.
The proposal of [V1], see also [DLMV2,L1], to reformulate non-ultralocal integrable field theories with twist functions as affine Gaudin models similarly provides a deeper origin, more algebraic in nature, of the integrable structure in these theories.
Both the gauge theoretic and algebraic approaches to 2-dimensional integrable field theories, which are of course intimately related [V2], have been used to construct many new examples of 2-dimensional classical integrable field theories in recent years; see for instance [DLMV1,DLMV2,BL,ABL] in the affine Gaudin model setting and the references above in the 4d Chern-Simons theory setting. Finite Gaudin models, or equivalently 3d mixed BF theory, could similarly be used to extend the list of known finite-dimensional integrable systems.
However, the main interest in both approaches lies in their potential to offer new perspectives on various long-standing open problems in quantum integrable field theory, such as the problem of quantisation of non-ultralocal integrable field theories or the search for a deeper understanding of the celebrated ODE/IM correspondence [DT,BLZ1,BLZ2,FFr,LZ]. Indeed, one of the main original motivations in [V1] for reformulating non-ultralocal integrable field theories as affine Gaudin models was the remarkable observation made in [FFr], based on the example of quantum KdV theory, that this may provide an explanation for the ODE/IM correspondence in terms of some suitable affine generalisation of the geometric Langlands correspondence.
By contrast with the affine case, however, the quantisation of the finite Gaudin model, and more generally of the Hitchin system, is extremely well understood; see e.g. [BD,FFR,Fr1,Fr2,MV1,MV2,MV3,MV4,MTV1,FFT,FFRy,Ry]. The connection between 3d mixed BF theory and finite Gaudin models should therefore provide a useful toy model for further developing our understanding of the gauge theoretic approach to integrable models and more generally integrable field theories in the sense of [CY]. In particular, it would be very desirable to understand the Bethe ansatz construction in quantum Gaudin models, and more generally the Gaudin/oper correspondence [Fr1,Fr2,MV1,MV2,FFT,FFRy], from the point of view of quantum 3d mixed BF theory.