Asymptotic four point functions

We initiate the study of four-point functions of large BPS operators at any value of the coupling. We do it by casting it as a sum over exchange of superconformal primaries and computing the structure constants using integrability. Along the way, we incorporate the nested Bethe ansatz structure to the hexagon formalism for the three-point functions and obtain a compact formula for the asymptotic structure constant of a non-BPS operator in a higher rank sector.

Four-point correlation functions are probably the most interesting entities in a conformal field theory. While two-and three-point functions are kinematically constrained by conformal symmetry, four point functions can depend on conformal cross-ratios and will be strikingly different for different conformal theories with different physics.
In principle, the spectrum and operator product expansion (OPE) coefficients of a conformal field theory entail a full non-perturbative solution of a conformal field theory since they can be put together to construct any higher point function. In practice, it is usually unpractical to compute all needed spectra and three point functions and then perform the sum over all possible exchanged operators appearing in the OPE to finally obtain the four point correlator.
In planar N = 4 Super Yang-Mills theory integrability comes to the rescue and renders this task feasible. In this paper, we will construct planar four point functions of large BPS operators at any value of the 't Hooft coupling from the knowledge of two-and threepoint functions which in turn can be computed by means of integrability. We shall be dealing with large enough external operators so that so-called wrapping corrections can be discarded; we denote such four point functions as Asymptotic Four Point Functions.
To compute these four point functions we need to compute the three point functions between two BPS operators and any non-BPS operator appearing in its OPE. These non-BPS operators are described using integrability by a set of (at most) seven different type of Nested Bethe roots [1]. Here we will show that this intimidating Nested Bethe ansatz can actually be described very simply within the hexagon formalism [2] leading to very compact expressions for the relevant three-point functions and hence for the asymptotic four point functions alluded to above.
It would be fascinating to take our final expressions for the four-point correlation functions and initiate a systematic exploration of their various interesting mathematical limits thus extracting various relevant physical regimes, many of which with a relevant holographic interpretation. We look forward to performing these analysis in the near future.
In section 2 we discuss four point functions, their operator product decomposition and the precise limits which allow one to discard finite size corrections. In section 3 we use the Hexagon approach to conjecture all loop expressions for those asymptotic correlators. In section 4 we check the integrability predictions against perturbative data and we conclude in section 5. Various appendices complement the main text.
Note. While this paper was being prepared, we learned of the forthcoming paper [3], which discusses similar subjects (the OPE of four-point functions and its relation to integrability) from a different perspective. We decided to coordinate the submissions to the arXiv.

Super OPE and finite Bethe roots
Defined as

JHEP07(2019)082
the reduced correlator is a nice conformal invariant quantify. It is a function of the SO (2,4) and SO (6) 2 and y ij = y i · y j with y j being the standard six-dimensional null vectors parametrizing the orientation of the external BPS operators which are inserted at four-dimensional positions x j .

Reservoir picture and asymptotic four-point functions
Let us recall some well-understood facts about the correlator (2.1). We will describe it through its infinite OPE series governing what flows from operators 1, 2 to operators 3, 4. In principle, all the multi-trace operators can show up in this OPE representation. However, at large N there is an important simplification: only single-and double-trace operators contribute. Then the four-point function can be expanded as where the conformal blocks F in the first line are fixed by super-conformal symmetry and are summarized in appendix A. Our main focus here is on the last term in the first line corresponding to the contribution of single-trace non-protected operators, whose threepoint functions can be computed by the hexagon approach [2]. A priori, it is non-trivial to disentangle the double-trace contribution from the single-trace contributions since, at finite coupling, they can have the same twist and mix with each other. 1 However, in perturbation theory the twist of each exchanged operator is close to its classical value and this allows us to neatly separate the single-and double-trace contributions -especially if we consider large external operators with p 1 -since the exchanged double traces will have classical twist τ ≥ 2p. Hence, in the OPE limit where z,z are small and the ratio JHEP07(2019)082 + + + · · · + + · · · +  z/z is fixed, for all twists τ ≤ 2p − 2 we can safely restrict our attention to single-trace operators as schematically depicted in figure 1.
As illustrated in that figure, we can think of operators in the OPE, organized by twist, as originating from a big "reservoir" of propagators at the bottom (and top). Operators with a small twist τ flowing in the OPE arise from opening up a few links at the bottom. As such, they will have small side bridges but very large bottom and top bridges. For these operators wrapping in the so-called opposed channel is greatly suppressed in perturbation theory [2,5]. (The adjacent wrapping does matter eventually, at τ /2 + 2 loops to be precise.) In the other extreme case we have the contribution of operators with twist close to the double trace threshold, τ = 2p − O(1). Those have huge side bridges which soak up the reservoirs almost completely. For these large twist operators it is thus the adjacent wrapping which is greatly suppressed. (On the other hand, the bottom and top bridges can now be small so that opposed wrapping eventually kicks in at p − τ /2 + 1 loops.) Finally we have the intermediate regime which is the most relevant one for the present paper. For operators whose twist is very large and yet far from emptying the reservoir, 1 τ 2p -as depicted in the middle of figure 1 -wrapping is suppressed in both the adjacent and the opposed channels. For such contributions we can thus ignore wrapping contributions altogether and use only the so-called asymptotic prediction for the three point coefficients in the OPE expansion.
By playing with the polarization vectors we can easily make sure the adjacent bridges are very large, see e.g. [6]. The basic idea is that if operators O 1 and O 2 have a large nonzero combined R-charge then by R-charge conservation the operators in their OPE must have a large twist, at least as large as the R-charge. For example, we could choose O 1 to be where y 1 = (1, i, β 1 , iβ 1 , 0, 0) , proceed similarly using the remaining complex scalars Y 's andȲ . 2 Combined, the total X U(1) charge of operators 1 and 2 cancels out but the Z U(1) charge does not. Instead there are 2p − 2q units of such R-charge. As such, in the OPE of such operators we have operators whose twist is at least τ = 2(p − q). Those leading twist operators would have side bridges of length l = p − q. Operators with subleading twists will have even larger side bridges. In sum, for the correlator (2. 6) with p and p − q both very large, the side wrapping effects in the OPE channel 12 can be delayed tremendously as they will only kick in at p − q + 2 loops. Furthermore, if q is also very large then the bottom wrapping is also very suppressed since there will be a huge bottom bridge connecting the X's andX's which requires a lot of twist to eat up. More precisely, for a flowing twist τ = 2p − 2q + 2n < 2p − 2 bottom wrapping corrections will only show up at q + 1 − n loops. Only for very subleading twist with n very large will these effect become relevant. To summarize: at weak coupling, for most practical purposes we can ignore all wrapping corrections when computing (2.6). Such four point functions are thus dubbed asymptotic four-point functions.

Super operator product expansion
In the OPE (2.3) we sum over super-conformal primaries only. The descendants are automatically taken into account by the super-conformal blocks F which we summarized out in appendix A.
In the integrability context each single trace operator is described by a set of seven kind of Bethe roots satisfying so called Beisert-Staudacher Bethe equations [1,7]. This description breaks down at some point due to so called wrapping or finite size corrections at which point one must switch to more sophisticated machinery such as the Y-system [8,9], the Thermodynamic Bethe ansatz [10][11][12], finite integral equations [13] or, the current spectrum problem Ferrari, the quantum spectral curve [14]. In this paper we can disregard finite size corrections as explained in the previous section so that the Beisert-Staudacher equations will suffice in what follows.
The notation for the Bethe roots is depicted in figure 2. We use u j to denote the middle node Bethe roots which obey the middle node equations with a spin-chain length L. Then we have a set of three type of Bethe roots v describing one of the su(2|2) wings and another set of roots w describing the other su(2|2) wing. 2 All in all, We use K to denote the number of middle node roots and K (a) andK (a) with a = 1, 2, 3 to indicate the number of Bethe roots in each of the wings. 3 Now, not all solutions to Bethe equations suit our purpose. Super-conformal primaries are solutions to Bethe equations where all Bethe roots are finite. Furthermore, we should exclude solutions where x(v (1) ) = x(w (1) ) = 0 which also correspond to super descendants 4 [1] unless these solutions are part of exact strings in which case the corresponding solutions are denoted as singular solutions and should a priori be considered. 5 The sum in (2.3) stand therefore for a sum over such finite Bethe roots configurations.
The number of Bethe roots of each kind can be read of from the quantum numbers of the exchanged operator. Since our external operators are all BPS, the three-point functions preserve a diagonal su(2|2) subgroup [2,15] which immediately implies that the occupation numbers of the wings must be identified to yield a non-zero result,K a = K a . The relation between the Bethe ansatz occupation numbers and length and the labels

Scaling dimension
Lorentz su(2) × su(2) so(6) R-charge The notation here differs from the one in [1] as {u1,j, u2,j, u3,j, u4,j, u5,j, u6,j, u7,j} there = {v Throughout the paper we will use the sl(2) grading which corresponds to η1 = η2 = −1 in [1]. 4 By means of a dynamic transformation introduced in [1] we can map a root x(v (1) ) = 0 to a root x(v (3) ) = ∞ which corresponds to a descendant of the state without the root at ∞. 5 Coincidentally or not we found out that -on all examples we checked -these singular solutions yield a vanishing three-point function.

JHEP07(2019)082
where the anomalous dimension δ∆ = j 2ig x + (u j ) − 2ig x − (u j ) up to higher loop finite size corrections which, as mentioned above, we are discarding throughout this work. The quantum numbers in (2.8)-(2.11) are all non-negative integers. This puts restrictions over the length L and the occupation numbers K and K (a) . Also, for the Bethe equations to admit finite solutions the occupation numbers usually need to decrease as we go from the middle node occupation K into the wing extremities K (3) , see e.g., [1,16,17].
To summarize, we should a priori find all finite solutions to Bethe equations with K (a) =K (a) , read of their quantum numbers from (2.8)-(2.11), compute their threepoint functions using the hexagon approach [2] and add them up as in (2.3) using the super-conformal blocks summarized in appendix A.
In the next section we will analyze the integrability computation in more detail. We will then observe a remarkable implication of integrability which dramatically simplifies even further the computation of the four-point correlator. It turns out that a version of Yangian symmetry actually implies much more than the global symmetry constraint K (a) =K (a) . To get a non-vanishing OPE contribution we must in fact have absolutely symmetrical wings root by root, that is v j . This is a very sharp and novel spacetime implication of the world-sheet integrability.

Nested Bethe wave function
The crux of the hexagon formalism lies in cutting a pair of pants, which represents the structure constant, into two hexagonal patches. Upon cutting, magnons in each operator are divided between two hexagons, and we sum over all such possibilities with appropriate weights, namely a propagation phase e ip and S-matrices. When magnons belong to a rank 1 sector, things are rather simple since the S-matrix is just a scalar phase. In general however, one has to deal with a complicated index structure: to understand this point, let us consider two-particle states with indices A and B. After being cut, it produces a complicated set of states as shown in figure 3. Each term in the figure can be computed using the hexagon proposal in [2], but the result after summing up different terms (and also the indices C and D) in general takes a complicated form. In addition to this difficulty, there are two extra complications: first, when we have several magnons with different indices, the contribution from each hexagon itself has a matrix part which can be quite complicated. Second, if we want to compute a structure constant, we need to require that the external operators are eigenstates of dilatation operators. To do so, we need to contract the indices in figure 3 against a specific wave function 6 for the indices Ψ AB so that the 6 Note that the "wave functions" that we discuss in this paper are all about the wave functions for the indices; in the nested Bethe ansatz language they are the wave functions at the nested level. This is in contrast to the perturbative computation at weak coupling in which one also needs to talk about the wave function of magnons themselves (in other words, the wave function for the middle node). The structure of the wave function of magnons themselves is already taken into account by the hexagon form factors, and we simply need to figure out how to efficiently deal with the index structure. This is the main theme of this paper. Figure 3. Splitting a two-particle state with indices. When the first particle passes through the second particle, it gets multiplied by a S-matrix S(u 1 , u 2 ). The resulting state is a complicated object which includes a summation over indices (C and D in the figure).

JHEP07(2019)082
resulting state is an eigenstate. The contraction with the wave function "couples" the contribution from two hexagons and makes the computation even more complicated.
The way to circumvent such complication of indices is, as is well-known, the Nested Bethe Ansatz. In the Nested Bethe Ansatz, we first make an ansatz for the wave function of the flavour indices, which depends on the order of momentum-carrying roots u i , and a set of roots at higher levels w. It has an important property that it "diagonalizes" the action of the S-matrix, with S(u, v) being an abelian phase. When w satisfy the Bethe equation for higher levels, the wave function has an additional "nested periodicity" property, where f is a theory-dependent phase factor. With these two properties, one can rewrite the right hand side of the periodicity condition of the full wave function, in the following way: This leads to the Bethe equation for momentum-carrying roots . . , u 1 Figure 4. Matrix part forᾱ = ∅: one can simply act the S-matrix to the right wave function ψ R and simplify the structure. Note we are using the normalization of the matrix part, in which the abelian part (sl(2) S-matrix) is unity.
Alternatively, we can use this relation to read off the phase factor f from the Bethe equation.

Nested hexagon
We now use the aforementioned properties to compute general asymptotic three-point functions for two BPS and one non-BPS operators. Bethe states in N = 4 SYM spin chain are characterized by seven sets of roots, which are split into two "wings" by the momentum-carrying roots u (see figure 2). Owing to this structure, the wave function at the nested level is given by a product of two wave functions Ψ (left) andΨ (right), which are identical if we relabeled the undotted by dotted indexes and trade v by w.
To apply the hexagon formalism, we first reorder the magnons (or equivalently the momentum-carrying roots) and split them into two subsets α andᾱ. Thanks to the property of the nested Bethe wave function, the state we get after reordering is as simple as the original one: Here S(u, v) is the S-matrix in the sl(2) sector. The next step is to contract the nested wave function with the hexagon form factor. Whenᾱ is empty, this can be done easily since the hexagon is essentially given by a product of S-matrices which are already diagonalized by the wave function. As a result, we obtain where ← − Ψ ← − u |Ψ u represents a contraction of rightΨ and left Ψ wings defined pictorially in figure 4. More precisely, for the following two states in the psu(2|2) spin chain with the JHEP07(2019)082 f (u j , w) Figure 5. Matrix part forᾱ = ∅: the magnons for the right wing are in a different order from those in the wave function. To contract the wave function with the hexagon, we have to rewrite it using the "nested periodicity" (3.2).
opposite orderings of the inhomogeneities, with h χχ being the one-particle hexagon form factor defined in [2]. By contrast, whenᾱ is not an empty set, things are a little bit more involved. In the diagram that computes the matrix part, the magnons for the right wing are in the orderᾱα (see how the lines enter into the leftmost blue dotted box in figure 5) while the wave functionΨ (in the same figure) is originally defined in the order αᾱ. To simplify the computation, we rewrite the wave functionΨ in the orderᾱα (rightmost dotted box in figure 5). This can be done by using the nested periodicity (3.2), and it produces a product of phase factors f , which can be read off from the asymptotic Bethe equation for the momentum-carrying root: . (3.10) Here x(u) is the Zhukowski variable u = g(x + 1/x) with g = √ λ/4π, and f ± (u) ≡ f (u ± i/2). Note that we should only take factors which depend on v's since f is the phase JHEP07(2019)082 factor coming just from the right wing. After doing so, we can straightforwardly contract the wave functions with the hexagons and act the S-matrices on the wave functions. This leads to an expression In order to compute the contraction ← − Ψ ← − αᾱ |Ψ αᾱ we need the explicit form of the nested su(2|2) wave-functions which can be found in [65,66]. By performing this computation in the psu(1, 1|2) subsector with corresponding wings in su(1|1), see appendix F, we learn that on the support of the Bethe equations (on-shell condition): where the right hand side is the usual scalar product of the spin chain and is independent of the ordering of the momentum-carrying roots. Hence it can be taken out of the sum over partitions as a partition-independent prefactor. We did not perform the same computation in the full sector but conjecture that equality (3.12) still holds when the full su(2|2) wings are excited. 7 Furthermore, because of the orthogonality of two different on-shell states, the scalar product (3.12) vanishes unless all the roots of left and right wings are identical, namely v (i) = w (i) . This suggests the existence of a hidden symmetry, which forces infinitely many structure constants to vanish. In the next subsection, we explicitly construct such a symmetry using the transfer matrix of psu(2|2).
Putting together all the elements, we obtain our main formula for the structure constant in higher rank sectors, Here µ(u) is the measure [2] and A is a higher-rank generalization of the sum over partitions, which reads . (3.14) The factor u|u denotes the Gaudin norm of the full psu(2, 2|4) spin chain whereas v|v , equivalent to (3.12), is the norm for the (left) wing. Their precise definitions are given in JHEP07(2019)082 Figure 6. Definitions of the Gaudin norms u|u and v|v : φ u,v,w is a logarithm of the nested Bethe equation, e iφu,v,w = 1. u|u is a determinant of the full matrix shown above whereas v|v is a determinant of the upper (or equivalently lower) diagonal matrix shown in red (green). The missing matrix elements are vanishing, as a result of the structure of the psu(2, 2|4) BAEs (e.g., there is no interaction between auxiliary roots w and v lying on different wings). figure 6. In the sl(2) sector there are not wings so we have v|v = 1 and f (u) = 1 such that (3.13) reduces to equation (31) in [2]. 8 As shown in appendix B, the ratio v|v 2 / u|u can be rewritten as a single determinant by eliminating the dependence on the roots at higher levels. It is also possible to express the sum over partition A as a Pfaffian of a 2K × 2K matrix. Combining these two expressions, we can express the square of the structure constant simply as a ratio of two determinants. See appendix D for details.
In section 4, we will use this formula to reproduce the data obtained by the OPE decomposition of the four-point functions.

"Yangian" symmetry
As we saw above, the structure constant vanishes unless the roots in the two wings are identical. Below we uncover the underlying symmetry responsible for such a super-selection rule.
Let us take a look again at the matrix part of C ••• 123 . Using the Yang-Baxter relation, we can show that the difference of transfer matrices acting on two wings must always vanish JHEP07(2019)082 if the state is contracted with the hexagon (see figure 7): Here r can be any representation of psu(2|2) and T and ← − T denote the forward and the backward transfer matrices acting on the left and the right psu(2|2) respectively. This property turns out to be true even more generally: it is in fact easy to see that this also holds for correlators with more than one non-BPS operators, and it even holds in the presence of wrapping corrections if we ignore the subtleties coming from double-pole singularities [18].
As is well-known, the expansion of transfer matrices yields mutually commuting charges. Thus, the relation (3.15) is manifestation of infinitely many conservation laws hidden in the three-point function. With a slight abuse of the word, we call it "Yangian symmetry" in this paper.
When the state we contract is the on-shell nested Bethe state, we can replace the symmetry generator T r (u) − ← − T r (u) by its eigenvalue. Then it follows from (3.15) that, unless is satisfied as an eigenvalue equation, the structure constant must vanish. The eigenvalues of these transfer matrices are expressed in terms of nested roots [19] and the only possible way to satisfy (3.16) is to set the rapidities of two wings to be identical. This is the symmetry origin 9 of our super-selection rule.
In integrable systems, an infinite number of commuting charges are often accompanied by a real Yangian symmetry, namely a set of non-commutative non-local charges. An explicit construction of such a symmetry for our case is an interesting open problem for the future.

Comparison with data
In this section we combine the previous two sections. Namely, we put the integrability predictions of section 3 to test by comparing them with the OPE expansion described in section 2 of perturbative four point functions.

JHEP07(2019)082
Integrability yields predictions for individual structure constants given a set of Bethe roots corresponding to the non-BPS operator at hand. These operators have different quantum dimensions ∆ as read off from their Bethe roots but classically there is a large number of operators with the same classical dimension ∆ classical = ∆ − δ∆ -the right hand side of (2.8). In perturbation theory, what shows up in the OPE of a four-point function are sum rules over these degenerate operator spaces. The summand is the square of the structure constants weighted by powers of the quantum anomalous dimension δ∆.
This was illustrated in detail for the simplest sl(2) sector in [6], see for example formulae (56-65) therein summarizing some of those sum rules. Here we are dealing with the full nested space so the sum rules are a bit more involved; they are sums over all finite solutions to Bethe equations whose occupation numbers K, K j and length L yield the same psu(2, 2|4) classical charges appearing in the right hand side of (2.8)-(2.11).
Up to one loop, for instance, we can easily use the perturbative results of [22] to extract predictions for the sum rules P (0,0) , P (0,1) and P (1,1) defined as 10 Bethe solutions with fixed r.h.s. of (2.8)-(2.11) These predictions are summarized in table 1 for one loop OPE data extracted from the four point function of 1/2 BPS operators of length p = 4. We provide the sum rules corresponding to non-BPS operators with so(6) charges 0 ≤ m ≤ n ≤ 2 and twist τ ≤ 2 p−2. At twist τ = 2p the OPE would be contaminated by double trace contributions and we could no longer cleanly test the single trace integrability predictions against it. 11 To generate this table we used a four point function with external operators of length p = 4 so that at one loop we can test integrability predictions in the OPE involving operators of twist 2, 4 and 6. Had we used a larger p and we could have tested those twists (the result would be the same at this loop order) and more.
Having predictions for the right hand side of (4.1) we turn to the left hand which we will now generate using integrability. 10 In other words, 11 Of course, we could simply consider larger external operators to delay the double trace contribution as much as we want. It would also be interesting to play with the OPE analysis varying the external dimensions in order to isolate the double trace contribution from the extremal one. This would provide valuable data for guiding an integrability based approach towards studying extremal or double trace correlation functions.
and spin s. The data in color has been checked against integrability and the data in red is for the sl (2) sector.

JHEP07(2019)082
To reproduce this OPE data from integrability we first find all (wing-symmetric) solutions 12 of Beisert-Staudacher Bethe equations with finite Bethe roots. First we set g = 0 and solve these equations to leading order at weak coupling. This part is hard. Then to get the quantum corrections to the Bethe roots we simply correct the Bethe roots perturbatively by linearizing the O(g 2 ) Bethe equations around each tree level seed solution. This part is straightforward. Once we get the Bethe roots we plug them into the Hexagon prediction (3.13) with l = L/2, sum over all solutions as indicated in the left hand side of (4.1). The result is then compared with the OPE predictions for the right hand side of (4.1) which we extracted from perturbative data and summarized in table 1.
To find all the Bethe ansatz solutions at tree level we resorted to various pieces of technology. The simplest Bethe equations correspond to the sl(2) sector where we only excite the middle node and Bethe solutions for operators of spin s are given by sets of real roots {u 1 , · · · , u s }. Solving Bethe equations in the sl(2) sector in Mathematica is absolutely straightforward, see for example [23]. Checks of OPE data against integrability conjectures were already performed in [2] and earlier in [6]. The sums P (a,b) for this sector are highlight in red in table 1.
For other sectors such as higher rank sectors with excited wing nodes, Bethe equations become more complicated and also admit complex solutions including at times so called exceptional solutions [24,25]. It is the existence of complex solutions which renders the problem of finding all solutions to the Bethe equations much more challenging in this case. One way to proceed which we found quite useful is to use a Baxter formulation of Bethe equations and solve directly for the Baxter polynomials and the transfer matrix eigenvalues rather than individual Bethe roots. Another useful numerical method is the socalled Homotopy continuation method [26] where one starts from some simpler equations and adiabatically deform them until they become the Beisert-Staudacher equations. For su (2) solutions this was proposed in [26] and its generalization to the nested case also works very well. The third method -and the one we found to be the most convenient -is however to use the very powerful recently proposed analytic solver of [27,28] based on the Q-system. 13 This provided us with the complete set of Bethe solutions needed to reproduce all the number in blue in tables 1 using the hexagon conjecture (3.13).
To illustrate what goes into these computations consider the following example.
− g 2 1.81175390266299357208448117539026629935720844 12 Bethe solutions with asymmetric wings give a vanishing structure constant C ••• . We exclude as well symmetric solutions with w 1 = v (1) = 0, as they do not render highest weights. 13 We are very grateful to C. Marboe and D. Volin for sharing a working code of the fast analytic solver for psu(2, 2|4).

JHEP07(2019)082
we then recognize this renders the rational numbers: In an attached Mathematica notebook the reader can find our conjecture (3.13) coded up to one loop order and how the twenty solutions beautifully add to this nice rational number which perfectly reproduces the perturbative OPE data. All other blue numbers in table 1 were checked in the same way.
Note in particular that there is no data in table 1 when s + n − m is odd although there exist Bethe solutions with these global charges. Their absence is imposed by the symmetry of the OPE under the exchange of the pair of identical external operators. It is nice to see how this selection rule comes about from our integrability construction. The sum over partitions (3.14) is written in terms of the bridge length 31 . However, nothing in the original problem singles out this particular adjacent channel. We can find an equivalent formula expressed in terms of the complementary bridge length 23 when the Bethe state is on shell and cyclic. Namely, using the ABA equations, e ipᾱL 3 f (uᾱ) 2 Sᾱ α = 1 and f (u α )f (uᾱ) = 1, the permutation property of the hexagon form factor h αᾱ Sᾱ α = hᾱ α , and the zero momentum condition e −ipᾱ = e ipα , one easily derives that where K = |α| + |ᾱ| is the total number of magnons. For two identical operators, the spin chain is split into two equal parts of length 13 = 23 = L 3 /2 and the previous relation turns into a selection rule: A = 0 for K odd. In terms of the quantum numbers of the superconformal primary, see equation ( Finally, it is worth stressing that while all the checks we performed worked like a charm, they do not exhaust the available perturbative data by any stretch. Even at tree level and one loop we only confirmed the predictions in blue in table 1. From a Bethe ansatz point of view, most of these solutions are not general enough as they do not excite roots associated with all psu(2, 2|4) Dynkin nodes. The only solutions which contain roots of type v (3) and w (3) are the ones which contribute to P n=0,m=0 τ =6,s=0 in table 1. These solutions have some peculiarities, such as the appearance of odd powers of g in the rapidities, which are explained in appendix C. It would be very interesting -even at this low loop orderto perform a few higher twist checks and probe various Bethe solutions in full generality. It would also be very interesting to expand Bethe ansatz further and compare the integrability predictions with the available data at two [29][30][31][32], three [33] or even four loops [34][35][36][37][38]. When going to higher loops we should either start including finite size corrections to the three-point functions [5,18,39] (hard) or increase the length p of the external operators as explained in section 2 (easy).

JHEP07(2019)082 5 Conclusions and outlook
We initiated the study of four-point functions of large BPS operators by combining the operator product expansion and integrability. We found a compact formula for structure constants in higher rank sectors and checked it against OPE data at tree level and 1 loop.
There are numerous physically interesting questions which can be addressed with the methods and the results in this paper. The most important among them 14 is, perhaps, to study the strong coupling limit and understand the emergence of the bulk locality. Since we are discarding the double-trace contributions, we need to take a careful double-scaling limit in which both the coupling g and the lengths of the external operators p i tend to infinity. In such a limit, the solutions to the Bethe equation become dense and the summation over the solutions would be effectively replaced by integrals. It would be extremely interesting to see if such integrals, combined with the Pfaffian/determinant representations presented in this paper, give us analytic control over the local physics in AdS.
Right now, we are solving the Bethe equation and computing individual structure constants. However, what we really want to know is the full four-point function, not the individual structure constants. Furthermore, solving the Bethe equation would become horrendously complicated when the length of the operators and the number of magnons are large. It is thus important to develop a method which allows us to compute the sum over the states without explicitly solving the Bethe equations. 15 Quite recently, there appeared an alternative approach to study higher-point functions called hexagonalization 16 [43]. In that approach, the hexagons are glued always along the mirror edges and the cross ratios appear as the weight factors for mirror particles. By contrast, in our approach, we glue physical edges and the cross-ratio dependence comes entirely from the conformal blocks.
Both approaches have its own advantages and drawbacks. For example, for large operators and for the leading OPE contributions, the approach proposed here extends to higher loops in a trivial way since after solving Bethe equations at tree level it is trivial to correct them perturbatively to any loop order. The hexagonalization approach, on the other hand leads directly to beautiful resummed expressions for the full four point function without ever solving any Bethe equations but the number of mirror particle integrals grows with the loop order. The approach here struggles when we reach maximal twist and start becoming contaminated by double trace operators while the hexagonalization method automatically incorporates these effects. Super-conformal symmetry is manifest here since 14 Other interesting regimes to study would be the Regge limit and the double light-cone limit. 15 In the context of the scattering amplitudes, there are representations of amplitudes given by summation over algebraic equations known as scattering equations [40] which are formally not so dissimilar from the Bethe equations encountered here. There, beautiful methods have been developed to perform to sum over the solutions to the scattering equation without even explicitly solving those equations, see e.g. [41,42]. Can we do something similar here? 16 For tree-level four-point functions involving non-BPS operators, an extensive study in the SU(2) sector was performed in [44], and expressions in terms of sums over partitions, akin to the ones that arise from the hexagon formalism, were derived. A similar correlator in special kinematics was revisited recently [45] and it was observed that the result can be expressed in terms of hexagon form factors, based on the idea akin to hexagonalization.

JHEP07(2019)082
we sum over super-conformal primaries with the help of super-conformal blocks but crossing symmetry is far from obvious while the converse is true in the hexagonalization approach. The list could go on and on. Clearly, understanding the relation between two approaches is a very interesting future problem. If we could take the best out of each of them we would call it a win!

Acknowledgments
We are grateful to Christian Marboe and Dmytro Volin for discussions and sharing a useful Mathematica code with us. We thank Carlo Meneghelli for comments on the draft. Research at the Perimeter Institute is supported in part by the Government of Canada through NSERC and by the Province of Ontario through MRI.

A Superconformal blocks
The super-OPE relies on the use of a superconformal block to resum the contributions of all the superconformal descendants of a given superconformal primary with weights ∆, s, m, n. In this appendix we record its expression when the superconformal primary that is flowing is either long or half-BPS and when the external operators are identical half-BPS superconformal primaries. The non-BPS block presented in (A.1) can be read off from various references, e.g. [46][47][48][49]. In contrast, less is written explicitly about the BPS block. Following [46,47] it has become conventional to decompose the protected part of the four point function over an OPE-like basis of (single variable hypergeometric) functions [46,48], which solve a SUSY version [49] of the Casimir eigenvalue equation [50]. These functions are however not identical to the BPS blocks we are after, as one can easily see by checking their content in usual conformal waves. We give in (A.5) the expression we have found for the BPS block by adding enough of the former functions together until we got the appropriate OPE content for an half-BPS multiplet.
Non-BPS blocks. These are given concisely by where P n are Legendre polynomials and where h λ (z) = z λ 2 F 1 (λ+2, λ+2; 2λ+4; z), with the shifts in red being an N = 4 SUSY shift. In AdS/CFT jargon, we can say that the first line is SUSY, the second is AdS and the third accounts for the sphere. The slightly unconventional normalization factor in the last line ensures that the R-charge blocks behave as Y n,m (α,ᾱ) = 1 × α −n−2ᾱ−m−2 (1 + O(α,ᾱ)). The bosonic blocks are JHEP07(2019)082 By sending z,z → 0, one reads out the su(4) block of the superconformal primary, with Dynkin labels [n − m, 2m, n − m], while by sending α,ᾱ → 0 one recovers the conformal block of a SUSY descendent with dimension ∆ + 4 and spin s, 17 Equivalently, the so(2, 4) block G ∆,s for a conformal primary with dimension ∆ and spin s is given [52] by F ∆,s without the shifts in red in the arguments of the hypergeometric function h below (A.1).
BPS blocks. In order to find the 1 2 -BPS block F ∆ we can start with an Ansatz given by a linear combination of the solutions of the super-Casimir equation. Schematically this is which includes the short solution G short ∆ defined in [49] 18 and the long blocks of the form (A.1) with a vanishing super-conformal Casimir eingenvalue. This is such that F ∆ satisfies the super-conformal Casimir equation with the same eigenvalue as G short . We fix the coefficients c ∆,k of the zero modes by demanding the correct OPE behaviour F ∆ = (zz) ∆ Z ∆,∆ (α,ᾱ)(1 + (z,z)).
Our final result can be written as linear combinations of six bosonic blocks, The formula (A.5) agrees with those given in [51] for ∆ = 2 (irrep 20 = [0, 2, 0]) and ∆ = 4 (irrep 105 = [0, 4, 0]), see equations (8.17) and (8.24) in [51]. Each of the 6 conformal waves in (A.5) corresponds to one bosonic conformal primary, with zero hypercharge Y = 0 and in a left-right symmetric irrep of so(3, 1) and su(4), on the middle line of the half-BPS supermultiplet given in table (B.1) of [51]. For the stress tensor multiplet ∆ = 2, only the first three terms survive in (A.5), corresponding to the dimension 2 chiral primary, the JHEP07(2019)082 dimension 3 R-symmetry current and the dimension 4 stress energy tensor, in agreement with the bosonic (hypercharge zero) components in the ∆ = 2 supermultiplet reviewed in table (2.15) of [51]. Notice that the dual and self-dual parts of the dilaton carry non zero hypercharges and thus decouple, in accord with the non-renormalization property of the BPS structure constant.

B Dualization, Gaudin norm and diagonal symmetry
The considerations in section 3 were based on the description of the Bethe states in the sl(2) grading, corresponding to η 1 = η 2 = −1 in the notations of [1]. This choice was motivated by the comparison with data. Conformal primaries in the sl(2) grading are indeed closer to the superconformal primaries (they share the same length in the spin chain description for instance). We would nonetheless get a fully equivalent description in the su(2) grading, corresponding to the choice η 1 = η 2 = +1, as shown in appendix F in a particular subsector. As well known, the two gradings are related by the (simultaneous) dualizations of the fermionic nodes in the two wings of the psu(2, 2|4) Dynkin diagram [1]. In this appendix, we will show that our main formula (3.13) transforms properly under this diagonal dualization. We shall also demonstrate that it is invariant under diagonal su(2|2) D transformations, including the length changing effect accompanying them [1]. To prove both properties, we shall find convenient to first show that the ratio of determinants appearing in (3.13) can be written as a single Gaudin determinant for the effective BAEs for the main roots u. The latter are obtained after implicitly integrating out the auxiliary rapidities along the wings, at given wing mode numbers m v = m w .
Cosmetic rewriting. To simplify the discussion, we assume that we can unite the fermionic roots of type 1 and type 3, on each wing of the diagram 2, by applying the dynamical transformation of [1]. It amounts to inverting the Zhukowski roots of type 3 such that they appear as roots of type 1, while simultaneously redefining the length of the operator, L → L − 2K (3) , in order to absorb the momentum factor spit out during the inversion. It is clear from the structure of the wing dependent factor (3.10) why we can do that in the sum over partitions (3.14). It will become clear, at the end of the appendix, why we can also do it at the level of the Gaudin determinants (once these ones are properly projected down to the subspace of cyclic states).
The formula for the higher rank structure constant (3.13) factorizes into two main blocks. The first one is the ratio of determinants r, with v = Ψ and w =Ψ the on-shell left and right wing wave functions and with G = u|u the Gaudin determinant of the full Bethe state, see figure 6. The next one is the partition dependent factor

JHEP07(2019)082
with implicit products over the elements of the various sets and with T = Ψ| T |Ψ / Ψ|Ψ the eigenvalue of the SU(2|2) transfer matrix in the diagonal state Ψ =Ψ. The latter transfer matrix is evaluated in (B.2) on the Bethe rootsᾱ ⊂ u, We normalize it such that T (u i ) = 1 for u = {u i } in the sl(2) sector. For a general Bethe state, one finds that T (u i ) = f (u i ) with the f factor (3.10), see e.g. [19]. There is no need to carry out the sum over the partitions, since our discussion will apply to each partition independently. We also ignore the remaining overall factors in (3.13) and (3.14), like the product over the measures, etc., which also play no role at all. The transfer matrix (B.3) responds promptly to all the manipulations we want to perform. Changing "its" grading, for instance, is straightforward, using the general expression given in [19]: it boils down to replacing the auxiliary roots y and their S-matrices S by their dual versionsỹ andS. (This is so up to an overall factor Aᾱ α that is used to convert h αᾱ in (B.2) from its sl(2) to its su(2) value [2].) What is less obvious is that the ratio (B.1) responds in the exact same way. The problem being that in (B.1) one needs to take derivatives of the BAEs, before dualizing. As we shall see, for the ratio (B.1) one can equivalently proceeds in the reverse order, that is dualize before taking the derivatives.
Induced Gaudin determinant. The separation between main roots u and auxiliary roots v is suggestive of a factorization into two determinants for the two subsystems of equations φ u = 2πm u and φ v = 2πm v . Were these two subsystems independent, we would of course immediately conclude that (B.4) (Note that we will not need to distinguish between the two types of auxiliary roots that we have at our disposal. This is why we unite them into a single set, v ∪ w → v. The cut off between real and auxiliary roots is actually immaterial and our discussion applies to a general decomposition into two or more non-overlapping subsystems.) The factorization (B.4) would also apply to triangular systems, that are such that the dynamics of one subsystem does not depend on the complementary subset of roots. If, for instance, the u system of equations φ u = 2πm u is such that ∂ v φ u = 0, then we would still have (B.4) except that G v → G v|u , with G v|u = det ∂φ v i /∂v j the determinant of the v system of equations φ v = 2πm v (with the roots u entering as external parameters). The point is that we can always bring the full system to a triangular form if we treat the v's as being the slaves of the u's. What we get in the more general case is that the factor G u is replaced by the determinant G u|φv = det dφ u i /u j , with the derivative d/du j being taken at fixed mode numbers φ v = 2πm v instead of fixed rapidities v.
The proof of this factorization is elementary. We simply need to recall the interpretation of the determinant G as the Jacobian for the mapping between rapidities and mode numbers,

JHEP07(2019)082
and evaluate it in two steps, We would arrive at the same conclusion starting with integrating over v, at fixed u, and then integrating over u, As alluded to before, G v|u is the minor of G, obtained by deleting the φ u -rows and ucolumns, while G u|φv is the induced Gaudin determinant for the set of effective equations for the u's. The latter are obtained by 1) solving the equations for v at given mode numbers m v and for a given set of rapidities u (assuming that solutions exist) and 2) plugging the solution v = φ −1 v (m v , u) into the equations for u (assuming that it is unique). The nice thing about G u|φv is that it eliminates the dependence on the variables used to parameterize the directions transverse to the subspace φ v = 2πm v . It indeed measures the density of states (per volume du) on a given "physical subspace" φ v = 2πm v .
In the case of interest we have a tripartite decomposition u ∪ v ∪ w where v and w do not interact with each other, see figure 6. Hence we can write and since v = w, as a result of the on-shell super selection rules, we also have G v|u = G w|u . Therefore 1 It shows that 1/r 2 is the induced Gaudin determinant for the roots u on the subspace φ v = φ w = 2πm v . As such, it should be clear that it does not depend on how we parameterize the higher levels of the wave function and, in particular, on which grading we choose; it is invariant under dualization of the auxiliary roots v →ṽ, since, on both sides of this equation, the total derivatives are taken along the same physical subspace. Put differently, when computing G u|φv,w we are allowed to dualize (or equivalently use the on-shell conditions for the v ∪ w's) prior to take the derivatives w.r.t. the u's.

JHEP07(2019)082
Diagonal symmetry. Let us comment finally on the diagonal su(2|2) D symmetry. This is a symmetry of the structure constant [2]. It should thus be reflected in our final expression (3.13). As well known, see e.g. [19], this symmetry can be phrased as the invariance under addition of auxiliary roots at the special points y = 0, ∞. 19 For y = 0 the symmetry is dynamical and related to the joining/splitting of free multiplets at the unitarity bound [1]. (This is literally true for states in the su(1, 1|2) sector.) It holds for cyclic states only and comes along with a redefinition of the length of the spin chain [1]. 20 Roughly speaking, adding or removing a root at y =ẏ = 0 soaks up or spits out two units of spin chain length. For example, a length L two derivative BMN operator and a length L + 2 two scalar BMN operator, sharing the same roots u 1 = −u 2 , fall in the same supermultiplet and differ by the adjunction of a root at y =ẏ = 0 (+ some at infinity). The diagonal symmetry is manifest in (B.2) since the transfer matrix is invariant under sending roots to infinity. It also transforms multiplicatively, in the spin chain frame, under addition of a root at y = 0, hence implementing the length changing effect L → L − 2 for the removal of the fermionic roots y = 0 andẏ = 0 from the state. (We work here in the non compact grading.) It is also not difficult to prove that the ratio (B.1) remains unchanged when sending auxiliary roots to infinity. It does not transform correctly when setting roots at y =ẏ = 0 however. The problem is that the Gaudin determinant requires to take derivatives of the BAEs prior to impose cyclicity, while we would need the reverse to make our point. The latter two operations do not quite commute, which is why the ratio (B.1) is not per se a diagonal invariant. The mismatch is not that big however and, as shown below, drops out of the full structure constant. Given our earlier discussion, it should be clear that the quantity that is invariant is the induced norm on the subspace of cyclic states. This is not quite the same as (B.1). The Gaudin determinant G u|φv,w | v=w measures the density of unconstrained spin chain states on the subspace φ v = φ w = 2πm v . The latter counts L too many states, as compared to the gauge theory, because of the L−1 subspaces of unrealized solutions to the equation e ipuL = 1 (which itself is a consequence of the BAEs). Hence, imposing the cyclic condition e ipu = 1 at the level of the determinant amounts to rescaling it by L, and the invariant quantity is and not r alone. This small extra factor makes a difference since the length L transforms in the multiplet splitting/joining. The product Lr 2 should not. In the hexagon construction, the factor √ L multiplying r is provided by the vacuum structure constant 19 One must also include w = ∞, with w = v (2) an su(2) root, as well asw = ∞ in the dual frame for the second su(2) ⊂ su(2|2). Combinations of the type (y, w) = (0, ∞) and (y, w) = (∞, ∞), with w − g(y + 1/y) held fixed, should also be considered to get all the supercharges.
where N is the rank of the gauge group. The (properly normalized) structure constant is thus an su(2|2) D invariant, as expected.
Technically, to prove the invariance of Lr 2 , we write where L/(L−2) emerges as the ratio of two Jacobians: dLP/dP = L, obtained by replacing the equation LP = 2πm by the cyclic constraint P = 2πm in the determinant for the initial spin chain of length L, 21 and dP/d(L−2)P = 1/(L−2), obtained by reversing the procedure for the determinant of the final spin chain of length L − 2. Inbetween we have applied and similarly forẏ, which are valid as soon as the derivative d/du k is taken along the cyclic subspace y =ẏ = 0. By the same token, one can demonstrate that the dynamic transformation mentioned at the very beginning is correctly implemented in our expressions.

C Comparison with data: a special case
In this appendix we analyze in detail how to obtain the OPE data sum P n=0,m=0 τ =6,s=0 in table 1 using the conjecture (3.13). Unlike the rest of the data that we have checked, this sum receives contributions from operators of different lengths. For these quantum numbers we have a total of 7 wing-symmetric Bethe solutions: Lenght L Field content at O(g 0 ) # Roots in sl(2)-grading # Wing-symmetric sols 4 T r(DDZZZZ) + · · · {1, 2, 3, 6, 3, 2, 1} 2 6 T r(ZZZZZZ) + · · · {0, 2, 4, 6, 4, 2, 0} 5 (C.1) The five solutions with L = 6 correspond to operators in the so(6) sector at O(g 0 ). At loop level, their roots receive corrections in even powers of the coupling g and can be used straightforwardly in (3.13) to obtain the corresponding structure constants. Theses solutions behave in a standard way so we do not review them in this appendix.
The two solutions with L = 4 have roots in all the 7 nodes of the psu(2, 2|4) Dynkin diagram. Hence it constitutes an interesting case that proves all the components of the conjecture (3.13). In the following we analyze in more detail these solutions.

JHEP07(2019)082
At g = 0 these solutions contain two vanishing fermionic roots v . These however are not associated to the action of supercharges. As we show in table 2 these zeros receive corrections at loop order and are lifted to take opposite non-zero values. Their corrections start at O(g 1 ), unlike the rest of roots that start at O(g 2 ), and have an unusual expansion in odd powers 22 of g. In terms of the Zhukovski variables the relation between the fermionic roots translates into: At loop level we can perform a dynamical transformation of the fermionic roots, as explained in appendix B. We can treat the roots of type (3) as of type (1), in both wings, by going through the cut (x → 1 x ) and increasing the length of the operator. In this way, at loop order, the operator with length L = 4 and excitations {1, 2, 3, 6, 3, 2, 1} has an alternative description with length L = 6 and excitations {0, 2, 4, 6, 4, 2, 0}. In this latter description the new fourth wing root v Zhukovsky: x(v

D Pfaffian representations
In this appendix, we show that the sum over partition (3.14) at finite coupling can be recast as a Pfaffian of a finite-dimensional matrix. This rewriting has two virtues: it reduces the cost of numerical computation, which is extremely heavy when the number of roots is large. In addition, the argument is potentially applicable to the mirror corrections as well.
The first step is to rewrite (3.14) as The next step is to use the identity, 24 which holds for even N , where pf denotes a Pfaffian of a matrix and k(x, y) is given by Using this identity, we can rewrite a product of H(u, v) over a set of rapidities s = {u 1 , u 2 , . . .} as We thus have The series representation (D.8) is akin to the expansion of the so-called Fredholm Pfaffian [54]. In fact, using the expansion formula for the Pfaffian, we can recast it into a Pfaffian of a 2K × 2K matrix as follows: Here J and E are given by .
(D.2) 24 We found this formula empirically using Mathematica. We then learned that it is a special case of the elliptic generalization of Schur's Pfaffian formula (see eq. (16) of [53]).

JHEP07(2019)082
with I M being the identity matrix of rank K. Using (D.9) and the well-known fact that a square of a Pfaffian is a determinant, we can express A 2 in (3.13) as a determinant. After rewriting a bit, the result reads where K has a block structure K ≡ K 11 K 12 .

(D.13)
Combining this result with the result in appendix B, we can express the square of the structure constant simply as a ratio of two determinants as where G u|φv,w is the induced Gaudin determinant defined in (B.11). Let us finally mention the potential application to the mirror corrections. The structure of the interaction term in the mirror-particle integrand given in [5] takes the same form as (D.3). Therefore, one can use the generalized Schur's Pfaffian identity also for the mirror particles and recast each term as a Pfaffian. Furthermore, in the case of mirror particles, the full expansion coincides exactly with the expansion of the Fredholm Pfaffian. It is still to be seen if the sum over bound-state indices leads to a further simplification, but in any case, such an expression would certainly be useful for resumming the finite size correction [55] at finite coupling.

E The so(6) structure constant at tree level
In this appendix we compute the tree-level three point function of operators in the so(6) sector of planar N = 4 SYM. We focus on the case of two 1 2 -BPS operators and one non-BPS single-trace operators: The protected operators are given by the BMN vacua spanned by the elementary field Z and by the rotated fieldZ ≡ Z +Z + X −X . The non-protected operator is given by an eigenstate of the one-loop dilatation operator. Such operator can be obtained by diagonalizing the Hamiltonian of the dual integrable spin chain [56], using the Algebraic Bethe Ansatz (ABA). We develop on such construction in section E.3.

JHEP07(2019)082
"l" bridge "L-l" bridge trivial bridge This three point function can be computed by Wick contractions as shown in figure 8. Following the tailoring procedure, introduced for the su(2) sector in [57], we can express the wick contractions as scalar products of spin chain states dual to the single trace operators. In our configuration -dubbed the reservoir picture in [2] -we have two trivial bridges which only feature propagators of the type Z-Z (blue lines). The only non-trivial wick contraction comes therefore from the bridge l = (L + L 1 − L 2 )/2 between operators O 1 and O 3 and is given by the spin chain scalar product 25 where Ψ l is a sub-chain of length l in the cyclic spin chain state Ψ, dual to the single trace operator i.e. O 3 ≡ |Ψ . In this way the computation of the tree level three-point function is reduced to finding an inner product of states in the dual so(6) spin chain. In the su(2) sector the relevant scalar product was computed with ABA techniques [57], obtaining sum formulas or more compact determinant expressions for some special cases, see [58] for a review. In [59] su(3) scalar products were computed and used in [60] for the computation of structure constants in this sector. While for the so(6) spin chain there are no available formulas for the scalar products 26 in the literature. One of the obstacles being the complexity of the ABA for so(2n) models [62]. To overcome this problem we 25 This renders an unnormalized structure constant. The normalized version includes the norms of the three operators involved. 26 An attempt to conjecture the so(6) scalar product, based on the results for su (2), was given in [61].

JHEP07(2019)082
construct an alternative version 27 of the so(6) ABA which allows for a simpler approach to the computation of scalar products. 28 In particular we use this machinery to find the scalar product (E.2) as a sum formula. This result makes direct contact with the conjecture (3.13), restricted to so(6) at tree level, showing the structure of a sum over partitions and a matrix part depending on the nested levels of the Bethe Ansatz. It would be interesting to generalize this Bethe Ansatz to address operators in the full sector psu(2, 2|4) and reproduce the conjecture (3.13) at tree level. The rest of this appendix is organized as follows: in section E.1 we introduce the so(6) integrable spin chain, the correspondent transfer matrix, as well as some notation. In section E.2 we present a novel so(6) vertex model, specifying the Boltzmann weights and the way to extract the Bethe states from the lattice. In section E.3 we develop a so(6) ABA for the diagonalization of the spin chain Hamiltonian and transfer matrix. In section E.4 we present the Yang-Baxter algebra, showing how to use it to derive the so called wanted and unwanted terms of the ABA. In section E.5 we present the coordinate Bethe Ansatz (CBA) which can be derived from our ABA and vertex model. Finally, in section E.6 we put in used the Yang-Baxter algebra to compute the tree level structure constant given by the scalar product (E.2). We show how to simplify the result to obtain the so(6) tree-level analog of the conjecture (3.13).

E.1 The so(6) spin chain
To obtain a basis of non-BPS operators we need to diagonalize the so(6) integrable spin chain Hamiltonian: where I,P and K are identity, permutation and trace operators respectively. This spin chain Hamiltonian is proportional to the one loop dilatation operator in the so(6) sector of N = 4 SYM. The basis of eigenstates of this Hamiltonian constitutes a basis for non-BPS operator of the one loop so(6) sector, which (partially) lifts the original tree level degeneracy. The single trace operators of this sector are mapped to cyclic states of the spin chain as: where the elementary fields are given by the complex scalars fields: These scalar fields form a multiplet of the antisymmetric 6 representation of su(4), isomorphic to the vector representation of so (6). This isomorphism is realized by the transforma-JHEP07(2019)082 tion.
In this appendix we stick to the basis of complex scalars (E.5).
The ABA finds the spectrum of the so(6) Hamiltonian by solving the eigenvalue problem of the transfer matrix, the trace of a monodromy operator. This operator is constructed by "scattering" a probe particle, in an auxiliary space Λ and spectral parameter (momentum)u, with all the spin chain sites. We can build various monodromies by choosing the auxiliary space to lie in any of the representations of the spin chain symmetry group: From all these possibilities one is distinguished and corresponds to the choice of auxiliary space in the same representation as the spin chain sites. The trace of this special choice, the transfer matrix, is a generating function of a family of local conserved charges including the nearest-neighbour Hamiltonian of the spin chain. For our so(6) spin chain the distinguished monodromy is T 6 , its trace generates local conserved charges such as the Hamiltonian (E.3). Other choices of auxiliary space do not generate the spin chain Hamiltonian, nevertheless their correspondent transfer matrices are in convolution with the distinguished one. This means that we can address the eigenvalue problem for the spin chain Hamiltonian and all the transfer matrices in convolution at once. So we can choose to solve the eigenvalue problem of the simplest transfer matrix. For our so(6) spin chain the simplest choice corresponds to T 4 with auxiliary space in the 4 fundamental representation of su(4) as: where Φ a k b k and Φ c k d k are the incoming and outgoing so(6) flavours of the k th spin chain site in the "scattering" with the auxiliary particle. The indexes a and b indicate the JHEP07(2019)082  with the nodes of the so(6) or su(4) Dynkin diagram: To each of the simple lines we associate a probe particle with rapidities: v for the vertical lines on the left wing , u for the horizontal lines and w for the vertical lines on the right wing of figure 10. These probe particles "scatter" among each other forming wing auxiliary lattices. This scattering or crossing of lines is controlled by the R-matrix: which also provides the Boltzmann weights of the vertex model in the wing lattices as: This R-matrix is composed of two terms: the identity I that keeps the flavours in their original direction and the permutation P that swaps the directions the flavours follow. On the right hand side of equation (E.10) the arrows indicate the flow of the flavour, while the auxiliary root u(v) is always attached to the horizontal(vertical) direction.
The novel part of our vertex model is given by the core of the lattice, where simple (4) and double (6) lines crossed. The Boltzmann weights associated to these crossed lines are given by the R 46 -matrix: Now the permutation term P 46 has two pieces corresponding to the two possibilities of swapping the flavour index a with the indexes of Φ bc .

E.2.2 The boundary conditions and Bethe state
In order to obtain the so(6) Bethe state from this vertex model we fix the boundaries of the wings as shown in figure 10. With these restrictions the auxiliary spaces on the wings become effectively two-dimensional and the wing lattices render two 6-vertex models: wing 3,4 and wing 1,2 , associated to su (2)  Following the rules of this vertex model it is possible to determine the general form of the Bethe state as a linear combination of the coordinate basis. We present this in section E.5 as a Coordinate Bethe Ansatz(CBA). In the following section we present the origin of this vertex model from the ABA. The "scattering" of a probe particle in the 4 representation with the spin chain sites in the 6 representation is given by a product of R 46 matrices and renders the T 4 monodromy JHEP07(2019)082 matrix: 6 (E.14) From the point of view of the auxiliary space the monodromy is a 4 × 4 matrix, whose elements are operators that act exclusively over the spin chain: 15) In order to obtain the elements of this matrix in the graphical representation in (E.14), we simply fix the boundaries of the horizontal line to take specific flavour values {1, 2, 3, 4}. The correspondent transfer matrix, given by the trace of the monodromy (E.15), is: Now our aim is to construct the ABA to find the eigenvalues and eigenstates of this transfer matrix. For this we start by identifying one of its trivial eigenstates, the pseudo-vacuum: We consider this pseudo-vacuum as a reference state to start the construction of the ABA. With this choice the action of the monodromy matrix organizes into 2 × 2 blocks. The pseudo-vacuum diagonalizes the A and D blocks: with (considering θ k = 0): and it is annihilated by the C-block elements:

JHEP07(2019)082
The action of the B-operators over (E.17) is much less trivial. They create magnon-states or plane waves when acting over the pseudo-vacuum. The operator B jk injects flavour k ∈ {3, 4} and absorbs flavour j ∈ {1, 2} from the state it acts over. When acting over the pseudovacuum (E.17) it creates a magnon of type Φ 1k when j = 2 or type Φ 2k when j = 1. For instance the operator B 13 creates a X ≡ Φ 23 magnon as: u−θn−i/2 as can be read off from the lattice with Boltzmann weights (E.11).
Similarly we can create other magnon-states with different flavours or so(6) charges by using operators B 23 , B 24 and B 14 . The relationship between these creation operators and the so(6) charges is summarized in the following figure:

E.3.2 The Bethe ansatz
The B creation operators play a key role in the construction of the spectrum of the transfer matrix. The states created by repeated action of B-operators over the reference state (E.17) serve as a basis to propose a general Ansatz for the eigenstates of the transfer matrix as: where the indexes a k andã k take on flavour values {3, 4} and {1, 2} respectively, and the unconstrained tensors ψ a 1 ···a M andψã 1 ···ȧ M weight the contributions of states constructed by different choices of Bã mam -operators. We can rewrite the ansatz (E.23) by using 2 × 2 B-blocks instead of individual Bã k a koperators. With this purpose we introduce the wing-auxiliary states |ψ and ψ | as: ψ a 1 a 2 ···a M = a 1 a 2 · · · a M |ψ (E.24) ψã 1ã2 ···ã M = ψ |ã 1ã2 · · ·ã M (E.25)
Using B-blocks and wing-auxiliary states we reformulate the Bethe Ansatz (E.23) as: Using the Ansatz (E.28) we now need to solve the eigenvalue problem: This means we need to find the restrictions over the auxiliary roots {u} and the wings states such equation (E. 29) holds. In what follows we sketch the steps to achieve this diagonalization. These will heavily rely on the Yang-Baxter algebra presented in section E.4. We first reexpress the transfer matrix (E.16) by defining block operators A and D as: with non-trivial action over the spaces: A a :2 a ⊗ 6 1 ⊗ · · · ⊗ 6 L −→2 a ⊗ 6 1 ⊗ · · · ⊗ 6 L (E.31) D a : 2 a ⊗ 6 1 ⊗ · · · ⊗ 6 L −→ 2 a ⊗ 6 1 ⊗ · · · ⊗ 6 L (E. 32) In this language the transfer matrix is now a sum of traces of A and D blocks and the eigenvalue problem has two pieces associated to these blocks: where "a" labels an auxiliary space2 for the A-block and 2 for the D-block. Now starting with equation (E.33) the strategy is to commute the A and D blocks through the product of B-blocks until we reach the pseudo-vacuum that satisfies (E.18). This is possible using the commutation relations provided by the Yang-Baxter algebra (see section E.4). Once we follow this procedure the result has two type of terms: wanted and unwanted. From the wanted terms we can reproduce the eigenvalue equation (E.29) and read off the correspondent transfer matrix eigenvalue Λ. On the other hand the unwanted terms spoil the eigenvalue equation and their vanishing is a necessary condition to satisfy (E.29). Here we only show the wanted terms: with: As a by-product of the commutations of A-B and D-B blocks in (E.34) we obtain two auxiliary nested su(2) monodromies 29 acting in the spaces of ψ | and |ψ respectively: with R-matrices given by (E.9) but now restricted to act over su(2) subspaces. Namely they read Furthermore we can directly act with A over the pseudovacuum in the first line of (E.34), since it does not act non-trivially over the wing states. Similarly in the second line we can commute D a and the nested monodromy T a in the presence of the trace and act over the pseudovacuum. Using the diagonalization properties (E.18) of pseudovacuum we can simplify (E.34) and obtain wing transfer matrices as traces of (E.36) and (E.37): We now see that to reproduce the eigenvalue equation (E.29) from the "wanted" terms, the wing states must be eigenstates of the corresponding nested su(2) transfer matrices of the monodromies (E.36) and (E.37). This auxiliary problem is simply solved by the standard su(2) ABA [64]: ψ | = Ω2| C(w 1 ) · · · C(wK) and |ψ = B(v 1 ) · · · B(v K )|Ω 2 (E. 40) where B and C are creation and annihilation operators extracted from the monodromies (E.37) and (E.36) respectively. They act over su(2) vacuum states of the wings given by: In addition the sets of auxiliary roots {v} and {w} must be on-shell, that is they must fulfil su(2) Bethe equations with the set {u} as inhomogeneities. Assuming these conditions the su(2) Bethe states (E.40) diagonalize the wing transfer matrices as: Figure 11. The so(6) Bethe state.
With the wing-states on-shell , (E.39) becomes the eigenvalue equation (E.29) up to unwanted terms: More explicitly the transfer matrix eigenvalue is given in terms of the spectral parameter λ and the sets of auxiliary roots {u} M , {v} K and {w}K: The vanishing of the unwanted terms puts constrains over the set of roots {u}. These constrains constitute the Bethe equations of the so(6) middle node. Alternatively we can arrive to the same conditions by imposing the vanishing of the spurious poles at λ = u 1···M of the transfer matrix eigenvalue (E.44). This latter method to obtain Bethe equations is the so called analytic Bethe Ansatz: 45) To obtain the standard form of so(6) Bethe equations we must perform the shifts: 30 In summary, the so (

E.4 The Yang-Baxter algebra
In this section we present the Yang-Baxter algebra of our so(6) model. This algebra provides the technical steps to find the wanted and unwanted terms of the Bethe Ansatz in section E.3, as well as for our final result (E.83) for the scalar product (E.2). The R matrices in (E.9) and (E.11) fulfil the Yang-Baxter equation: where a and b label two different spaces in the 4 fundamental representation. This can be straightforwardly generalized to a Yang-Baxter relation for the monodromy in (E.14), the so called RT T relation: 31 Furthermore, taking the trace over the 4-dimensional auxiliary spaces in (E.48) we obtain the commutation relation for the transfer matrices: The RT T relation also provides the algebra of the monodromy elements in (E.15), known as the Yang-Baxter algebra. This is a set of commutation relations that can be obtained by specifying the boundary conditions in the four-dimensional auxiliary spaces: a : (i) → (k) and b : (j) → (l) as follows: where the lower indexes in parenthesis indicate the initial flavours and the upper indexes correspond to the final flavours. We leave implicit the intermediate flavours indexes over which we must sum over.
Expanding the R-matrices in (E.50) into identity and permutation as in (E.9) we obtain the following algebra of operators: In our graphical representations the order of action of operators should be read from left to right. While in our equations we respect the usual order of operator action, that is operators on the right act first. In this way the left-hand side in figure 12 represents the right-hand side on equation (E.48).

JHEP07(2019)082
where T 13 ≡ B 13 and likewise for other operators in (E.15). The Yang-Baxter algebra (E.51) plays a key role in the construction of the ABA in section E.3 and also in the computation of the scalar product that gives the tree level structure constant in section E. 6. In what follows we provide some of the details involved in these calculations.

E.4.1 The wanted and unwanted terms of the Bethe ansatz
In the ABA construction we use the Yang-Baxter algebra organized into blocks. For this we restrict the four-dimensional auxiliary spaces in (E.48) to the subspaces 2 (flavours {3, 4}) or2 (flavours {1, 2}), instead of strictly fixing the boundary conditions. For instance to derive the D-block and B-block commutation relation we restrict the auxiliary spaces as: a : 2 (3,4) → 2 (3,4) and b : 2 (3,4) →2 (1,2) : (3,4) (3,4) (E.52) expanding the R-matrix of the left hand side we obtain: Under these restrictions now a and b label two-dimensional spaces. Similarly, by making another appropriate choice of boundary conditions: a :2 (1,2) → 2 (1,2) and b : 2 (3,4) →2 (1,2) , we obtain the A-B commutation relation: Using these commutation relations, (E.53) and (E.54), we can compute the wanted and unwanted terms as a result of commuting A and D through the product of B-blocks in the ansatz (E.28). These results are given by: where we use the short-hand notation: as well as:

JHEP07(2019)082
The Φ k coefficients are: Similarly we commute a D-block through a product of B-blocks as: with coefficients:

E.4.2 The C-commutation relations for the scalar product
When computing the scalar product (E.2), using the ABA, we will need of the commutation relations of the annihilation C-block and the other elements of the monodromy. These can be obtained from the Yang-Baxter algebra (E.51) as: where |Z X · · · stands for an element of the coordinate basis and the coefficient Ψ ZX··· is its correspondent wave-function. This wavefunction is obtained as a partition function in the lattice in figure 10, with top boundary conditions imposed by the correspondent element of the coordinate basis.

E.5.1 The coordinate basis as strings of auxiliary roots
To introduce the wave-function of a given element in the coordinate basis we first define a one to one map between the elementary fields (E.5) and a set of rapidities:

JHEP07(2019)082
where θ's are the inhomogeneities defined for each spin chain site. To make manifest the structure of the nested Bethe ansatz, here we represented the fields by stacking the roots: u is the root at the middle node whereas w and v are the nested roots associated with the left and the right nodes respectively.
Using this representation, we can re-express the coordinate basis as a collection of (sets of) rapidities, for instance: Here we assigned a numeration to the auxiliary roots in the order of appearance. 32 In what follows, we call such a collection of rapidities a string.

E.5.2 The wave-function
The wavefunction Ψ s for the string s is given by a sum over weighted permutations over all the auxiliary roots: (E.68) where the notation { * } σ denotes that the set * is permuted according to the permutation σ. The multiparticle S-matrix S({u} σ ) brings the ordered momenta {u} to the ordering {u} σ and is given by a factorized product of two-body S-matrices as in the examples: In a spin chain with L sites, the correspondent "bare" wavefunction Ψ bare s is given by

JHEP07(2019)082
where the individual wave function Φ(s n ), defined for the n-th element of the string s, is given by where we include the label of the rapidities (u m ) only when this is necessary to express the correspondent wavefunction. The factors ϕ andφ are one-particle wave functions and are given by: Needless to say, when the rapidities are permuted in (E.70), we should also permute the rapidities in the definitions of the wave functions, which are given by the right hand sides of (E.71), accordingly.
E.6 The scalar product: tree level structure constant As we saw in the introduction of this appendix, the tree level planar three-point function in figure 8 is given by the scalar product between a rotated BMN vacuum and a Bethe state: where Bethe state is given by the Ansatz in figure 11 but with the number of sites or elementary fields equal to the bridge "l".

E.6.1 A global rotation
In order to compute the scalar product (E.73) using the machinery of the ABA, we first need to express the rotated vacuum in this language. This is achieved by means of a global rotation of the original vacuum:

JHEP07(2019)082
The generator "b" of this global rotation can be simply obtained from the B-block by taking the trace and sending to ∞ the spectral parameter "u"; this lowering generator is composed by the elements: When b 13 acts over the vacuum generates a X excitation, when b 24 acts generates -X and when both act over the same site aZ excitation is generated. In this way we generate the rotated vacuum:Z In the scalar product (E.73) we must use instead the bra state for which we use the C-block as:

E.6.2 The scalar product
Now we outline the steps we take to compute the scalar product (E.73). The first step is to notice that the Bethe state has a defined so(6) charge determined by the number of (finite) auxiliary roots {u} M ,{v} and {w}. While in the expansion of global rotation e c = 1 + c + · · · only the term c M matches this so(6) charge . So the scalar product (E.73) can be simplify to: The next step is to commute all the c operators through the B-blocks, such as we can use their annihilation properties: c |Ω = 0 and Ω| B(u) = 0 (E.80) The commutator of c and B can be found from the Yang-Baxter algebra. Taking the limit (E.78) of (E.61) we obtain: Since A and D blocks are generated we also need of their commutators with c, which can be extracted in a similar way from (E.62) and (E.63) as:  All in all the result of commuting c M through a product of M B-block operators is given by a sum over bi-partitions: where we use the short-hand notation: and the C and c-terms are products of operators that annihilate the so(6) pseudo-vacuum. The matrix operator R α,ᾱ is a product of su(2) R-matrices (E.38) that changes the order of the roots {u 1 · · · u M } to the order {α,ᾱ}, while the operator Rᾱ ,α takes the roots from the ordering {ᾱ, α} to the ordering {u 1 · · · u M }.
For instance, when α = {u 2 , u 4 } andᾱ = {u 1 , u 3 , u 5 }, the multiparticle scattering operators are: where we use the short-hand notation R ab = R ab (u a , u b ) Finally we can compute the scalar product (E.79) using (E.83) and the eigenstate relations of the pseudovacuum (E.18): The matrix part can be further simplified considering that the wing-states are (on-shell) su(2) Bethe states. In this case the action of the R-matrices has the simple effect of reshuffling the inhomogeneities {u} of the wing-states so we obtain: where |ψ α,ᾱ is the wing state |ψ with the inhomogeneities order as in the top of figure 13a and ψᾱ ,α | is the wing ψ | with the ordering as in the bottom of figure 13b. In this way the matrix part is simply given by the scalar product of the states (E.87). To simplify this scalar product it is necessary to place the inhomogeneities of the two states in the same ordering. This can be achieved by using again the multi-scattering R matrix: The special feature of this reordering scattering matrix is that it can be expressed as a product of nested su(2) transfer matrices with spectral parameters u ∈ᾱ and inhomogeneities {u}: Since our wing states are on-shell we can replace the transfer matrices by the correspondent eigenvalue and obtain the periodicity relation: So for on-shell Bethe states the cost of reordering is just a phase. Then the matrix part is given by: We now have the scalar product of two on-shell Bethe states with inhomogeneities in the same ordering. Considering their orthogonality property we know this scalar product vanishes unless the set of wing roots are identical {v} = {w}. 33 Under this condition the scalar product is given by the Gaudin-determinant. This determinant is invariant under permutations of the inhomogeneities so it can be taken out of the sum over partitions. The final expression for the unnormalized tree level structure constant is: where Λ su(2) (u) = K k=1 u−v k +i/2 u−v k −i/2 , after performing the shift (E.46). 33 Here we refer to the shifted wing roots v → v + i/2 and w → w − i/2 which appear in the standard form of the so(6) Bethe equations.

JHEP07(2019)082
F Nested hexagons in the su(1, 1|2) sector In this appendix, we study the hexagon amplitudes in the su(1, 1|2) sector and provide a detailed derivation of the formulas given in section 3. This is a large sector of operators, the largest closed subsector containing both the sl(2) and su(2) diagonal sectors, see e.g. [16]. Still, the analysis in this sector turns out to be remarkably simple. In the following, we review the Bethe ansatz wave functions for these operators [65] and employ them to the computation of the hexagon amplitudes.
Wave functions. The states of interest are usual BMN operators above the BPS vacuum |0 = tr Z L , with each excitation χ lying in a (1|1) ⊗ (1|1) irrep of the residual symmetry subalgebra su(1|1) ⊕ṡu(1|1). Equivalently, a magnon χ can take 4 possible values, out of the 16 ones, with Y a complex scalar, Ψ andΨ two gauginos, and D a lightcone covariant derivative. The dynamics factorizes along the two wings and a generic scattering eigenstate can be written as a tensor product of a left and a right wave function, built out of φ|ψ andφ|ψ respectively. We recall below the structure of these wave functions, focusing on the left part of the state. The right part defines an isomorphic problem and all the formulae hereafter apply to it after "dotting" whatever can be dotted. (One must useu = u, for the main roots, since these are inhomogeneities common to the two factorized auxiliary spin chains.) Off-shell Bethe states for the inhomogeneous su(1|1) spin chain are constructed in the usual manner, using the restriction of the full su(2|2) S-matrix [19] to the φ|ψ subspace as a fundamental R-matrix. The analysis is quite similar to the standard algebraic Bethe ansatz for the su(2) spin chain with the difference that the B operator here is fermionic. As a starting point, one must choose a pseudovacuum state. We shall work in the su(2) grading for which the reference state, the so-called level II vacuum, is chosen to be made of scalars only The subscript u reminds us that this vacuum state depends on the ordering of the lattice. Acting with a B operator produces a fermionic "spin wave" along the chain 34 The B(y) operator is obtained by scattering a probe particle with rapidity y = x − 0 through the chain, with the boundary condition that it starts as a fermion and ends as a boson. When defining B in this way, using the S-matrix of [19], we also strip out an inessential overall factor. The latter factor is invariant under permutation of the spin chain inhomogeneities and thus does not affect the fundamental property of the B operator.

JHEP07(2019)082
with the wave function [65] Ψ n (y) = a n y − x + n n−1 j=1 S II,I (y, x j ) . (F.5) Here a n = i(x − n − x + n ) is a free parameter of the representation (for the relative normalization between boson and fermion), which we have fixed to its unitary value, for convenience. The phase can be seen as the S-matrix for bringing the fermion through the lattice site x j . The scattering among the fermionic waves happens to be trivial S II,II (y 1 , y 2 ) = 1, up to the statistics. Hence the multiparticle wave function is the totally antisymmetric product of the individual spin waves, Note that all these wave functions fulfill, as a consequence of the Yang-Baxter equation, the so-called compatibility condition [65] S π |Ψ = A π |Ψ π , (F. 8) where π is an arbitrary permutation of the labels u and |Ψ π = B π(u) (y 1 ) . . . B π(u) (y F )|0 II π(u) is the state obtained through this relabelling. The overall factor in (F.8) is the eigenvalue of the scattering matrix S π on the spin chain vacuum, with the product running over the pairwise scattering events in the permutation π and with A 12 the scattering phase for two identical bosons, see [19]. 35 Lastly, imposing periodic boundary conditions for the spin waves generates the Bethe ansatz equations of the level II rapidities y = {y j }. Looking back at the wave function (F.5) and recalling that fermions do not interact in this model, one easily infers that this choice corresponds to 36

JHEP07(2019)082
Hexagon amplitudes. Equipped with the wave functions, we can attack the problem of evaluating the hexagon amplitudes in the su(1, 1|2) subsector. We shall focus on configurations with an equal number of left and right fermions on each hexagon, since these are the sole configurations for which the form factors are nonzero. Symmetrywise a left fermion adds 1/2 to the left Lorentz spin of the operator and similarly for a right fermion. Hence the overall condition that F =Ḟ is just saying that operators in the OPE of two scalars are in symmetric (traceless) Lorentz representations.
We begin with the simplest set-up where all the magnons are sitting on the same hexagon, i.e., α = u. In this case there is only one hexagon form factor to compute, the one for the full state, where (−1) f is a grading factor for bringing the original state (F.1) to the above factorized form, with all the left magnons standing on the left of the right magnons. When F =Ḟ , we can just write f = F (F − 1)/2. The rule for evaluating the hexagon form factor h Ψ |Ψ is to first scatter the left part from its ingoing to its outgoing configuration and then contract the outcome with the right part using the left-right inner product (χ K . . . χ 1 ,χ 1 . . .χ K ) = K j=1 δ χ jχj , (F.12) where δ χχ = 1 if the two excitations are the same, and δ χχ = 0 otherwise. The latter orthogonality condition between bosons and fermions, together with the fact that the number of fermions is preserved by the S-matrix (in the su(1, 1|2) subsector), imply that the amplitude (F.11) is zero if the numbers of left and right fermions are different. As alluded before, this is a consequence of the diagonal symmetry preserved by the hexagon, and the argument applies to each partition separately in a generic split configurations with α,ᾱ = ∅.
To avoid possible confusion, we shall use round brackets (., .), as in (F.12), to denote the left-right overlap. This one is a real bilinear form on the tensor product of the left and right Hilbert spaces H ⊗Ḣ. It is different from the usual Hermitian inner product .|. , defined separately on each Hilbert space, which is of course sesquilinear. Since the left and right Hilbert spaces are isomorphic we can also overlap left and right states using the Hermitian inner product. The relation to the round product is then given by Ψ|Ψ = (Ψ † ,Ψ) , (F. 13) where Ψ † is the Hermitian conjugate of the state Ψ, described by the complex conjugate wave function on the transposed lattice.
Introducing the permutation π that reverses the ordering of the spin chain, we can write the hexagon amplitude as h Ψ |Ψ = (S π Ψ,Ψ) . (F.14)

(F.18)
It follows from it that we can write the hexagon form factor as the overlap with the last equality holding for a real state, that is such that ← − Ψ n (y) = Ψ n (y) * . Notice that the prefactor S Ψ has a simple interpretation. If we think of our state as being made of the union of the y-and u-"diagonalized excitations", then S Ψ is the S-matrix obtained upon reversing the ordered set y ∪ u. It is not really significant for the computation of the structure constant however. The reason is that the hexagon amplitude only determines the structure constant up to an overall factor that implements the change of normalization between infinite-and finite-volume Bethe states. 37 This one contains, in particular, the factor 1/ S Ψ SΨ that removes the dependence on the ordering of the rapidities and cancels (F.18) when Ψ =Ψ. The grading factor (−1) f also cancels out in (F.11). Hence, in the end, we are left with the product ( ← − Ψ |Ψ) for the partition (α,ᾱ) = (u, ∅), in line with our main formula. The proof of (F.17) is immediate. Permuting the u labels in the wave function (F.5), with π(i) = K + 1 − i, yields In the general case, we have to partition the chain into two non-empty subsets, α ∪ᾱ = u. For simplicity, but with no loss of generality, we consider a partition that preserves the JHEP07(2019)082 ordering of the full set u, that is, α = {x 1 , . . . , x l } andᾱ = {x l+1 , . . . , x K }, with l = |α|. The action of the S-matrix S π = S πα S πᾱ now depends on whether the fermion is on the α or on theᾱ partition, since π α reverses the order of α only, and similarly for πᾱ. If 1 n l we find the same result as in (F.20), Ψ π n (y) = l j=1 S II,I (y, x j ) × ← − Ψ π(n) (y) , (F. 21) with the product being restricted to the x's in the subset α, but if n > l we get an extra phase in front, Ψ π n (y) = Remarkably, when the roots y are on shell, this difference becomes immaterial and the flipped wave function Ψ π is found to be proportional to the outgoing one ← − Ψ again. After taking care of the way the determinant changes under the reordering of the y's, we get for a multi-fermion state on the split configuration. Here, A πα = i<j∈α A ij and f α = F α (F α − 1)/2, with F α the number of fermions on the subchain α, and similarly forᾱ. The grading factors disappear in the full amplitude, giving A αᾱ = A πα A πᾱ × S I,II (ᾱ, y) × ( ← − Ψ |Ψ) , (F.24) after using the BAEs (F.10) and the unitarity of the S-matrix, S I,II (y, x) = 1/S II,I (x, y), to rewrite the y-dependent factor. Formula (F.24) is the su(2) counterpart of the one given in section 3.2 in the sl(2) grading. They both convey the same message, namely, that the amplitude for the split configuration is proportional to the inner product Ψ|Ψ up to the phase S I,II (ᾱ, y) for the scattering of the roots inᾱ with the higher level rapidities y. The overall A factors in (F.24) accounts for the difference of gradings between (F.24) and (3.14). They can be absorbed in the dynamical parts of the hexagon form factors using the relation between pure sl(2) and pure su(2) hexagon amplitudes [2] h ij | su(2) = A ij h ij | sl (2) . where in the last equality we used that the sum telescopes. Hence, when the two rapidities are on-shell the term in bracket vanishes, see (F.10), and so does M ij , ifẏ i = y j . When the two states are the same, Ψ =Ψ (or y =ẏ, with the rapidities ordered in the same way), the matrix M is diagonal and the norm of the state is the product of the individual norms, We can use the freedom we have to rescale the individual wave functions, (F.5) and (F.16), by some function of y to convert ∂ y into ∂ v in (F.30), with v = g(y + 1/y), and match the convention used in the bulk of the paper, see (6). This choice is irrelevant for the normalized structure constant, which is divided by the square root of the Gaudin determinant G for the full wave function (which includes the level I and all the lower ones). 38 Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited. 38 Implicit here is the assumption that the norm of the level II wave function (F.29) matches with a minor of G, as obtained by freezing the level I rapidities u.